1 /*
2 * CDDL HEADER START
3 *
4 * The contents of this file are subject to the terms of the
5 * Common Development and Distribution License, Version 1.0 only
6 * (the "License"). You may not use this file except in compliance
7 * with the License.
8 *
9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
10 * or http://www.opensolaris.org/os/licensing.
11 * See the License for the specific language governing permissions
12 * and limitations under the License.
13 *
14 * When distributing Covered Code, include this CDDL HEADER in each
15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
16 * If applicable, add the following below this CDDL HEADER, with the
17 * fields enclosed by brackets "[]" replaced with your own identifying
18 * information: Portions Copyright [yyyy] [name of copyright owner]
19 *
20 * CDDL HEADER END
21 */
22 /*
23 * Copyright 2003 Sun Microsystems, Inc. All rights reserved.
24 * Use is subject to license terms.
25 */
26
27 #include "quad.h"
28
29 static const double C[] = {
30 0.0,
31 0.5,
32 1.0,
33 68719476736.0,
34 536870912.0,
35 48.0,
36 16.0,
37 1.52587890625000000000e-05,
38 2.86102294921875000000e-06,
39 5.96046447753906250000e-08,
40 3.72529029846191406250e-09,
41 1.70530256582424044609e-13,
42 7.10542735760100185871e-15,
43 8.67361737988403547206e-19,
44 2.16840434497100886801e-19,
45 1.27054942088145050860e-21,
46 1.21169035041947413311e-27,
47 9.62964972193617926528e-35,
48 4.70197740328915003187e-38
49 };
50
51 #define zero C[0]
52 #define half C[1]
53 #define one C[2]
54 #define two36 C[3]
55 #define two29 C[4]
56 #define three2p4 C[5]
57 #define two4 C[6]
58 #define twom16 C[7]
59 #define three2m20 C[8]
60 #define twom24 C[9]
61 #define twom28 C[10]
62 #define three2m44 C[11]
63 #define twom47 C[12]
64 #define twom60 C[13]
65 #define twom62 C[14]
66 #define three2m71 C[15]
67 #define three2m91 C[16]
68 #define twom113 C[17]
69 #define twom124 C[18]
70
71 static const unsigned
72 fsr_re = 0x00000000u,
73 fsr_rn = 0xc0000000u;
74
75 #ifdef __sparcv9
76
77 /*
78 * _Qp_sqrt(pz, x) sets *pz = sqrt(*x).
79 */
80 void
_Qp_sqrt(union longdouble * pz,const union longdouble * x)81 _Qp_sqrt(union longdouble *pz, const union longdouble *x)
82
83 #else
84
85 /*
86 * _Q_sqrt(x) returns sqrt(*x).
87 */
88 union longdouble
89 _Q_sqrt(const union longdouble *x)
90
91 #endif /* __sparcv9 */
92
93 {
94 union longdouble z;
95 union xdouble u;
96 double c, d, rr, r[2], tt[3], xx[4], zz[5];
97 unsigned int xm, fsr, lx, wx[3];
98 unsigned int msw, frac2, frac3, frac4, rm;
99 int ex, ez;
100
101 if (QUAD_ISZERO(*x)) {
102 Z = *x;
103 QUAD_RETURN(Z);
104 }
105
106 xm = x->l.msw;
107
108 __quad_getfsrp(&fsr);
109
110 /* handle nan and inf cases */
111 if ((xm & 0x7fffffff) >= 0x7fff0000) {
112 if ((x->l.msw & 0xffff) | x->l.frac2 | x->l.frac3 |
113 x->l.frac4) {
114 if (!(x->l.msw & 0x8000)) {
115 /* snan, signal invalid */
116 if (fsr & FSR_NVM) {
117 __quad_fsqrtq(x, &Z);
118 } else {
119 Z = *x;
120 Z.l.msw |= 0x8000;
121 fsr = (fsr & ~FSR_CEXC) | FSR_NVA |
122 FSR_NVC;
123 __quad_setfsrp(&fsr);
124 }
125 QUAD_RETURN(Z);
126 }
127 Z = *x;
128 QUAD_RETURN(Z);
129 }
130 if (x->l.msw & 0x80000000) {
131 /* sqrt(-inf), signal invalid */
132 if (fsr & FSR_NVM) {
133 __quad_fsqrtq(x, &Z);
134 } else {
135 Z.l.msw = 0x7fffffff;
136 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff;
137 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC;
138 __quad_setfsrp(&fsr);
139 }
140 QUAD_RETURN(Z);
141 }
142 /* sqrt(inf), return inf */
143 Z = *x;
144 QUAD_RETURN(Z);
145 }
146
147 /* handle negative numbers */
148 if (xm & 0x80000000) {
149 if (fsr & FSR_NVM) {
150 __quad_fsqrtq(x, &Z);
151 } else {
152 Z.l.msw = 0x7fffffff;
153 Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff;
154 fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC;
155 __quad_setfsrp(&fsr);
156 }
157 QUAD_RETURN(Z);
158 }
159
160 /* now x is finite, positive */
161 __quad_setfsrp((unsigned *)&fsr_re);
162
163 /* get the normalized significand and exponent */
164 ex = (int)(xm >> 16);
165 lx = xm & 0xffff;
166 if (ex) {
167 lx |= 0x10000;
168 wx[0] = x->l.frac2;
169 wx[1] = x->l.frac3;
170 wx[2] = x->l.frac4;
171 } else {
172 if (lx | (x->l.frac2 & 0xfffe0000)) {
173 wx[0] = x->l.frac2;
174 wx[1] = x->l.frac3;
175 wx[2] = x->l.frac4;
176 ex = 1;
177 } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) {
178 lx = x->l.frac2;
179 wx[0] = x->l.frac3;
180 wx[1] = x->l.frac4;
181 wx[2] = 0;
182 ex = -31;
183 } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) {
184 lx = x->l.frac3;
185 wx[0] = x->l.frac4;
186 wx[1] = wx[2] = 0;
187 ex = -63;
188 } else {
189 lx = x->l.frac4;
190 wx[0] = wx[1] = wx[2] = 0;
191 ex = -95;
192 }
193 while ((lx & 0x10000) == 0) {
194 lx = (lx << 1) | (wx[0] >> 31);
195 wx[0] = (wx[0] << 1) | (wx[1] >> 31);
196 wx[1] = (wx[1] << 1) | (wx[2] >> 31);
197 wx[2] <<= 1;
198 ex--;
199 }
200 }
201 ez = ex - 0x3fff;
202 if (ez & 1) {
203 /* make exponent even */
204 lx = (lx << 1) | (wx[0] >> 31);
205 wx[0] = (wx[0] << 1) | (wx[1] >> 31);
206 wx[1] = (wx[1] << 1) | (wx[2] >> 31);
207 wx[2] <<= 1;
208 ez--;
209 }
210
211 /* extract the significands into doubles */
212 c = twom16;
213 xx[0] = (double)((int)lx) * c;
214
215 c *= twom24;
216 xx[0] += (double)((int)(wx[0] >> 8)) * c;
217
218 c *= twom24;
219 xx[1] = (double)((int)(((wx[0] << 16) | (wx[1] >> 16)) &
220 0xffffff)) * c;
221
222 c *= twom24;
223 xx[2] = (double)((int)(((wx[1] << 8) | (wx[2] >> 24)) &
224 0xffffff)) * c;
225
226 c *= twom24;
227 xx[3] = (double)((int)(wx[2] & 0xffffff)) * c;
228
229 /* approximate the divisor for the Newton iteration */
230 c = xx[0] + xx[1];
231 c = __quad_dp_sqrt(&c);
232 rr = half / c;
233
234 /* compute the first five "digits" of the square root */
235 zz[0] = (c + two29) - two29;
236 tt[0] = zz[0] + zz[0];
237 r[0] = (xx[0] - zz[0] * zz[0]) + xx[1];
238
239 zz[1] = (rr * (r[0] + xx[2]) + three2p4) - three2p4;
240 tt[1] = zz[1] + zz[1];
241 r[0] -= tt[0] * zz[1];
242 r[1] = xx[2] - zz[1] * zz[1];
243 c = (r[1] + three2m20) - three2m20;
244 r[0] += c;
245 r[1] = (r[1] - c) + xx[3];
246
247 zz[2] = (rr * (r[0] + r[1]) + three2m20) - three2m20;
248 tt[2] = zz[2] + zz[2];
249 r[0] -= tt[0] * zz[2];
250 r[1] -= tt[1] * zz[2];
251 c = (r[1] + three2m44) - three2m44;
252 r[0] += c;
253 r[1] = (r[1] - c) - zz[2] * zz[2];
254
255 zz[3] = (rr * (r[0] + r[1]) + three2m44) - three2m44;
256 r[0] = ((r[0] - tt[0] * zz[3]) + r[1]) - tt[1] * zz[3];
257 r[1] = -tt[2] * zz[3];
258 c = (r[1] + three2m91) - three2m91;
259 r[0] += c;
260 r[1] = (r[1] - c) - zz[3] * zz[3];
261
262 zz[4] = (rr * (r[0] + r[1]) + three2m71) - three2m71;
263
264 /* reduce to three doubles, making sure zz[1] is positive */
265 zz[0] += zz[1] - twom47;
266 zz[1] = twom47 + zz[2] + zz[3];
267 zz[2] = zz[4];
268
269 /* if the third term might lie on a rounding boundary, perturb it */
270 if (zz[2] == (twom62 + zz[2]) - twom62) {
271 /* here we just need to get the sign of the remainder */
272 c = (((((r[0] - tt[0] * zz[4]) - tt[1] * zz[4]) + r[1])
273 - tt[2] * zz[4]) - (zz[3] + zz[3]) * zz[4]) - zz[4] * zz[4];
274 if (c < zero)
275 zz[2] -= twom124;
276 else if (c > zero)
277 zz[2] += twom124;
278 }
279
280 /*
281 * propagate carries/borrows, using round-to-negative-infinity mode
282 * to make all terms nonnegative (note that we can't encounter a
283 * borrow so large that the roundoff is unrepresentable because
284 * we took care to make zz[1] positive above)
285 */
286 __quad_setfsrp(&fsr_rn);
287 c = zz[1] + zz[2];
288 zz[2] += (zz[1] - c);
289 zz[1] = c;
290 c = zz[0] + zz[1];
291 zz[1] += (zz[0] - c);
292 zz[0] = c;
293
294 /* adjust exponent and strip off integer bit */
295 ez = (ez >> 1) + 0x3fff;
296 zz[0] -= one;
297
298 /* the first 48 bits of fraction come from zz[0] */
299 u.d = d = two36 + zz[0];
300 msw = u.l.lo;
301 zz[0] -= (d - two36);
302
303 u.d = d = two4 + zz[0];
304 frac2 = u.l.lo;
305 zz[0] -= (d - two4);
306
307 /* the next 32 come from zz[0] and zz[1] */
308 u.d = d = twom28 + (zz[0] + zz[1]);
309 frac3 = u.l.lo;
310 zz[0] -= (d - twom28);
311
312 /* condense the remaining fraction; errors here won't matter */
313 c = zz[0] + zz[1];
314 zz[1] = ((zz[0] - c) + zz[1]) + zz[2];
315 zz[0] = c;
316
317 /* get the last word of fraction */
318 u.d = d = twom60 + (zz[0] + zz[1]);
319 frac4 = u.l.lo;
320 zz[0] -= (d - twom60);
321
322 /* keep track of what's left for rounding; note that the error */
323 /* in computing c will be non-negative due to rounding mode */
324 c = zz[0] + zz[1];
325
326 /* get the rounding mode */
327 rm = fsr >> 30;
328
329 /* round and raise exceptions */
330 fsr &= ~FSR_CEXC;
331 if (c != zero) {
332 fsr |= FSR_NXC;
333
334 /* decide whether to round the fraction up */
335 if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 ||
336 (c == twom113 && ((frac4 & 1) || (c - zz[0] != zz[1])))))) {
337 /* round up and renormalize if necessary */
338 if (++frac4 == 0)
339 if (++frac3 == 0)
340 if (++frac2 == 0)
341 if (++msw == 0x10000) {
342 msw = 0;
343 ez++;
344 }
345 }
346 }
347
348 /* stow the result */
349 z.l.msw = (ez << 16) | msw;
350 z.l.frac2 = frac2;
351 z.l.frac3 = frac3;
352 z.l.frac4 = frac4;
353
354 if ((fsr & FSR_CEXC) & (fsr >> 23)) {
355 __quad_setfsrp(&fsr);
356 __quad_fsqrtq(x, &Z);
357 } else {
358 Z = z;
359 fsr |= (fsr & 0x1f) << 5;
360 __quad_setfsrp(&fsr);
361 }
362 QUAD_RETURN(Z);
363 }
364