1*34125Sbostic# Copyright (c) 1985 Regents of the University of California. 2*34125Sbostic# All rights reserved. 3*34125Sbostic# 4*34125Sbostic# Redistribution and use in source and binary forms are permitted 5*34125Sbostic# provided that this notice is preserved and that due credit is given 6*34125Sbostic# to the University of California at Berkeley. The name of the University 7*34125Sbostic# may not be used to endorse or promote products derived from this 8*34125Sbostic# software without specific prior written permission. This software 9*34125Sbostic# is provided ``as is'' without express or implied warranty. 10*34125Sbostic# 11*34125Sbostic# All recipients should regard themselves as participants in an ongoing 12*34125Sbostic# research project and hence should feel obligated to report their 13*34125Sbostic# experiences (good or bad) with these elementary function codes, using 14*34125Sbostic# the sendbug(8) program, to the authors. 15*34125Sbostic# 16*34125Sbostic# @(#)sqrt.s 5.2 (Berkeley) 04/29/88 17*34125Sbostic# 1824729Selefunt .data 1924729Selefunt .align 2 2024729Selefunt_sccsid: 21*34125Sbostic.asciz "@(#)sqrt.s 1.1 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) 04/29/88" 2224729Selefunt 2324729Selefunt/* 2424571Szliu * double sqrt(arg) revised August 15,1982 2524571Szliu * double arg; 2624571Szliu * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 2724571Szliu * if arg is a reserved operand it is returned as it is 2824571Szliu * W. Kahan's magic square root 2924571Szliu * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 3024571Szliu * 3124571Szliu * entry points:_d_sqrt address of double arg is on the stack 3224571Szliu * _sqrt double arg is on the stack 3324571Szliu */ 3424571Szliu .text 3524571Szliu .align 1 3624571Szliu .globl _sqrt 3724571Szliu .globl _d_sqrt 3824571Szliu .globl libm$dsqrt_r5 3924571Szliu .set EDOM,33 4024571Szliu 4124571Szliu_d_sqrt: 4224571Szliu .word 0x003c # save r5,r4,r3,r2 4324571Szliu movq *4(ap),r0 4424571Szliu jmp dsqrt2 4524571Szliu_sqrt: 4624571Szliu .word 0x003c # save r5,r4,r3,r2 4724571Szliu movq 4(ap),r0 4824571Szliudsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input 4924571Szliu jeql noexp # biased exponent is zero -> 0.0 or reserved 5024571Szliu bsbb libm$dsqrt_r5 5124571Szliunoexp: ret 5224571Szliu 5324571Szliu/* **************************** internal procedure */ 5424571Szliu 5524571Szliulibm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt 5624571Szliu # returns double square root scaled by 5724571Szliu # 2^r6 5824571Szliu 5924571Szliu movd r0,r4 6024571Szliu jleq nonpos # argument is not positive 6124571Szliu movzwl r4,r2 6224571Szliu ashl $-1,r2,r0 6324571Szliu addw2 $0x203c,r0 # r0 has magic initial approximation 6424571Szliu/* 6524571Szliu * Do two steps of Heron's rule 6624571Szliu * ((arg/guess) + guess) / 2 = better guess 6724571Szliu */ 6824571Szliu divf3 r0,r4,r2 6924571Szliu addf2 r2,r0 7024571Szliu subw2 $0x80,r0 # divide by two 7124571Szliu 7224571Szliu divf3 r0,r4,r2 7324571Szliu addf2 r2,r0 7424571Szliu subw2 $0x80,r0 # divide by two 7524571Szliu 7624571Szliu/* Scale argument and approximation to prevent over/underflow */ 7724571Szliu 7824571Szliu bicw3 $0x807f,r4,r1 7924571Szliu subw2 $0x4080,r1 # r1 contains scaling factor 8024571Szliu subw2 r1,r4 8124571Szliu movl r0,r2 8224571Szliu subw2 r1,r2 8324571Szliu 8424571Szliu/* Cubic step 8524571Szliu * 8624571Szliu * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, 8724571Szliu * a is approximation, and n is the original argument. 8824571Szliu * (let s be scale factor in the following comments) 8924571Szliu */ 9024571Szliu clrl r1 9124571Szliu clrl r3 9224571Szliu muld2 r0,r2 # r2:r3 = a*a/s 9324571Szliu subd2 r2,r4 # r4:r5 = n/s - a*a/s 9424571Szliu addw2 $0x100,r2 # r2:r3 = 4*a*a/s 9524571Szliu addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s 9624571Szliu muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s 9724571Szliu divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 9824571Szliu addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 9924571Szliu addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) 10024571Szliu rsb # DONE! 10124571Szliunonpos: 10224571Szliu jneq negarg 10324571Szliu ret # argument and root are zero 10424571Szliunegarg: 10524571Szliu pushl $EDOM 10624571Szliu calls $1,_infnan # generate the reserved op fault 10724571Szliu ret 108