1 /* e_sqrtf.c -- float version of e_sqrt.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
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 $
9 * ====================================================
10 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
16 * ====================================================
20 #include "math_private.h"
23 static const float one = 1.0, tiny=1.0e-30;
25 static float one = 1.0, tiny=1.0e-30;
29 float __ieee754_sqrtf(float x)
31 float __ieee754_sqrtf(x)
36 int32_t sign = (int)0x80000000;
42 /* take care of Inf and NaN */
43 if((ix&0x7f800000)==0x7f800000) {
44 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
47 /* take care of zero */
49 if((ix&(~sign))==0) return x;/* sqrt(+-0) = +-0 */
51 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
55 if(m==0) { /* subnormal x */
56 for(i=0;(ix&0x00800000)==0;i++) ix<<=1;
59 m -= 127; /* unbias exponent */
60 ix = (ix&0x007fffff)|0x00800000;
61 if(m&1) /* odd m, double x to make it even */
63 m >>= 1; /* m = [m/2] */
65 /* generate sqrt(x) bit by bit */
67 q = s = 0; /* q = sqrt(x) */
68 r = 0x01000000; /* r = moving bit from right to left */
81 /* use floating add to find out rounding direction */
83 z = one-tiny; /* trigger inexact flag */
92 ix = (q>>1)+0x3f000000;