![[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] / 43BSDTahoe / usr.lib / libm / tahoe|
Return to sqrt.s CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / usr.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 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) 6/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
56: 1: 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: */
63: libm$dsqrt_r5: # returns double square root scaled by 2^r6
64: .word 0x0000 # save nothing
65: 2: 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
122: nonpos:
123: bneq negarg
124: ret # argument and root are zero
125: negarg:
126: pushl $EDOM
127: callf $8,_infnan # generate the reserved op fault
128: ret
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.