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