1/* 2 * Copyright (c) 1987 Regents of the University of California. 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms are permitted 6 * provided that the above copyright notice and this paragraph are 7 * duplicated in all such forms and that any documentation, 8 * advertising materials, and other materials related to such 9 * distribution and use acknowledge that the software was developed 10 * by the University of California, Berkeley. The name of the 11 * University may not be used to endorse or promote products derived 12 * from this software without specific prior written permission. 13 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR 14 * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED 15 * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. 16 * 17 * All recipients should regard themselves as participants in an ongoing 18 * research project and hence should feel obligated to report their 19 * experiences (good or bad) with these elementary function codes, using 20 * the sendbug(8) program, to the authors. 21 */ 22 .data 23 .align 2 24_sccsid: 25.asciz "@(#)sqrt.s 5.4 (ucb.elefunt) 06/30/88" 26 27/* 28 * double sqrt(arg) revised August 15,1982 29 * double arg; 30 * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 31 * if arg is a reserved operand it is returned as it is 32 * W. Kahan's magic square root 33 * Coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82. 34 * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87. 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 2 41 .globl _sqrt 42 .globl _d_sqrt 43 .globl libm$dsqrt_r5 44 .set EDOM,33 45 46_d_sqrt: 47 .word 0x003c # save r2-r5 48 movl 4(fp),r2 49 movl (r2),r0 50 movl 4(r2),r1 # r0:r1 = x 51 brb 1f 52_sqrt: 53 .word 0x003c # save r2-r5 54 movl 4(fp),r0 55 movl 8(fp),r1 # r0:r1 = x 561: andl3 $0x7f800000,r0,r2 # r2 = biased exponent 57 bneq 2f 58 ret # biased exponent is zero -> 0 or reserved op. 59/* 60 * # internal procedure 61 * # ENTRY POINT FOR cdabs and cdsqrt 62 */ 63libm$dsqrt_r5: # returns double square root scaled by 2^r6 64 .word 0x0000 # save nothing 652: ldd r0 66 std r4 67 bleq nonpos # argument is not positive 68 andl3 $0xfffe0000,r4,r2 69 shar $1,r2,r0 70 addl2 $0x203c0000,r0 # r0 has magic initial approximation 71/* 72 * # Do two steps of Heron's rule 73 * # ((arg/guess)+guess)/2 = better guess 74 */ 75 ldf r4 76 divf r0 77 addf r0 78 stf r0 79 subl2 $0x800000,r0 # divide by two 80 ldf r4 81 divf r0 82 addf r0 83 stf r0 84 subl2 $0x800000,r0 # divide by two 85/* 86 * # Scale argument and approximation 87 * # to prevent over/underflow 88 */ 89 andl3 $0x7f800000,r4,r1 90 subl2 $0x40800000,r1 # r1 contains scaling factor 91 subl2 r1,r4 # r4:r5 = n/s 92 movl r0,r2 93 subl2 r1,r2 # r2 = a/s 94/* 95 * # Cubic step 96 * # b = a+2*a*(n-a*a)/(n+3*a*a) where 97 * # b is better approximation, a is approximation 98 * # and n is the original argument. 99 * # s := scale factor. 100 */ 101 clrl r1 # r0:r1 = a 102 clrl r3 # r2:r3 = a/s 103 ldd r0 # acc = a 104 muld r2 # acc = a*a/s 105 std r2 # r2:r3 = a*a/s 106 negd # acc = -a*a/s 107 addd r4 # acc = n/s-a*a/s 108 std r4 # r4:r5 = n/s-a*a/s 109 addl2 $0x1000000,r2 # r2:r3 = 4*a*a/s 110 ldd r2 # acc = 4*a*a/s 111 addd r4 # acc = n/s+3*a*a/s 112 std r2 # r2:r3 = n/s+3*a*a/s 113 ldd r0 # acc = a 114 muld r4 # acc = a*n/s-a*a*a/s 115 divd r2 # acc = a*(n-a*a)/(n+3*a*a) 116 std r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 117 addl2 $0x800000,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 118 ldd r4 # acc = 2*a*(n-a*a)/(n+3*a*a) 119 addd r0 # acc = a+2*a*(n-a*a)/(n+3*a*a) 120 std r0 # r0:r1 = a+2*a*(n-a*a)/(n+3*a*a) 121 ret # rsb 122nonpos: 123 bneq negarg 124 ret # argument and root are zero 125negarg: 126 pushl $EDOM 127 callf $8,_infnan # generate the reserved op fault 128 ret 129