OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/k_tanf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/k_tanf.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 * Optimized by Bruce D. Evans. | 4 * Optimized by Bruce D. Evans. |
5 */ | 5 */ |
6 /* | 6 /* |
7 * ==================================================== | 7 * ==================================================== |
8 * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. | 8 * Copyright 2004 Sun Microsystems, Inc. All Rights Reserved. |
9 * | 9 * |
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 /* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */ | 18 /* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */ |
19 static const double T[] = { | 19 static const double T[] = { |
20 0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */ | 20 0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */ |
21 0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */ | 21 0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */ |
22 0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */ | 22 0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */ |
23 0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */ | 23 0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */ |
24 0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */ | 24 0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */ |
25 0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */ | 25 0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */ |
26 }; | 26 }; |
27 | 27 |
28 float __tandf(double x, int odd) | 28 float __tandf(double x, int odd) { |
29 { | 29 double_t z, r, w, s, t, u; |
30 » double_t z,r,w,s,t,u; | |
31 | 30 |
32 » z = x*x; | 31 z = x * x; |
33 » /* | 32 /* |
34 » * Split up the polynomial into small independent terms to give | 33 * Split up the polynomial into small independent terms to give |
35 » * opportunities for parallel evaluation. The chosen splitting is | 34 * opportunities for parallel evaluation. The chosen splitting is |
36 » * micro-optimized for Athlons (XP, X64). It costs 2 multiplications | 35 * micro-optimized for Athlons (XP, X64). It costs 2 multiplications |
37 » * relative to Horner's method on sequential machines. | 36 * relative to Horner's method on sequential machines. |
38 » * | 37 * |
39 » * We add the small terms from lowest degree up for efficiency on | 38 * We add the small terms from lowest degree up for efficiency on |
40 » * non-sequential machines (the lowest degree terms tend to be ready | 39 * non-sequential machines (the lowest degree terms tend to be ready |
41 » * earlier). Apart from this, we don't care about order of | 40 * earlier). Apart from this, we don't care about order of |
42 » * operations, and don't need to to care since we have precision to | 41 * operations, and don't need to to care since we have precision to |
43 » * spare. However, the chosen splitting is good for accuracy too, | 42 * spare. However, the chosen splitting is good for accuracy too, |
44 » * and would give results as accurate as Horner's method if the | 43 * and would give results as accurate as Horner's method if the |
45 » * small terms were added from highest degree down. | 44 * small terms were added from highest degree down. |
46 » */ | 45 */ |
47 » r = T[4] + z*T[5]; | 46 r = T[4] + z * T[5]; |
48 » t = T[2] + z*T[3]; | 47 t = T[2] + z * T[3]; |
49 » w = z*z; | 48 w = z * z; |
50 » s = z*x; | 49 s = z * x; |
51 » u = T[0] + z*T[1]; | 50 u = T[0] + z * T[1]; |
52 » r = (x + s*u) + (s*w)*(t + w*r); | 51 r = (x + s * u) + (s * w) * (t + w * r); |
53 » return odd ? -1.0/r : r; | 52 return odd ? -1.0 / r : r; |
54 } | 53 } |
OLD | NEW |