1*61318Sbostic# Copyright (c) 1985, 1993 2*61318Sbostic# The Regents of the University of California. All rights reserved. 334125Sbostic# 445308Sbostic# %sccs.include.redist.sh% 534125Sbostic# 6*61318Sbostic# @(#)sqrt.s 8.1 (Berkeley) 06/04/93 734125Sbostic# 824729Selefunt .data 924729Selefunt .align 2 1024729Selefunt_sccsid: 11*61318Sbostic.asciz "@(#)sqrt.s 1.1 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 06/04/93" 1224729Selefunt 1324729Selefunt/* 1424571Szliu * double sqrt(arg) revised August 15,1982 1524571Szliu * double arg; 1624571Szliu * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 1724571Szliu * if arg is a reserved operand it is returned as it is 1824571Szliu * W. Kahan's magic square root 1924571Szliu * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 2024571Szliu * 2124571Szliu * entry points:_d_sqrt address of double arg is on the stack 2224571Szliu * _sqrt double arg is on the stack 2324571Szliu */ 2424571Szliu .text 2524571Szliu .align 1 2624571Szliu .globl _sqrt 2724571Szliu .globl _d_sqrt 2824571Szliu .globl libm$dsqrt_r5 2924571Szliu .set EDOM,33 3024571Szliu 3124571Szliu_d_sqrt: 3224571Szliu .word 0x003c # save r5,r4,r3,r2 3324571Szliu movq *4(ap),r0 3424571Szliu jmp dsqrt2 3524571Szliu_sqrt: 3624571Szliu .word 0x003c # save r5,r4,r3,r2 3724571Szliu movq 4(ap),r0 3824571Szliudsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input 3924571Szliu jeql noexp # biased exponent is zero -> 0.0 or reserved 4024571Szliu bsbb libm$dsqrt_r5 4124571Szliunoexp: ret 4224571Szliu 4324571Szliu/* **************************** internal procedure */ 4424571Szliu 4524571Szliulibm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt 4624571Szliu # returns double square root scaled by 4724571Szliu # 2^r6 4824571Szliu 4924571Szliu movd r0,r4 5024571Szliu jleq nonpos # argument is not positive 5124571Szliu movzwl r4,r2 5224571Szliu ashl $-1,r2,r0 5324571Szliu addw2 $0x203c,r0 # r0 has magic initial approximation 5424571Szliu/* 5524571Szliu * Do two steps of Heron's rule 5624571Szliu * ((arg/guess) + guess) / 2 = better guess 5724571Szliu */ 5824571Szliu divf3 r0,r4,r2 5924571Szliu addf2 r2,r0 6024571Szliu subw2 $0x80,r0 # divide by two 6124571Szliu 6224571Szliu divf3 r0,r4,r2 6324571Szliu addf2 r2,r0 6424571Szliu subw2 $0x80,r0 # divide by two 6524571Szliu 6624571Szliu/* Scale argument and approximation to prevent over/underflow */ 6724571Szliu 6824571Szliu bicw3 $0x807f,r4,r1 6924571Szliu subw2 $0x4080,r1 # r1 contains scaling factor 7024571Szliu subw2 r1,r4 7124571Szliu movl r0,r2 7224571Szliu subw2 r1,r2 7324571Szliu 7424571Szliu/* Cubic step 7524571Szliu * 7624571Szliu * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, 7724571Szliu * a is approximation, and n is the original argument. 7824571Szliu * (let s be scale factor in the following comments) 7924571Szliu */ 8024571Szliu clrl r1 8124571Szliu clrl r3 8224571Szliu muld2 r0,r2 # r2:r3 = a*a/s 8324571Szliu subd2 r2,r4 # r4:r5 = n/s - a*a/s 8424571Szliu addw2 $0x100,r2 # r2:r3 = 4*a*a/s 8524571Szliu addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s 8624571Szliu muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s 8724571Szliu divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 8824571Szliu addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 8924571Szliu addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) 9024571Szliu rsb # DONE! 9124571Szliunonpos: 9224571Szliu jneq negarg 9324571Szliu ret # argument and root are zero 9424571Szliunegarg: 9524571Szliu pushl $EDOM 9624571Szliu calls $1,_infnan # generate the reserved op fault 9724571Szliu ret 98