home *** CD-ROM | disk | FTP | other *** search
/ Between Heaven & Hell 2 / BetweenHeavenHell.cdr / 500 / 471 / rccl103 < prev    next >
Text File  |  1987-03-02  |  1KB  |  69 lines
  1. # double sqrt(arg):revised July 18,1980
  2. # double arg
  3. # if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) }
  4. # W. Kahan's magic sqrt
  5. # coded by Heidi Stettner
  6.  
  7.     .set    EDOM,98
  8.     .text
  9.     .align    1
  10.     .globl    _sqrt
  11.     .globl    dsqrt_r5
  12.     .globl    _errno
  13. _sqrt:
  14.     .word    0x003c          # save r5,r4,r3,r2
  15.     movd    4(ap),r0
  16.     bsbb    dsqrt_r5
  17.     ret
  18. dsqrt_r5:
  19.     movd    r0,r4
  20.     jleq    nonpos        #argument is not positive
  21.     movzwl    r4,r2
  22.     ashl    $-1,r2,r0
  23.     addw2    $0x203c,r0    #r0 has magic initial appx
  24.  
  25. # Do two steps of Heron's rule
  26.  
  27.     divf3    r0,r4,r2    
  28.     addf2    r2,r0
  29.     subw2    $0x80,r0
  30.  
  31.     divf3    r0,r4,r2
  32.     addf2    r2,r0
  33.     subw2    $0x80,r0
  34.  
  35.  
  36. # Scale argument and approximation to prevent over/underflow
  37. # NOTE: The following four steps would not be necessary if underflow
  38. #       were gentle.
  39.  
  40.     bicw3    $0xffff807f,r4,r1
  41.     subw2    $0x4080,r1        # r1 contains scaling factor
  42.     subw2    r1,r4
  43.     movl    r0,r2
  44.     subw2    r1,r2
  45.  
  46. # Cubic step
  47.  
  48.     clrl    r1
  49.     clrl    r3
  50.     muld2    r0,r2
  51.     subd2    r2,r4
  52.     addw2    $0x100,r2
  53.     addd2    r4,r2
  54.     muld2    r0,r4
  55.     divd2    r2,r4
  56.     addw2    $0x80,r4
  57.     addd2    r4,r0
  58.     rsb
  59. nonpos:
  60.     jneq    negarg
  61.     clrd    r0        #argument is zero
  62.     rsb
  63. negarg:
  64.     movl    $EDOM,_errno
  65.     mnegd    r4,-(sp)
  66.     calls    $2,_sqrt
  67.     mnegd    r0,r0        # returns -sqrt(-arg)
  68.     ret
  69.