OLD | NEW |
1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */ | 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_expl.c */ |
2 /* | 2 /* |
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> |
4 * | 4 * |
5 * Permission to use, copy, modify, and distribute this software for any | 5 * Permission to use, copy, modify, and distribute this software for any |
6 * purpose with or without fee is hereby granted, provided that the above | 6 * purpose with or without fee is hereby granted, provided that the above |
7 * copyright notice and this permission notice appear in all copies. | 7 * copyright notice and this permission notice appear in all copies. |
8 * | 8 * |
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES |
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF |
(...skipping 50 matching lines...) Expand 10 before | Expand all | Expand 10 after Loading... |
61 * | 61 * |
62 * message condition value returned | 62 * message condition value returned |
63 * exp underflow x < MINLOG 0.0 | 63 * exp underflow x < MINLOG 0.0 |
64 * exp overflow x > MAXLOG MAXNUM | 64 * exp overflow x > MAXLOG MAXNUM |
65 * | 65 * |
66 */ | 66 */ |
67 | 67 |
68 #include "libm.h" | 68 #include "libm.h" |
69 | 69 |
70 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 70 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
71 long double expl(long double x) | 71 long double expl(long double x) { |
72 { | 72 return exp(x); |
73 » return exp(x); | |
74 } | 73 } |
75 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 | 74 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 |
76 | 75 |
77 static const long double P[3] = { | 76 static const long double P[3] = { |
78 1.2617719307481059087798E-4L, | 77 1.2617719307481059087798E-4L, 3.0299440770744196129956E-2L, |
79 3.0299440770744196129956E-2L, | 78 9.9999999999999999991025E-1L, |
80 9.9999999999999999991025E-1L, | |
81 }; | 79 }; |
82 static const long double Q[4] = { | 80 static const long double Q[4] = { |
83 3.0019850513866445504159E-6L, | 81 3.0019850513866445504159E-6L, 2.5244834034968410419224E-3L, |
84 2.5244834034968410419224E-3L, | 82 2.2726554820815502876593E-1L, 2.0000000000000000000897E0L, |
85 2.2726554820815502876593E-1L, | |
86 2.0000000000000000000897E0L, | |
87 }; | 83 }; |
88 static const long double | 84 static const long double LN2HI = 6.9314575195312500000000E-1L, |
89 LN2HI = 6.9314575195312500000000E-1L, | 85 LN2LO = 1.4286068203094172321215E-6L, |
90 LN2LO = 1.4286068203094172321215E-6L, | 86 LOG2E = 1.4426950408889634073599E0L; |
91 LOG2E = 1.4426950408889634073599E0L; | |
92 | 87 |
93 long double expl(long double x) | 88 long double expl(long double x) { |
94 { | 89 long double px, xx; |
95 » long double px, xx; | 90 int k; |
96 » int k; | |
97 | 91 |
98 » if (isnan(x)) | 92 if (isnan(x)) |
99 » » return x; | 93 return x; |
100 » if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */ | 94 if (x > 11356.5234062941439488L) /* x > ln(2^16384 - 0.5) */ |
101 » » return x * 0x1p16383L; | 95 return x * 0x1p16383L; |
102 » if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */ | 96 if (x < -11399.4985314888605581L) /* x < ln(2^-16446) */ |
103 » » return -0x1p-16445L/x; | 97 return -0x1p-16445L / x; |
104 | 98 |
105 » /* Express e**x = e**f 2**k | 99 /* Express e**x = e**f 2**k |
106 » * = e**(f + k ln(2)) | 100 * = e**(f + k ln(2)) |
107 » */ | 101 */ |
108 » px = floorl(LOG2E * x + 0.5); | 102 px = floorl(LOG2E * x + 0.5); |
109 » k = px; | 103 k = px; |
110 » x -= px * LN2HI; | 104 x -= px * LN2HI; |
111 » x -= px * LN2LO; | 105 x -= px * LN2LO; |
112 | 106 |
113 » /* rational approximation of the fractional part: | 107 /* rational approximation of the fractional part: |
114 » * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2)) | 108 * e**x = 1 + 2x P(x**2)/(Q(x**2) - x P(x**2)) |
115 » */ | 109 */ |
116 » xx = x * x; | 110 xx = x * x; |
117 » px = x * __polevll(xx, P, 2); | 111 px = x * __polevll(xx, P, 2); |
118 » x = px/(__polevll(xx, Q, 3) - px); | 112 x = px / (__polevll(xx, Q, 3) - px); |
119 » x = 1.0 + 2.0 * x; | 113 x = 1.0 + 2.0 * x; |
120 » return scalbnl(x, k); | 114 return scalbnl(x, k); |
121 } | 115 } |
122 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 | 116 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 |
123 // TODO: broken implementation to make things compile | 117 // TODO: broken implementation to make things compile |
124 long double expl(long double x) | 118 long double expl(long double x) { |
125 { | 119 return exp(x); |
126 » return exp(x); | |
127 } | 120 } |
128 #endif | 121 #endif |
OLD | NEW |