1# Copyright (c) 1985, 1993 2# The Regents of the University of California. All rights reserved. 3# 4# %sccs.include.redist.sh% 5# 6# @(#)sqrt.s 8.1 (Berkeley) 06/04/93 7# 8 .data 9 .align 2 10_sccsid: 11.asciz "@(#)sqrt.s 1.1 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 06/04/93" 12 13/* 14 * double sqrt(arg) revised August 15,1982 15 * double arg; 16 * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 17 * if arg is a reserved operand it is returned as it is 18 * W. Kahan's magic square root 19 * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 20 * 21 * entry points:_d_sqrt address of double arg is on the stack 22 * _sqrt double arg is on the stack 23 */ 24 .text 25 .align 1 26 .globl _sqrt 27 .globl _d_sqrt 28 .globl libm$dsqrt_r5 29 .set EDOM,33 30 31_d_sqrt: 32 .word 0x003c # save r5,r4,r3,r2 33 movq *4(ap),r0 34 jmp dsqrt2 35_sqrt: 36 .word 0x003c # save r5,r4,r3,r2 37 movq 4(ap),r0 38dsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input 39 jeql noexp # biased exponent is zero -> 0.0 or reserved 40 bsbb libm$dsqrt_r5 41noexp: ret 42 43/* **************************** internal procedure */ 44 45libm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt 46 # returns double square root scaled by 47 # 2^r6 48 49 movd r0,r4 50 jleq nonpos # argument is not positive 51 movzwl r4,r2 52 ashl $-1,r2,r0 53 addw2 $0x203c,r0 # r0 has magic initial approximation 54/* 55 * Do two steps of Heron's rule 56 * ((arg/guess) + guess) / 2 = better guess 57 */ 58 divf3 r0,r4,r2 59 addf2 r2,r0 60 subw2 $0x80,r0 # divide by two 61 62 divf3 r0,r4,r2 63 addf2 r2,r0 64 subw2 $0x80,r0 # divide by two 65 66/* Scale argument and approximation to prevent over/underflow */ 67 68 bicw3 $0x807f,r4,r1 69 subw2 $0x4080,r1 # r1 contains scaling factor 70 subw2 r1,r4 71 movl r0,r2 72 subw2 r1,r2 73 74/* Cubic step 75 * 76 * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, 77 * a is approximation, and n is the original argument. 78 * (let s be scale factor in the following comments) 79 */ 80 clrl r1 81 clrl r3 82 muld2 r0,r2 # r2:r3 = a*a/s 83 subd2 r2,r4 # r4:r5 = n/s - a*a/s 84 addw2 $0x100,r2 # r2:r3 = 4*a*a/s 85 addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s 86 muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s 87 divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 88 addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 89 addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) 90 rsb # DONE! 91nonpos: 92 jneq negarg 93 ret # argument and root are zero 94negarg: 95 pushl $EDOM 96 calls $1,_infnan # generate the reserved op fault 97 ret 98