4.3BSD - /usr/src/usr.lib/libm/VAX/sqrt.s

   1: /*
   2:  * Copyright (c) 1985 Regents of the University of California.
   3:  *
   4:  * Use and reproduction of this software are granted  in  accordance  with
   5:  * the terms and conditions specified in  the  Berkeley  Software  License
   6:  * Agreement (in particular, this entails acknowledgement of the programs'
   7:  * source, and inclusion of this notice) with the additional understanding
   8:  * that  all  recipients  should regard themselves as participants  in  an
   9:  * ongoing  research  project and hence should  feel  obligated  to report
  10:  * their  experiences (good or bad) with these elementary function  codes,
  11:  * using "sendbug 4bsd-bugs@BERKELEY", to the authors.
  12:  *
  13:  *
  14:  * @(#)sqrt.s	1.1 (Berkeley) 8/21/85
  15:  *
  16:  * double sqrt(arg)   revised August 15,1982
  17:  * double arg;
  18:  * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
  19:  * if arg is a reserved operand it is returned as it is
  20:  * W. Kahan's magic square root
  21:  * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
  22:  *
  23:  * entry points:_d_sqrt		address of double arg is on the stack
  24:  *		_sqrt		double arg is on the stack
  25:  */
  26:         .text
  27:         .align  1
  28:         .globl  _sqrt
  29:         .globl  _d_sqrt
  30:         .globl  libm$dsqrt_r5
  31:         .set    EDOM,33
  32: 
  33: _d_sqrt:
  34:         .word   0x003c          # save r5,r4,r3,r2
  35:         movq    *4(ap),r0
  36:         jmp     dsqrt2
  37: _sqrt:
  38:         .word   0x003c          # save r5,r4,r3,r2
  39:         movq    4(ap),r0
  40: dsqrt2: bicw3   $0x807f,r0,r2   # check exponent of input
  41:         jeql    noexp           # biased exponent is zero -> 0.0 or reserved
  42:         bsbb    libm$dsqrt_r5
  43: noexp:  ret
  44: 
  45: /* **************************** internal procedure */
  46: 
  47: libm$dsqrt_r5:                  # ENTRY POINT FOR cdabs and cdsqrt
  48:                                 # returns double square root scaled by
  49:                                 # 2^r6
  50: 
  51:         movd    r0,r4
  52:         jleq    nonpos          # argument is not positive
  53:         movzwl  r4,r2
  54:         ashl    $-1,r2,r0
  55:         addw2   $0x203c,r0      # r0 has magic initial approximation
  56: /*
  57:  * Do two steps of Heron's rule
  58:  * ((arg/guess) + guess) / 2 = better guess
  59:  */
  60:         divf3   r0,r4,r2
  61:         addf2   r2,r0
  62:         subw2   $0x80,r0        # divide by two
  63: 
  64:         divf3   r0,r4,r2
  65:         addf2   r2,r0
  66:         subw2   $0x80,r0        # divide by two
  67: 
  68: /* Scale argument and approximation to prevent over/underflow */
  69: 
  70:         bicw3   $0x807f,r4,r1
  71:         subw2   $0x4080,r1              # r1 contains scaling factor
  72:         subw2   r1,r4
  73:         movl    r0,r2
  74:         subw2   r1,r2
  75: 
  76: /* Cubic step
  77:  *
  78:  * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
  79:  * a is approximation, and n is the original argument.
  80:  * (let s be scale factor in the following comments)
  81:  */
  82:         clrl    r1
  83:         clrl    r3
  84:         muld2   r0,r2                   # r2:r3 = a*a/s
  85:         subd2   r2,r4                   # r4:r5 = n/s - a*a/s
  86:         addw2   $0x100,r2               # r2:r3 = 4*a*a/s
  87:         addd2   r4,r2                   # r2:r3 = n/s + 3*a*a/s
  88:         muld2   r0,r4                   # r4:r5 = a*n/s - a*a*a/s
  89:         divd2   r2,r4                   # r4:r5 = a*(n-a*a)/(n+3*a*a)
  90:         addw2   $0x80,r4                # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
  91:         addd2   r4,r0                   # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a)
  92:         rsb                             # DONE!
  93: nonpos:
  94:         jneq    negarg
  95:         ret                     # argument and root are zero
  96: negarg:
  97:         pushl   $EDOM
  98:         calls   $1,_infnan      # generate the reserved op fault
  99:         ret

Defined functions