1 /* ef_sqrtf.c -- float version of e_sqrt.c.
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3 */
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 "fdlibm.h"
17
18 float
sqrtf(float x)19 sqrtf(float x)
20 {
21 float z;
22 __uint32_t r, hx;
23 __int32_t ix, s, q, m, t, i;
24
25 GET_FLOAT_WORD(ix, x);
26 hx = ix & 0x7fffffff;
27
28 /* take care of Inf and NaN */
29 if (!FLT_UWORD_IS_FINITE(hx)) {
30 if (ix < 0 && !isnanf(x))
31 return __math_invalidf(x); /* sqrt(-inf)=sNaN */
32 return x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf */
33 }
34
35 /* take care of zero and -ves */
36 if (FLT_UWORD_IS_ZERO(hx))
37 return x; /* sqrt(+-0) = +-0 */
38 if (ix < 0)
39 return __math_invalidf(x); /* sqrt(-ve) = sNaN */
40
41 /* normalize x */
42 m = (ix >> 23);
43 if (FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
44 for (i = 0; (ix & 0x00800000L) == 0; i++)
45 ix <<= 1;
46 m -= i - 1;
47 }
48 m -= 127; /* unbias exponent */
49 ix = (ix & 0x007fffffL) | 0x00800000L;
50 if (m & 1) /* odd m, double x to make it even */
51 ix += ix;
52 m >>= 1; /* m = [m/2] */
53
54 /* generate sqrt(x) bit by bit */
55 ix += ix;
56 q = s = 0; /* q = sqrt(x) */
57 r = 0x01000000L; /* r = moving bit from right to left */
58
59 while (r != 0) {
60 t = s + r;
61 if (t <= ix) {
62 s = t + r;
63 ix -= t;
64 q += r;
65 }
66 ix += ix;
67 r >>= 1;
68 }
69
70 if (ix != 0) {
71 FE_DECL_ROUND(rnd);
72 if (__is_nearest(rnd))
73 q += (q & 1);
74 else if (__is_upward(rnd))
75 q += 2;
76 }
77 ix = (q >> 1) + 0x3f000000L;
78 ix += (m << 23);
79 SET_FLOAT_WORD(z, ix);
80 return z;
81 }
82
83 _MATH_ALIAS_f_f(sqrt)
84
