* Copyright (c) 1985 Regents of the University of California.
* Use and reproduction of this software are granted in accordance with
* the terms and conditions specified in the Berkeley Software License
* Agreement (in particular, this entails acknowledgement of the programs'
* source, and inclusion of this notice) with the additional understanding
* that all recipients should regard themselves as participants in an
* ongoing research project and hence should feel obligated to report
* their experiences (good or bad) with these elementary function codes,
* using "sendbug 4bsd-bugs@BERKELEY", to the authors.
* @(#)sqrt.s 1.1 (Berkeley) 8/21/85; 1.2 (ucb.elefunt) %G%
* double sqrt(arg) revised August 15,1982
* if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
* if arg is a reserved operand it is returned as it is
* W. Kahan's magic square root
* coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82
* entry points:_d_sqrt address of double arg is on the stack
* _sqrt double arg is on the stack
.word 0x003c # save r5,r4,r3,r2
.word 0x003c # save r5,r4,r3,r2
dsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input
jeql noexp # biased exponent is zero -> 0.0 or reserved
/* **************************** internal procedure */
libm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt
# returns double square root scaled by
jleq nonpos # argument is not positive
addw2 $0x203c,r0 # r0 has magic initial approximation
* Do two steps of Heron's rule
* ((arg/guess) + guess) / 2 = better guess
subw2 $0x80,r0 # divide by two
subw2 $0x80,r0 # divide by two
/* Scale argument and approximation to prevent over/underflow */
subw2 $0x4080,r1 # r1 contains scaling factor
* b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation,
* a is approximation, and n is the original argument.
* (let s be scale factor in the following comments)
muld2 r0,r2 # r2:r3 = a*a/s
subd2 r2,r4 # r4:r5 = n/s - a*a/s
addw2 $0x100,r2 # r2:r3 = 4*a*a/s
addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s
muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s
divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a)
addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a)
ret # argument and root are zero
calls $1,_infnan # generate the reserved op fault