![[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 / national|
Return to sqrt.s CVS log | Up to [CSRG BSD Unix] / 43BSDReno / lib / libm / national |
1.1 root 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) 6/30/88
22: ;
23:
24: ; double sqrt(x)
25: ; double x;
26: ; IEEE double precision sqrt
27: ; code in NSC assembly by K.C. Ng
28: ; 12/13/85
29: ;
30: ; Method:
31: ; Use Kahan's trick to get 8 bits initial approximation
32: ; by integer shift and add/subtract. Then three Newton
33: ; iterations to get down error to within one ulp. Finally
34: ; twiddle the last bit to make to correctly rounded
35: ; according to the rounding mode.
36: ;
37: .vers 2
38: .text
39: .align 2
40: .globl _sqrt
41: _sqrt:
42: enter [r3,r4,r5,r6,r7],44
43: movl f4,tos
44: movl f6,tos
45: movd 2146435072,r2 ; r2 = 0x7ff00000
46: movl 8(fp),f0 ; f2 = x
47: movd 12(fp),r3 ; r3 = high part of x
48: movd r3,r4 ; make a copy of high part of x in r4
49: andd r2,r3 ; r3 become the bias exponent of x
50: cmpd r2,r3 ; if r3 = 0x7ff00000 then x is INF or NAN
51: bne L22
52: ; to see if x is INF
53: movd 8(fp),r0 ; r0 = low part of x
54: movd r4,r1 ; r1 is high part of x again
55: andd 0xfff00000,r1 ; mask off the sign and exponent of x
56: ord r0,r1 ; or with low part, if 0 then x is INF
57: cmpqd 0,r1 ;
58: bne L1 ; not 0; therefore x is NaN; return x.
59: cmpqd 0,r4 ; now x is Inf, is it +inf?
60: blt L1 ; +INF, return x
61: ; -INF, return NaN by doing 0/0
62: nan: movl 0f0.0,f0 ;
63: divl f0,f0
64: br L1
65: L22: ; now x is finite
66: cmpl 0f0.0,f0 ; x = 0 ?
67: beq L1 ; return x if x is +0 or -0
68: cmpqd 0,r4 ; Is x < 0 ?
69: bgt nan ; if x < 0 return NaN
70: movqd 0,r5 ; r5 == scalx initialize to zero
71: cmpqd 0,r3 ; is x is subnormal ? (r3 is the exponent)
72: bne L21 ; if x is normal, goto L21
73: movl L30,f2 ; f2 = 2**54
74: mull f2,f0 ; scale up x by 2**54
75: subd 0x1b00000,r5 ; off set the scale factor by -27 in exponent
76: L21:
77: ; now x is normal
78: ; notations:
79: ; r1 == copy of fsr
80: ; r2 == mask of e inexact enable flag
81: ; r3 == mask of i inexact flag
82: ; r4 == mask of r rounding mode
83: ; r5 == x's scale factor (already defined)
84:
85: movd 0x20,r2
86: movd 0x40,r3
87: movd 0x180,r4
88: sfsr r0 ; store fsr to r0
89: movd r0,r1 ; make a copy of fsr to r1
90: bicd [5,6,7,8],r0 ; clear e,i, and set r to round to nearest
91: lfsr r0
92:
93: ; begin to compute sqrt(x)
94: movl f0,-8(fp)
95: movd -4(fp),r0 ; r0 the high part of modified x
96: lshd -1,r0 ; r0 >> 1
97: addd 0x1ff80000,r0 ; add correction to r0 ...got 5 bits approx.
98: movd r0,r6
99: lshd -13,r6 ; r6 = r0>>-15
100: andd 0x7c,r6 ; obtain 4*leading 5 bits of r0
101: addrd L29,r7 ; r7 = address of L29 = table[0]
102: addd r6,r7 ; r6 = address of L29[r6] = table[r6]
103: subd 0(r7),r0 ; r0 = r0 - table[r6]
104: movd r0,-4(fp)
105: movl -8(fp),f2 ; now f2 = y approximate sqrt(f0) to 8 bits
106:
107: movl 0f0.5,f6 ; f6 = 0.5
108: movl f0,f4
109: divl f2,f4 ; t = x/y
110: addl f4,f2 ; y = y + x/y
111: mull f6,f2 ; y = 0.5(y+x/y) got 17 bits approx.
112: movl f0,f4
113: divl f2,f4 ; t = x/y
114: addl f4,f2 ; y = y + x/y
115: mull f6,f2 ; y = 0.5(y+x/y) got 35 bits approx.
116: movl f0,f4
117: divl f2,f4 ; t = x/y
118: subl f2,f4 ; t = x/y - y
119: mull f6,f4 ; t = 0.5(x/y-y)
120: addl f4,f2 ; y = y + 0.5(x/y -y)
121: ; now y approx. sqrt(x) to within 1 ulp
122:
123: ; twiddle last bit to force y correctly rounded
124: movd r1,r0 ; restore the old fsr
125: bicd [6,7,8],r0 ; clear inexact bit but retain inexact enable
126: ord 0x80,r0 ; set rounding mode to round to zero
127: lfsr r0
128:
129: movl f0,f4
130: divl f2,f4 ; f4 = x/y
131: sfsr r0
132: andd r3,r0 ; get the inexact flag
133: cmpqd 0,r0
134: bne L18
135: ; if x/y exact, then ...
136: cmpl f2,f4 ; if y == x/y
137: beq L2
138: movl f4,-8(fp)
139: subd 1,-8(fp)
140: subcd 0,-4(fp)
141: movl -8(fp),f4 ; f4 = f4 - ulp
142: L18:
143: bicd [6],r1
144: ord r3,r1 ; set inexact flag in r1
145:
146: andd r1,r4 ; r4 = the old rounding mode
147: cmpqd 0,r4 ; round to nearest?
148: bne L17
149: movl f4,-8(fp)
150: addd 1,-8(fp)
151: addcd 0,-4(fp)
152: movl -8(fp),f4 ; f4 = f4 + ulp
153: br L16
154: L17:
155: cmpd 0x100,r4 ; round to positive inf ?
156: bne L16
157: movl f4,-8(fp)
158: addd 1,-8(fp)
159: addcd 0,-4(fp)
160: movl -8(fp),f4 ; f4 = f4 + ulp
161:
162: movl f2,-8(fp)
163: addd 1,-8(fp)
164: addcd 0,-4(fp)
165: movl -8(fp),f2 ; f2 = f2 + ulp
166: L16:
167: addl f4,f2 ; y = y + t
168: subd 0x100000,r5 ; scalx = scalx - 1
169: L2:
170: movl f2,-8(fp)
171: addd r5,-4(fp)
172: movl -8(fp),f0
173: lfsr r1
174: L1:
175: movl tos,f6
176: movl tos,f4
177: exit [r3,r4,r5,r6,r7]
178: ret 0
179: .data
180: L28: .byte 64,40,35,41,115,113,114,116,46,99
181: .byte 9,49,46,49,32,40,117,99,98,46
182: .byte 101,108,101,102,117,110,116,41,32,57
183: .byte 47,49,57,47,56,53,0
184: L29: .blkb 4
185: .double 1204
186: .double 3062
187: .double 5746
188: .double 9193
189: .double 13348
190: .double 18162
191: .double 23592
192: .double 29598
193: .double 36145
194: .double 43202
195: .double 50740
196: .double 58733
197: .double 67158
198: .double 75992
199: .double 85215
200: .double 83599
201: .double 71378
202: .double 60428
203: .double 50647
204: .double 41945
205: .double 34246
206: .double 27478
207: .double 21581
208: .double 16499
209: .double 12183
210: .double 8588
211: .double 5674
212: .double 3403
213: .double 1742
214: .double 661
215: .double 130
216: L30: .blkb 4
217: .double 1129316352 ;L30: .double 0,0x43500000
218: L31: .blkb 4
219: .double 0x1ff00000
220: L32: .blkb 4
221: .double 0x5ff00000
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.