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