home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 471 / rccl113 < prev    next >
Text File  |  1987-03-02  |  1KB  |  54 lines
  1. # double cbrt(arg)
  2. # double arg
  3. # no error exits
  4. #method: range reduction to [1/8,1], poly appox, newtons method
  5. # J F Jarvis, August 10,1978
  6. .globl    _cbrt
  7. .text
  8. .align    1
  9. _cbrt:
  10.     .word    0x00c0
  11.     clrl    r3
  12.     movd    4(ap),r4
  13.     jgtr    range
  14.     jeql    retz
  15.     mnegd    r4,r4    # |arg| in r0,r1
  16.     movl    $0x100,r3    # sign bit of result
  17. range:
  18.     extzv    $7,$8,r4,r6
  19.     insv    $128,$7,$8,r4    # 0.5<= frac: r4,r5 <1.0
  20.     clrl    r7
  21.     ediv    $3,r6,r6,r7    # r6= expnt/3; r7= expnt%3
  22.     addb2    $86,r6
  23.     bisl2    r3,r6    # sign,exponent of result
  24.     polyf    r4,$3,pcoef    # initial estimate is Hart&Cheney CBRT 0642
  25.                         # D=4.1
  26.     muld3    r0,r0,r2    # Newtons method, iteration 1, H&C 6.1.10
  27.     divd3    r2,r4,r2
  28.     subd3    r2,r0,r2
  29.     muld2    third,r2
  30.     subd2    r2,r0    # D=8.2
  31.     muld3    r0,r0,r2    # iteration 2
  32.     divd3    r2,r4,r2
  33.     subd3    r2,r0,r2
  34.     muld2    third,r2
  35.     subd2    r2,r0
  36.     muld2    hc[r7],r0    # set range
  37.     insv    r6,$7,$9,r0    # set sign,exponent
  38.     ret
  39. retz:
  40.     clrd    r0
  41.     ret
  42. .data
  43. .align    2
  44. third: .double 0d0.33333333333333333333e+0
  45. hc:
  46.     .double 0d1.25992104989487316476e+0
  47.     .double 0d1.58740105196819947475e+0
  48.     .double    0d1.0e+0
  49. pcoef:
  50.     .float 0f0.1467073818e+0
  51.     .float 0f-0.5173964673e+0
  52.     .float 0f0.9319858515e+0
  53.     .float 0f0.4387762363e+0
  54.