lib/msun/src/e_sqrtf.c

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