| 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 |