4.3BSD - /usr/src/usr.lib/liboldnm/sqrt.s

   1: #
   2: # Copyright (c) 1980 Regents of the University of California.
   3: # All rights reserved.  The Berkeley software License Agreement
   4: # specifies the terms and conditions for redistribution.
   5: #
   6: #	@(#)sqrt.s	5.1 (Berkeley) 5/8/85
   7: #
   8: #
   9: # double sqrt(arg):revised July 18,1980
  10: # double arg
  11: # if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) }
  12: # W. Kahan's magic sqrt
  13: # coded by Heidi Stettner
  14: 
  15:         .set    EDOM,98
  16:         .text
  17:         .align  1
  18:         .globl  _sqrt
  19:         .globl  dsqrt_r5
  20:         .globl  _errno
  21: _sqrt:
  22:         .word   0x003c          # save r5,r4,r3,r2
  23:         bispsw  $0xe0
  24:         movd    4(ap),r0
  25:         bsbb    dsqrt_r5
  26:         ret
  27: dsqrt_r5:
  28:         movd    r0,r4
  29:         jleq    nonpos          #argument is not positive
  30:         movzwl  r4,r2
  31:         ashl    $-1,r2,r0
  32:         addw2   $0x203c,r0      #r0 has magic initial appx
  33: 
  34: # Do two steps of Heron's rule
  35: 
  36:         divf3   r0,r4,r2
  37:         addf2   r2,r0
  38:         subw2   $0x80,r0
  39: 
  40:         divf3   r0,r4,r2
  41:         addf2   r2,r0
  42:         subw2   $0x80,r0
  43: 
  44: 
  45: # Scale argument and approximation to prevent over/underflow
  46: # NOTE: The following four steps would not be necessary if underflow
  47: #       were gentle.
  48: 
  49:         bicw3   $0xffff807f,r4,r1
  50:         subw2   $0x4080,r1              # r1 contains scaling factor
  51:         subw2   r1,r4
  52:         movl    r0,r2
  53:         subw2   r1,r2
  54: 
  55: # Cubic step
  56: 
  57:         clrl    r1
  58:         clrl    r3
  59:         muld2   r0,r2
  60:         subd2   r2,r4
  61:         addw2   $0x100,r2
  62:         addd2   r4,r2
  63:         muld2   r0,r4
  64:         divd2   r2,r4
  65:         addw2   $0x80,r4
  66:         addd2   r4,r0
  67:         rsb
  68: nonpos:
  69:         jneq    negarg
  70:         clrd    r0              #argument is zero
  71:         rsb
  72: negarg:
  73:         movl    $EDOM,_errno
  74:         mnegd   r4,-(sp)
  75:         calls   $2,_sqrt
  76:         mnegd   r0,r0           # returns -sqrt(-arg)
  77:         ret

Defined functions