Back to index

glibc  2.9
e_sqrtf.c
Go to the documentation of this file.
00001 /* e_sqrtf.c -- float version of e_sqrt.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice 
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_sqrtf.c,v 1.4 1995/05/10 20:46:19 jtc Exp $";
00018 #endif
00019 
00020 #include "math.h"
00021 #include "math_private.h"
00022 
00023 #ifdef __STDC__
00024 static const float   one    = 1.0, tiny=1.0e-30;
00025 #else
00026 static float  one    = 1.0, tiny=1.0e-30;
00027 #endif
00028 
00029 #ifdef __STDC__
00030        float __ieee754_sqrtf(float x)
00031 #else
00032        float __ieee754_sqrtf(x)
00033        float x;
00034 #endif
00035 {
00036        float z;
00037        int32_t sign = (int)0x80000000; 
00038        int32_t ix,s,q,m,t,i;
00039        u_int32_t r;
00040 
00041        GET_FLOAT_WORD(ix,x);
00042 
00043     /* take care of Inf and NaN */
00044        if((ix&0x7f800000)==0x7f800000) {                
00045            return x*x+x;           /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
00046                                       sqrt(-inf)=sNaN */
00047        } 
00048     /* take care of zero */
00049        if(ix<=0) {
00050            if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
00051            else if(ix<0)
00052               return (x-x)/(x-x);         /* sqrt(-ve) = sNaN */
00053        }
00054     /* normalize x */
00055        m = (ix>>23);
00056        if(m==0) {                         /* subnormal x */
00057            for(i=0;(ix&0x00800000)==0;i++) ix<<=1;
00058            m -= i-1;
00059        }
00060        m -= 127;     /* unbias exponent */
00061        ix = (ix&0x007fffff)|0x00800000;
00062        if(m&1)       /* odd m, double x to make it even */
00063            ix += ix;
00064        m >>= 1;      /* m = [m/2] */
00065 
00066     /* generate sqrt(x) bit by bit */
00067        ix += ix;
00068        q = s = 0;           /* q = sqrt(x) */
00069        r = 0x01000000;             /* r = moving bit from right to left */
00070 
00071        while(r!=0) {
00072            t = s+r; 
00073            if(t<=ix) { 
00074               s    = t+r; 
00075               ix  -= t; 
00076               q   += r; 
00077            } 
00078            ix += ix;
00079            r>>=1;
00080        }
00081 
00082     /* use floating add to find out rounding direction */
00083        if(ix!=0) {
00084            z = one-tiny; /* trigger inexact flag */
00085            if (z>=one) {
00086                z = one+tiny;
00087               if (z>one)
00088                   q += 2;
00089               else
00090                   q += (q&1);
00091            }
00092        }
00093        ix = (q>>1)+0x3f000000;
00094        ix += (m <<23);
00095        SET_FLOAT_WORD(z,ix);
00096        return z;
00097 }