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