124571Szliu/* 224571Szliu * Copyright (c) 1985 Regents of the University of California. 324571Szliu * 424571Szliu * Use and reproduction of this software are granted in accordance with 524571Szliu * the terms and conditions specified in the Berkeley Software License 624571Szliu * Agreement (in particular, this entails acknowledgement of the programs' 724571Szliu * source, and inclusion of this notice) with the additional understanding 824571Szliu * that all recipients should regard themselves as participants in an 924571Szliu * ongoing research project and hence should feel obligated to report 1024571Szliu * their experiences (good or bad) with these elementary function codes, 1124571Szliu * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12*24729Selefunt */ 13*24729Selefunt .data 14*24729Selefunt .align 2 15*24729Selefunt_sccsid: 16*24729Selefunt.asciz "@(#)sqrt.s 1.1 (Berkeley) 8/21/85; 1.3 (ucb.elefunt) 09/12/85" 17*24729Selefunt 18*24729Selefunt/* 1924571Szliu * double sqrt(arg) revised August 15,1982 2024571Szliu * double arg; 2124571Szliu * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 2224571Szliu * if arg is a reserved operand it is returned as it is 2324571Szliu * W. Kahan's magic square root 2424571Szliu * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 2524571Szliu * 2624571Szliu * entry points:_d_sqrt address of double arg is on the stack 2724571Szliu * _sqrt double arg is on the stack 2824571Szliu */ 2924571Szliu .text 3024571Szliu .align 1 3124571Szliu .globl _sqrt 3224571Szliu .globl _d_sqrt 3324571Szliu .globl libm$dsqrt_r5 3424571Szliu .set EDOM,33 3524571Szliu 3624571Szliu_d_sqrt: 3724571Szliu .word 0x003c # save r5,r4,r3,r2 3824571Szliu movq *4(ap),r0 3924571Szliu jmp dsqrt2 4024571Szliu_sqrt: 4124571Szliu .word 0x003c # save r5,r4,r3,r2 4224571Szliu movq 4(ap),r0 4324571Szliudsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input 4424571Szliu jeql noexp # biased exponent is zero -> 0.0 or reserved 4524571Szliu bsbb libm$dsqrt_r5 4624571Szliunoexp: ret 4724571Szliu 4824571Szliu/* **************************** internal procedure */ 4924571Szliu 5024571Szliulibm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt 5124571Szliu # returns double square root scaled by 5224571Szliu # 2^r6 5324571Szliu 5424571Szliu movd r0,r4 5524571Szliu jleq nonpos # argument is not positive 5624571Szliu movzwl r4,r2 5724571Szliu ashl $-1,r2,r0 5824571Szliu addw2 $0x203c,r0 # r0 has magic initial approximation 5924571Szliu/* 6024571Szliu * Do two steps of Heron's rule 6124571Szliu * ((arg/guess) + guess) / 2 = better guess 6224571Szliu */ 6324571Szliu divf3 r0,r4,r2 6424571Szliu addf2 r2,r0 6524571Szliu subw2 $0x80,r0 # divide by two 6624571Szliu 6724571Szliu divf3 r0,r4,r2 6824571Szliu addf2 r2,r0 6924571Szliu subw2 $0x80,r0 # divide by two 7024571Szliu 7124571Szliu/* Scale argument and approximation to prevent over/underflow */ 7224571Szliu 7324571Szliu bicw3 $0x807f,r4,r1 7424571Szliu subw2 $0x4080,r1 # r1 contains scaling factor 7524571Szliu subw2 r1,r4 7624571Szliu movl r0,r2 7724571Szliu subw2 r1,r2 7824571Szliu 7924571Szliu/* Cubic step 8024571Szliu * 8124571Szliu * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, 8224571Szliu * a is approximation, and n is the original argument. 8324571Szliu * (let s be scale factor in the following comments) 8424571Szliu */ 8524571Szliu clrl r1 8624571Szliu clrl r3 8724571Szliu muld2 r0,r2 # r2:r3 = a*a/s 8824571Szliu subd2 r2,r4 # r4:r5 = n/s - a*a/s 8924571Szliu addw2 $0x100,r2 # r2:r3 = 4*a*a/s 9024571Szliu addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s 9124571Szliu muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s 9224571Szliu divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 9324571Szliu addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 9424571Szliu addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) 9524571Szliu rsb # DONE! 9624571Szliunonpos: 9724571Szliu jneq negarg 9824571Szliu ret # argument and root are zero 9924571Szliunegarg: 10024571Szliu pushl $EDOM 10124571Szliu calls $1,_infnan # generate the reserved op fault 10224571Szliu ret 103