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