OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.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 * Debugged and optimized by Bruce D. Evans. | 4 * Debugged and optimized by Bruce D. Evans. |
5 */ | 5 */ |
6 /* | 6 /* |
7 * ==================================================== | 7 * ==================================================== |
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
9 * | 9 * |
10 * Developed at SunPro, a Sun Microsystems, Inc. business. | 10 * Developed at SunPro, a Sun Microsystems, Inc. business. |
11 * Permission to use, copy, modify, and distribute this | 11 * Permission to use, copy, modify, and distribute this |
12 * software is freely granted, provided that this notice | 12 * software is freely granted, provided that this notice |
13 * is preserved. | 13 * is preserved. |
14 * ==================================================== | 14 * ==================================================== |
15 */ | 15 */ |
16 /* __rem_pio2f(x,y) | 16 /* __rem_pio2f(x,y) |
17 * | 17 * |
18 * return the remainder of x rem pi/2 in *y | 18 * return the remainder of x rem pi/2 in *y |
19 * use double precision for everything except passing x | 19 * use double precision for everything except passing x |
20 * use __rem_pio2_large() for large x | 20 * use __rem_pio2_large() for large x |
21 */ | 21 */ |
22 | 22 |
23 #include "libm.h" | 23 #include "libm.h" |
24 | 24 |
25 #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 | 25 #if FLT_EVAL_METHOD == 0 || FLT_EVAL_METHOD == 1 |
26 #define EPS DBL_EPSILON | 26 #define EPS DBL_EPSILON |
27 #elif FLT_EVAL_METHOD==2 | 27 #elif FLT_EVAL_METHOD == 2 |
28 #define EPS LDBL_EPSILON | 28 #define EPS LDBL_EPSILON |
29 #endif | 29 #endif |
30 | 30 |
31 /* | 31 /* |
32 * invpio2: 53 bits of 2/pi | 32 * invpio2: 53 bits of 2/pi |
33 * pio2_1: first 25 bits of pi/2 | 33 * pio2_1: first 25 bits of pi/2 |
34 * pio2_1t: pi/2 - pio2_1 | 34 * pio2_1t: pi/2 - pio2_1 |
35 */ | 35 */ |
36 static const double | 36 static const double toint = 1.5 / EPS, |
37 toint = 1.5/EPS, | 37 invpio2 = |
38 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ | 38 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ |
39 pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ | 39 pio2_1 = 1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */ |
40 pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ | 40 pio2_1t = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */ |
41 | 41 |
42 int __rem_pio2f(float x, double *y) | 42 int __rem_pio2f(float x, double* y) { |
43 { | 43 union { |
44 » union {float f; uint32_t i;} u = {x}; | 44 float f; |
45 » double tx[1],ty[1]; | 45 uint32_t i; |
46 » double_t fn; | 46 } u = {x}; |
47 » uint32_t ix; | 47 double tx[1], ty[1]; |
48 » int n, sign, e0; | 48 double_t fn; |
| 49 uint32_t ix; |
| 50 int n, sign, e0; |
49 | 51 |
50 » ix = u.i & 0x7fffffff; | 52 ix = u.i & 0x7fffffff; |
51 » /* 25+53 bit pi is good enough for medium size */ | 53 /* 25+53 bit pi is good enough for medium size */ |
52 » if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ | 54 if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ |
53 » » /* Use a specialized rint() to get fn. Assume round-to-nearest.
*/ | 55 /* Use a specialized rint() to get fn. Assume round-to-nearest. */ |
54 » » fn = (double_t)x*invpio2 + toint - toint; | 56 fn = (double_t)x * invpio2 + toint - toint; |
55 » » n = (int32_t)fn; | 57 n = (int32_t)fn; |
56 » » *y = x - fn*pio2_1 - fn*pio2_1t; | 58 *y = x - fn * pio2_1 - fn * pio2_1t; |
57 » » return n; | 59 return n; |
58 » } | 60 } |
59 » if(ix>=0x7f800000) { /* x is inf or NaN */ | 61 if (ix >= 0x7f800000) { /* x is inf or NaN */ |
60 » » *y = x-x; | 62 *y = x - x; |
61 » » return 0; | 63 return 0; |
62 » } | 64 } |
63 » /* scale x into [2^23, 2^24-1] */ | 65 /* scale x into [2^23, 2^24-1] */ |
64 » sign = u.i>>31; | 66 sign = u.i >> 31; |
65 » e0 = (ix>>23) - (0x7f+23); /* e0 = ilogb(|x|)-23, positive */ | 67 e0 = (ix >> 23) - (0x7f + 23); /* e0 = ilogb(|x|)-23, positive */ |
66 » u.i = ix - (e0<<23); | 68 u.i = ix - (e0 << 23); |
67 » tx[0] = u.f; | 69 tx[0] = u.f; |
68 » n = __rem_pio2_large(tx,ty,e0,1,0); | 70 n = __rem_pio2_large(tx, ty, e0, 1, 0); |
69 » if (sign) { | 71 if (sign) { |
70 » » *y = -ty[0]; | 72 *y = -ty[0]; |
71 » » return -n; | 73 return -n; |
72 » } | 74 } |
73 » *y = ty[0]; | 75 *y = ty[0]; |
74 » return n; | 76 return n; |
75 } | 77 } |
OLD | NEW |