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