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 sqrtf = __sqrtf
  30 
  31 #include "libm.h"
  32 
  33 #ifdef __INLINE
  34 
  35 extern float __inline_sqrtf(float);
  36 
  37 float
  38 sqrtf(float x) {
  39         return (__inline_sqrtf(x));
  40 }
  41 
  42 #else   /* defined(__INLINE) */
  43 
  44 static const float huge = 1.0e35F, tiny = 1.0e-35F, zero = 0.0f;
  45 
  46 float
  47 sqrtf(float x) {
  48         float   dz, w;
  49         int     *pw = (int *)&w;
  50         int     ix, j, r, q, m, n, s, t;
  51 
  52         w = x;
  53         ix = pw[0];
  54         if (ix <= 0) {
  55                 /* x is <= 0 or nan */
  56                 j = ix & 0x7fffffff;
  57                 if (j == 0)
  58                         return (w);
  59                 return ((w * zero) / zero);
  60         }
  61 
  62         if ((ix & 0x7f800000) == 0x7f800000) {
  63                 /* x is +inf or nan */
  64                 return (w * w);
  65         }
  66 
  67         m = ir_ilogb_(&w);
  68         n = -m;
  69         w = r_scalbn_(&w, (int *)&n);
  70         ix = (pw[0] & 0x007fffff) | 0x00800000;
  71         n = m / 2;
  72         if ((n + n) != m) {
  73                 ix = ix + ix;
  74                 m -= 1;
  75                 n = m / 2;
  76         }
  77 
  78         /* generate sqrt(x) bit by bit */
  79         ix <<= 1;
  80         q = s = 0;
  81         r = 0x01000000;
  82         for (j = 1; j <= 25; j++) {
  83                 t = s + r;
  84                 if (t <= ix) {
  85                         s = t + r;
  86                         ix -= t;
  87                         q += r;
  88                 }
  89                 ix <<= 1;
  90                 r >>= 1;
  91         }
  92         if (ix == 0)
  93                 goto done;
  94 
  95         /* raise inexact and determine the ambient rounding mode */
  96         dz = huge - tiny;
  97         if (dz < huge)
  98                 goto done;
  99         dz = huge + tiny;
 100         if (dz > huge)
 101                 q += 1;
 102         q += (q & 1);
 103 
 104 done:
 105         pw[0] = (q >> 1) + 0x3f000000;
 106         return (r_scalbn_(&w, (int *)&n));
 107 }
 108 
 109 #endif  /* defined(__INLINE) */