![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to sqrt.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / tahoe|
Return to sqrt.s CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / tahoe |
1.1 root 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: (1) source distributions retain this entire copyright
7: * notice and comment, and (2) distributions including binaries display
8: * the following acknowledgement: ``This product includes software
9: * developed by the University of California, Berkeley and its contributors''
10: * in the documentation or other materials provided with the distribution
11: * and in all advertising materials mentioning features or use of this
12: * software. Neither the name of the University nor the names of its
13: * contributors may be used to endorse or promote products derived
14: * from this software without specific prior written permission.
15: * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
16: * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
17: * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
18: *
19: * All recipients should regard themselves as participants in an ongoing
20: * research project and hence should feel obligated to report their
21: * experiences (good or bad) with these elementary function codes, using
22: * the sendbug(8) program, to the authors.
23: */
24: .data
25: .align 2
26: _sccsid:
27: .asciz "@(#)sqrt.s 5.5 (ucb.elefunt) 6/1/90"
28:
29: /*
30: * double sqrt(arg) revised August 15,1982
31: * double arg;
32: * if(arg<0.0) { _errno = EDOM; return(<a reserved operand>); }
33: * if arg is a reserved operand it is returned as it is
34: * W. Kahan's magic square root
35: * Coded by Heidi Stettner and revised by Emile LeBlanc 8/18/82.
36: * Re-coded in tahoe assembly language by Z. Alex Liu 7/13/87.
37: *
38: * entry points:_d_sqrt address of double arg is on the stack
39: * _sqrt double arg is on the stack
40: */
41: .text
42: .align 2
43: .globl _sqrt
44: .globl _d_sqrt
45: .globl libm$dsqrt_r5
46: .set EDOM,33
47:
48: _d_sqrt:
49: .word 0x003c # save r2-r5
50: movl 4(fp),r2
51: movl (r2),r0
52: movl 4(r2),r1 # r0:r1 = x
53: brb 1f
54: _sqrt:
55: .word 0x003c # save r2-r5
56: movl 4(fp),r0
57: movl 8(fp),r1 # r0:r1 = x
58: 1: andl3 $0x7f800000,r0,r2 # r2 = biased exponent
59: bneq 2f
60: ret # biased exponent is zero -> 0 or reserved op.
61: /*
62: * # internal procedure
63: * # ENTRY POINT FOR cdabs and cdsqrt
64: */
65: libm$dsqrt_r5: # returns double square root scaled by 2^r6
66: .word 0x0000 # save nothing
67: 2: ldd r0
68: std r4
69: bleq nonpos # argument is not positive
70: andl3 $0xfffe0000,r4,r2
71: shar $1,r2,r0
72: addl2 $0x203c0000,r0 # r0 has magic initial approximation
73: /*
74: * # Do two steps of Heron's rule
75: * # ((arg/guess)+guess)/2 = better guess
76: */
77: ldf r4
78: divf r0
79: addf r0
80: stf r0
81: subl2 $0x800000,r0 # divide by two
82: ldf r4
83: divf r0
84: addf r0
85: stf r0
86: subl2 $0x800000,r0 # divide by two
87: /*
88: * # Scale argument and approximation
89: * # to prevent over/underflow
90: */
91: andl3 $0x7f800000,r4,r1
92: subl2 $0x40800000,r1 # r1 contains scaling factor
93: subl2 r1,r4 # r4:r5 = n/s
94: movl r0,r2
95: subl2 r1,r2 # r2 = a/s
96: /*
97: * # Cubic step
98: * # b = a+2*a*(n-a*a)/(n+3*a*a) where
99: * # b is better approximation, a is approximation
100: * # and n is the original argument.
101: * # s := scale factor.
102: */
103: clrl r1 # r0:r1 = a
104: clrl r3 # r2:r3 = a/s
105: ldd r0 # acc = a
106: muld r2 # acc = a*a/s
107: std r2 # r2:r3 = a*a/s
108: negd # acc = -a*a/s
109: addd r4 # acc = n/s-a*a/s
110: std r4 # r4:r5 = n/s-a*a/s
111: addl2 $0x1000000,r2 # r2:r3 = 4*a*a/s
112: ldd r2 # acc = 4*a*a/s
113: addd r4 # acc = n/s+3*a*a/s
114: std r2 # r2:r3 = n/s+3*a*a/s
115: ldd r0 # acc = a
116: muld r4 # acc = a*n/s-a*a*a/s
117: divd r2 # acc = a*(n-a*a)/(n+3*a*a)
118: std r4 # r4:r5 = a*(n-a*a)/(n+3*a*a)
119: addl2 $0x800000,r4 # r4:r5 = 2*a*(n-a*a)/(n+3*a*a)
120: ldd r4 # acc = 2*a*(n-a*a)/(n+3*a*a)
121: addd r0 # acc = a+2*a*(n-a*a)/(n+3*a*a)
122: std r0 # r0:r1 = a+2*a*(n-a*a)/(n+3*a*a)
123: ret # rsb
124: nonpos:
125: bneq negarg
126: ret # argument and root are zero
127: negarg:
128: pushl $EDOM
129: callf $8,_infnan # generate the reserved op fault
130: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.