src/math/sqrtf.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
   2 /*
   3  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   4  */
   5 /*
   6  * ====================================================
   7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   8  *
   9  * Developed at SunPro, a Sun Microsystems, Inc. business.
  10  * Permission to use, copy, modify, and distribute this
  11  * software is freely granted, provided that this notice
  12  * is preserved.
  13  * ====================================================
  14  */
  15
  16 #include "libm.h"
  17
  18 static const float tiny = 1.0e-30;
  19
  20 float sqrtf(float x)
  21 {
  22         float z;
  23         int32_t sign = (int)0x80000000;
  24         int32_t ix,s,q,m,t,i;
  25         uint32_t r;
  26
  27         GET_FLOAT_WORD(ix, x);
  28
  29         /* take care of Inf and NaN */
  30         if ((ix&0x7f800000) == 0x7f800000)
  31                 return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
  32
  33         /* take care of zero */
  34         if (ix <= 0) {
  35                 if ((ix&~sign) == 0)
  36                         return x;  /* sqrt(+-0) = +-0 */
  37                 if (ix < 0)
  38                         return (x-x)/(x-x);  /* sqrt(-ve) = sNaN */
  39         }
  40         /* normalize x */
  41         m = ix>>23;
  42         if (m == 0) {  /* subnormal x */
  43                 for (i = 0; (ix&0x00800000) == 0; i++)
  44                         ix<<=1;
  45                 m -= i - 1;
  46         }
  47         m -= 127;  /* unbias exponent */
  48         ix = (ix&0x007fffff)|0x00800000;
  49         if (m&1)  /* odd m, double x to make it even */
  50                 ix += ix;
  51         m >>= 1;  /* m = [m/2] */
  52
  53         /* generate sqrt(x) bit by bit */
  54         ix += ix;
  55         q = s = 0;       /* q = sqrt(x) */
  56         r = 0x01000000;  /* r = moving bit from right to left */
  57
  58         while (r != 0) {
  59                 t = s + r;
  60                 if (t <= ix) {
  61                         s = t+r;
  62                         ix -= t;
  63                         q += r;
  64                 }
  65                 ix += ix;
  66                 r >>= 1;
  67         }
  68
  69         /* use floating add to find out rounding direction */
  70         if (ix != 0) {
  71                 z = 1.0f - tiny; /* raise inexact flag */
  72                 if (z >= 1.0f) {
  73                         z = 1.0f + tiny;
  74                         if (z > 1.0f)
  75                                 q += 2;
  76                         else
  77                                 q += q & 1;
  78                 }
  79         }
  80         ix = (q>>1) + 0x3f000000;
  81         ix += m << 23;
  82         SET_FLOAT_WORD(z, ix);
  83         return z;
  84 }