1# Copyright (c) 1985 Regents of the University of California. 2# All rights reserved. 3# 4# Redistribution and use in source and binary forms are permitted 5# provided that the above copyright notice and this paragraph are 6# duplicated in all such forms and that any documentation, 7# advertising materials, and other materials related to such 8# distribution and use acknowledge that the software was developed 9# by the University of California, Berkeley. The name of the 10# University may not be used to endorse or promote products derived 11# from this software without specific prior written permission. 12# THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR 13# IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED 14# WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. 15# 16# All recipients should regard themselves as participants in an ongoing 17# research project and hence should feel obligated to report their 18# experiences (good or bad) with these elementary function codes, using 19# the sendbug(8) program, to the authors. 20# 21# @(#)sqrt.s 5.3 (Berkeley) 06/30/88 22# 23 .data 24 .align 2 25_sccsid: 26.asciz "@(#)sqrt.s 1.1 (Berkeley) 8/21/85; 5.3 (ucb.elefunt) 06/30/88" 27 28/* 29 * double sqrt(arg) revised August 15,1982 30 * double arg; 31 * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 32 * if arg is a reserved operand it is returned as it is 33 * W. Kahan's magic square root 34 * coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82 35 * 36 * entry points:_d_sqrt address of double arg is on the stack 37 * _sqrt double arg is on the stack 38 */ 39 .text 40 .align 1 41 .globl _sqrt 42 .globl _d_sqrt 43 .globl libm$dsqrt_r5 44 .set EDOM,33 45 46_d_sqrt: 47 .word 0x003c # save r5,r4,r3,r2 48 movq *4(ap),r0 49 jmp dsqrt2 50_sqrt: 51 .word 0x003c # save r5,r4,r3,r2 52 movq 4(ap),r0 53dsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input 54 jeql noexp # biased exponent is zero -> 0.0 or reserved 55 bsbb libm$dsqrt_r5 56noexp: ret 57 58/* **************************** internal procedure */ 59 60libm$dsqrt_r5: # ENTRY POINT FOR cdabs and cdsqrt 61 # returns double square root scaled by 62 # 2^r6 63 64 movd r0,r4 65 jleq nonpos # argument is not positive 66 movzwl r4,r2 67 ashl $-1,r2,r0 68 addw2 $0x203c,r0 # r0 has magic initial approximation 69/* 70 * Do two steps of Heron's rule 71 * ((arg/guess) + guess) / 2 = better guess 72 */ 73 divf3 r0,r4,r2 74 addf2 r2,r0 75 subw2 $0x80,r0 # divide by two 76 77 divf3 r0,r4,r2 78 addf2 r2,r0 79 subw2 $0x80,r0 # divide by two 80 81/* Scale argument and approximation to prevent over/underflow */ 82 83 bicw3 $0x807f,r4,r1 84 subw2 $0x4080,r1 # r1 contains scaling factor 85 subw2 r1,r4 86 movl r0,r2 87 subw2 r1,r2 88 89/* Cubic step 90 * 91 * b = a + 2*a*(n-a*a)/(n+3*a*a) where b is better approximation, 92 * a is approximation, and n is the original argument. 93 * (let s be scale factor in the following comments) 94 */ 95 clrl r1 96 clrl r3 97 muld2 r0,r2 # r2:r3 = a*a/s 98 subd2 r2,r4 # r4:r5 = n/s - a*a/s 99 addw2 $0x100,r2 # r2:r3 = 4*a*a/s 100 addd2 r4,r2 # r2:r3 = n/s + 3*a*a/s 101 muld2 r0,r4 # r4:r5 = a*n/s - a*a*a/s 102 divd2 r2,r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 103 addw2 $0x80,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 104 addd2 r4,r0 # r0:r1 = a + 2*a*(n-a*a)/(n+3*a*a) 105 rsb # DONE! 106nonpos: 107 jneq negarg 108 ret # argument and root are zero 109negarg: 110 pushl $EDOM 111 calls $1,_infnan # generate the reserved op fault 112 ret 113