1/* 2 * Copyright (c) 1987 Regents of the University of California. 3 * 4 * Use and reproduction of this software are granted in accordance with 5 * the terms and conditions specified in the Berkeley Software License 6 * Agreement (in particular, this entails acknowledgement of the programs' 7 * source, and inclusion of this notice) with the additional understanding 8 * that all recipients should regard themselves as participants in an 9 * ongoing research project and hence should feel obligated to report 10 * their experiences (good or bad) with these elementary function codes, 11 * using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12 */ 13 .data 14 .align 2 15_sccsid: 16.asciz "@(#)sqrt.s 1.3 (ucb.elefunt) 07/13/87" 17 18/* 19 * double sqrt(arg) revised August 15,1982 20 * double arg; 21 * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); } 22 * if arg is a reserved operand it is returned as it is 23 * W. Kahan's magic square root 24 * Coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82. 25 * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87. 26 * 27 * entry points:_d_sqrt address of double arg is on the stack 28 * _sqrt double arg is on the stack 29 */ 30 .text 31 .align 2 32 .globl _sqrt 33 .globl _d_sqrt 34 .globl libm$dsqrt_r5 35 .set EDOM,33 36 37_d_sqrt: 38 .word 0x003c # save r2-r5 39 movl 4(fp),r2 40 movl (r2),r0 41 movl 4(r2),r1 # r0:r1 = x 42 brb 1f 43_sqrt: 44 .word 0x003c # save r2-r5 45 movl 4(fp),r0 46 movl 8(fp),r1 # r0:r1 = x 471: andl3 $0x7f800000,r0,r2 # r2 = biased exponent 48 bneq 2f 49 ret # biased exponent is zero -> 0 or reserved op. 50/* 51 * # internal procedure 52 * # ENTRY POINT FOR cdabs and cdsqrt 53 */ 54libm$dsqrt_r5: # returns double square root scaled by 2^r6 55 .word 0x0000 # save nothing 562: ldd r0 57 std r4 58 bleq nonpos # argument is not positive 59 andl3 $0xfffe0000,r4,r2 60 shar $1,r2,r0 61 addl2 $0x203c0000,r0 # r0 has magic initial approximation 62/* 63 * # Do two steps of Heron's rule 64 * # ((arg/guess)+guess)/2 = better guess 65 */ 66 ldf r4 67 divf r0 68 addf r0 69 stf r0 70 subl2 $0x800000,r0 # divide by two 71 ldf r4 72 divf r0 73 addf r0 74 stf r0 75 subl2 $0x800000,r0 # divide by two 76/* 77 * # Scale argument and approximation 78 * # to prevent over/underflow 79 */ 80 andl3 $0x7f800000,r4,r1 81 subl2 $0x40800000,r1 # r1 contains scaling factor 82 subl2 r1,r4 # r4:r5 = n/s 83 movl r0,r2 84 subl2 r1,r2 # r2 = a/s 85/* 86 * # Cubic step 87 * # b = a+2*a*(n-a*a)/(n+3*a*a) where 88 * # b is better approximation, a is approximation 89 * # and n is the original argument. 90 * # s := scale factor. 91 */ 92 clrl r1 # r0:r1 = a 93 clrl r3 # r2:r3 = a/s 94 ldd r0 # acc = a 95 muld r2 # acc = a*a/s 96 std r2 # r2:r3 = a*a/s 97 negd # acc = -a*a/s 98 addd r4 # acc = n/s-a*a/s 99 std r4 # r4:r5 = n/s-a*a/s 100 addl2 $0x1000000,r2 # r2:r3 = 4*a*a/s 101 ldd r2 # acc = 4*a*a/s 102 addd r4 # acc = n/s+3*a*a/s 103 std r2 # r2:r3 = n/s+3*a*a/s 104 ldd r0 # acc = a 105 muld r4 # acc = a*n/s-a*a*a/s 106 divd r2 # acc = a*(n-a*a)/(n+3*a*a) 107 std r4 # r4:r5 = a*(n-a*a)/(n+3*a*a) 108 addl2 $0x800000,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a) 109 ldd r4 # acc = 2*a*(n-a*a)/(n+3*a*a) 110 addd r0 # acc = a+2*a*(n-a*a)/(n+3*a*a) 111 std r0 # r0:r1 = a+2*a*(n-a*a)/(n+3*a*a) 112 ret # rsb 113nonpos: 114 bneq negarg 115 ret # argument and root are zero 116negarg: 117 pushl $EDOM 118 callf $8,_infnan # generate the reserved op fault 119 ret 120