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) */ 151