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 (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21 /*
22 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak sqrt = __sqrt
30
31 #include "libm.h"
32
33 #ifdef __INLINE
34
35 extern double __inline_sqrt(double);
36
37 double
38 sqrt(double x) {
39 double z = __inline_sqrt(x);
40
41 if (isnan(x))
42 return (z);
43 return ((x < 0.0)? _SVID_libm_err(x, x, 26) : z);
44 }
45
46 #else /* defined(__INLINE) */
47
48 /*
49 * Warning: This correctly rounded sqrt is extremely slow because it computes
50 * the sqrt bit by bit using integer arithmetic.
51 */
52
53 static const double big = 1.0e30, small = 1.0e-30;
54
55 double
56 sqrt(double x)
57 {
58 double z;
59 unsigned r, t1, s1, ix1, q1;
60 int ix0, s0, j, q, m, n, t;
61 int *px = (int *)&x, *pz = (int *)&z;
62
63 ix0 = px[HIWORD];
64 ix1 = px[LOWORD];
65 if ((ix0 & 0x7ff00000) == 0x7ff00000) { /* x is inf or NaN */
66 if (ix0 == 0xfff00000 && ix1 == 0)
67 return (_SVID_libm_err(x, x, 26));
68 return (x + x);
69 }
70 if (((ix0 & 0x7fffffff) | ix1) == 0) /* x is zero */
71 return (x);
72
73 /* extract exponent and significand */
74 m = ilogb(x);
75 z = scalbn(x, -m);
76 ix0 = (pz[HIWORD] & 0x000fffff) | 0x00100000;
77 ix1 = pz[LOWORD];
78 n = m >> 1;
79 if (n + n != m) {
80 ix0 = (ix0 << 1) | (ix1 >> 31);
81 ix1 <<= 1;
82 m -= 1;
83 }
84
85 /* generate sqrt(x) bit by bit */
86 ix0 = (ix0 << 1) | (ix1 >> 31);
87 ix1 <<= 1;
88 q = q1 = s0 = s1 = 0;
89 r = 0x00200000;
90
91 for (j = 1; j <= 22; j++) {
92 t = s0 + r;
93 if (t <= ix0) {
94 s0 = t + r;
95 ix0 -= t;
96 q += r;
97 }
98 ix0 = (ix0 << 1) | (ix1 >> 31);
99 ix1 <<= 1;
100 r >>= 1;
101 }
102
103 r = 0x80000000;
104 for (j = 1; j <= 32; j++) {
105 t1 = s1 + r;
106 t = s0;
107 if (t < ix0 || (t == ix0 && t1 <= ix1)) {
108 s1 = t1 + r;
109 if ((t1 & 0x80000000) == 0x80000000 &&
110 (s1 & 0x80000000) == 0)
111 s0 += 1;
112 ix0 -= t;
113 if (ix1 < t1)
114 ix0 -= 1;
115 ix1 -= t1;
116 q1 += r;
117 }
118 ix0 = (ix0 << 1) | (ix1 >> 31);
119 ix1 <<= 1;
120 r >>= 1;
121 }
122
123 /* round */
124 if ((ix0 | ix1) == 0)
125 goto done;
126 z = big - small; /* trigger inexact flag */
127 if (z < big)
128 goto done;
129 if (q1 == 0xffffffff) {
130 q1 = 0;
131 q += 1;
132 goto done;
133 }
134 z = big + small;
135 if (z > big) {
136 if (q1 == 0xfffffffe)
137 q += 1;
138 q1 += 2;
139 goto done;
140 }
141 q1 += (q1 & 1);
142 done:
143 pz[HIWORD] = (q >> 1) + 0x3fe00000;
144 pz[LOWORD] = q1 >> 1;
145 if ((q & 1) == 1)
146 pz[LOWORD] |= 0x80000000;
147 return (scalbn(z, n));
148 }
149
150 #endif /* defined(__INLINE) */