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