home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 471 / rccl106 < prev    next >
Text File  |  1987-03-02  |  1KB  |  47 lines
  1. # double sqrt(arg)
  2. # double arg
  3. # if(arg<0.0) { _errno=EDOM; return(0.0) }
  4. # J. F. Jarvis August 2, 1978
  5. .set    EDOM,98
  6. .text
  7. .align    1
  8. .globl    _sqrt
  9. .globl    _errno
  10. _sqrt:
  11.     .word    0x0c0
  12.     movd    4(ap),r4
  13.     jgtr    range
  14.     jeql    retz
  15.     movl    $EDOM,_errno
  16. retz:
  17.     clrd    r0
  18.     ret
  19. #
  20. range:
  21.     extzv    $7,$8,r4,r6    # r6 is exponent of arg
  22.     insv    $128,$7,$8,r4    # r2: 0.5<=fraction(arg)<1.0
  23.     incl    r6
  24.     clrl    r7
  25.     ediv    $2,r6,r6,r7    # r6=(exp+1)/2; r7=(exp+1)%2
  26.     addb2    $64,r6    # r6 is correct exponent for result
  27.     polyf    r4,$4,pcoef    # init estimate of sqrt(frac)
  28.                         # Hart&Cheney SQRT 0132 D=5.1
  29.     divd3    r0,r4,r2    # Newtons method, 2 iterations
  30.     addd2    r2,r0
  31.     muld2    $0d0.5e+0,r0
  32.     divd3    r0,r4,r2    # Hart&Cheney 6.1.7
  33.     addd2    r2,r0    #d=21 at exit
  34.     muld2    hc[r7],r0    # *sqrt(2) requ for even org exp.
  35.     insv    r6,$7,$8,r0    # insert correct exp.
  36.     ret
  37. .data
  38. .align    2
  39. pcoef:
  40.     .float 0f-0.1214683825e+0
  41.     .float 0f0.5010420763e+0
  42.     .float 0f-0.9093210498e+0
  43.     .float 0f0.1300669049e+1
  44.     .float 0f0.2290699453e+0
  45. hc:        .double 0d0.35355339059327376220e+0    # sqrt(2)/4
  46.         .double 0d0.5e+0
  47.