OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/ld80/k_cosl.c */ |
2 /* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */ | 2 /* origin: FreeBSD /usr/src/lib/msun/ld128/k_cosl.c */ |
3 /* | 3 /* |
4 * ==================================================== | 4 * ==================================================== |
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
6 * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. | 6 * Copyright (c) 2008 Steven G. Kargl, David Schultz, Bruce D. Evans. |
7 * | 7 * |
8 * Developed at SunSoft, a Sun Microsystems, Inc. business. | 8 * Developed at SunSoft, a Sun Microsystems, Inc. business. |
9 * Permission to use, copy, modify, and distribute this | 9 * Permission to use, copy, modify, and distribute this |
10 * software is freely granted, provided that this notice | 10 * software is freely granted, provided that this notice |
11 * is preserved. | 11 * is preserved. |
12 * ==================================================== | 12 * ==================================================== |
13 */ | 13 */ |
14 | 14 |
15 | |
16 #include "libm.h" | 15 #include "libm.h" |
17 | 16 |
18 #if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 17 #if (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
19 #if LDBL_MANT_DIG == 64 | 18 #if LDBL_MANT_DIG == 64 |
20 /* | 19 /* |
21 * ld80 version of __cos.c. See __cos.c for most comments. | 20 * ld80 version of __cos.c. See __cos.c for most comments. |
22 */ | 21 */ |
23 /* | 22 /* |
24 * Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]: | 23 * Domain [-0.7854, 0.7854], range ~[-2.43e-23, 2.425e-23]: |
25 * |cos(x) - c(x)| < 2**-75.1 | 24 * |cos(x) - c(x)| < 2**-75.1 |
(...skipping 10 matching lines...) Expand all Loading... |
36 * a mediocre approximation for the coefficient of x^4, but a rounding | 35 * a mediocre approximation for the coefficient of x^4, but a rounding |
37 * error of 0.5 ulps for this coefficient would only contribute ~0.01 | 36 * error of 0.5 ulps for this coefficient would only contribute ~0.01 |
38 * ulps to the final error, so this is unimportant. Rounding errors in | 37 * ulps to the final error, so this is unimportant. Rounding errors in |
39 * higher coefficients are even less important. | 38 * higher coefficients are even less important. |
40 * | 39 * |
41 * In fact, coefficients above the x^4 one only need to have 53-bit | 40 * In fact, coefficients above the x^4 one only need to have 53-bit |
42 * precision, and this is more efficient. We get this optimization | 41 * precision, and this is more efficient. We get this optimization |
43 * almost for free from the complications needed to search for the best | 42 * almost for free from the complications needed to search for the best |
44 * higher coefficients. | 43 * higher coefficients. |
45 */ | 44 */ |
46 static const long double | 45 static const long double C1 = |
47 C1 = 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ | 46 0.0416666666666666666136L; /* 0xaaaaaaaaaaaaaa9b.0p-68 */ |
48 static const double | 47 static const double C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ |
49 C2 = -0.0013888888888888874, /* -0x16c16c16c16c10.0p-62 */ | 48 C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ |
50 C3 = 0.000024801587301571716, /* 0x1a01a01a018e22.0p-68 */ | 49 C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ |
51 C4 = -0.00000027557319215507120, /* -0x127e4fb7602f22.0p-74 */ | 50 C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ |
52 C5 = 0.0000000020876754400407278, /* 0x11eed8caaeccf1.0p-81 */ | 51 C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ |
53 C6 = -1.1470297442401303e-11, /* -0x19393412bd1529.0p-89 */ | 52 C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */ |
54 C7 = 4.7383039476436467e-14; /* 0x1aac9d9af5c43e.0p-97 */ | 53 #define POLY(z) \ |
55 #define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*C7))))))) | 54 (z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * (C6 + z * C7))))))) |
56 #elif LDBL_MANT_DIG == 113 | 55 #elif LDBL_MANT_DIG == 113 |
57 /* | 56 /* |
58 * ld128 version of __cos.c. See __cos.c for most comments. | 57 * ld128 version of __cos.c. See __cos.c for most comments. |
59 */ | 58 */ |
60 /* | 59 /* |
61 * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]: | 60 * Domain [-0.7854, 0.7854], range ~[-1.80e-37, 1.79e-37]: |
62 * |cos(x) - c(x))| < 2**-122.0 | 61 * |cos(x) - c(x))| < 2**-122.0 |
63 * | 62 * |
64 * 113-bit precision requires more care than 64-bit precision, since | 63 * 113-bit precision requires more care than 64-bit precision, since |
65 * simple methods give a minimax polynomial with coefficient for x^2 | 64 * simple methods give a minimax polynomial with coefficient for x^2 |
66 * that is 1 ulp below 0.5, but we want it to be precisely 0.5. See | 65 * that is 1 ulp below 0.5, but we want it to be precisely 0.5. See |
67 * above for more details. | 66 * above for more details. |
68 */ | 67 */ |
69 static const long double | 68 static const long double C1 = 0.04166666666666666666666666666666658424671L, |
70 C1 = 0.04166666666666666666666666666666658424671L, | 69 C2 = -0.001388888888888888888888888888863490893732L, |
71 C2 = -0.001388888888888888888888888888863490893732L, | 70 C3 = 0.00002480158730158730158730158600795304914210L, |
72 C3 = 0.00002480158730158730158730158600795304914210L, | 71 C4 = -0.2755731922398589065255474947078934284324e-6L, |
73 C4 = -0.2755731922398589065255474947078934284324e-6L, | 72 C5 = 0.2087675698786809897659225313136400793948e-8L, |
74 C5 = 0.2087675698786809897659225313136400793948e-8L, | 73 C6 = -0.1147074559772972315817149986812031204775e-10L, |
75 C6 = -0.1147074559772972315817149986812031204775e-10L, | 74 C7 = 0.4779477332386808976875457937252120293400e-13L; |
76 C7 = 0.4779477332386808976875457937252120293400e-13L; | 75 static const double C8 = -0.1561920696721507929516718307820958119868e-15, |
77 static const double | 76 C9 = 0.4110317413744594971475941557607804508039e-18, |
78 C8 = -0.1561920696721507929516718307820958119868e-15, | 77 C10 = -0.8896592467191938803288521958313920156409e-21, |
79 C9 = 0.4110317413744594971475941557607804508039e-18, | 78 C11 = 0.1601061435794535138244346256065192782581e-23; |
80 C10 = -0.8896592467191938803288521958313920156409e-21, | 79 #define POLY(z) \ |
81 C11 = 0.1601061435794535138244346256065192782581e-23; | 80 (z * \ |
82 #define POLY(z) (z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*(C6+z*(C7+ \ | 81 (C1 + \ |
83 » z*(C8+z*(C9+z*(C10+z*C11))))))))))) | 82 z * (C2 + \ |
| 83 z * (C3 + \ |
| 84 z * (C4 + \ |
| 85 z * (C5 + \ |
| 86 z * (C6 + \ |
| 87 z * (C7 + \ |
| 88 z * (C8 + \ |
| 89 z * (C9 + z * (C10 + z * C11))))))))))) |
84 #endif | 90 #endif |
85 | 91 |
86 long double __cosl(long double x, long double y) | 92 long double __cosl(long double x, long double y) { |
87 { | 93 long double hz, z, r, w; |
88 » long double hz,z,r,w; | |
89 | 94 |
90 » z = x*x; | 95 z = x * x; |
91 » r = POLY(z); | 96 r = POLY(z); |
92 » hz = 0.5*z; | 97 hz = 0.5 * z; |
93 » w = 1.0-hz; | 98 w = 1.0 - hz; |
94 » return w + (((1.0-w)-hz) + (z*r-x*y)); | 99 return w + (((1.0 - w) - hz) + (z * r - x * y)); |
95 } | 100 } |
96 #endif | 101 #endif |
OLD | NEW |