| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */ |
| 2 /*- | 2 /*- |
| 3 * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> | 3 * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> |
| 4 * All rights reserved. | 4 * All rights reserved. |
| 5 * | 5 * |
| 6 * Redistribution and use in source and binary forms, with or without | 6 * Redistribution and use in source and binary forms, with or without |
| 7 * modification, are permitted provided that the following conditions | 7 * modification, are permitted provided that the following conditions |
| 8 * are met: | 8 * are met: |
| 9 * 1. Redistributions of source code must retain the above copyright | 9 * 1. Redistributions of source code must retain the above copyright |
| 10 * notice, this list of conditions and the following disclaimer. | 10 * notice, this list of conditions and the following disclaimer. |
| 11 * 2. Redistributions in binary form must reproduce the above copyright | 11 * 2. Redistributions in binary form must reproduce the above copyright |
| 12 * notice, this list of conditions and the following disclaimer in the | 12 * notice, this list of conditions and the following disclaimer in the |
| 13 * documentation and/or other materials provided with the distribution. | 13 * documentation and/or other materials provided with the distribution. |
| 14 * | 14 * |
| 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | 15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND |
| 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | 16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| 17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | 17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| 18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | 18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
| 19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | 19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
| 20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | 20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
| 21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | 21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
| 22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | 22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
| 23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | 23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
| 24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | 24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
| 25 * SUCH DAMAGE. | 25 * SUCH DAMAGE. |
| 26 */ | 26 */ |
| 27 | 27 |
| 28 #include "libm.h" | 28 #include "libm.h" |
| 29 | 29 |
| 30 static const uint32_t k = 1799; /* constant for reduction */ | 30 static const uint32_t k = 1799; /* constant for reduction */ |
| 31 static const double kln2 = 1246.97177782734161156; /* k * ln2 */ | 31 static const double kln2 = 1246.97177782734161156; /* k * ln2 */ |
| 32 | 32 |
| 33 /* | 33 /* |
| 34 * Compute exp(x), scaled to avoid spurious overflow. An exponent is | 34 * Compute exp(x), scaled to avoid spurious overflow. An exponent is |
| 35 * returned separately in 'expt'. | 35 * returned separately in 'expt'. |
| 36 * | 36 * |
| 37 * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 | 37 * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 |
| 38 * Output: 2**1023 <= y < 2**1024 | 38 * Output: 2**1023 <= y < 2**1024 |
| 39 */ | 39 */ |
| 40 static double __frexp_exp(double x, int *expt) | 40 static double __frexp_exp(double x, int* expt) { |
| 41 { | 41 double exp_x; |
| 42 » double exp_x; | 42 uint32_t hx; |
| 43 » uint32_t hx; | |
| 44 | 43 |
| 45 » /* | 44 /* |
| 46 » * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to | 45 * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to |
| 47 » * minimize |exp(kln2) - 2**k|. We also scale the exponent of | 46 * minimize |exp(kln2) - 2**k|. We also scale the exponent of |
| 48 » * exp_x to MAX_EXP so that the result can be multiplied by | 47 * exp_x to MAX_EXP so that the result can be multiplied by |
| 49 » * a tiny number without losing accuracy due to denormalization. | 48 * a tiny number without losing accuracy due to denormalization. |
| 50 » */ | 49 */ |
| 51 » exp_x = exp(x - kln2); | 50 exp_x = exp(x - kln2); |
| 52 » GET_HIGH_WORD(hx, exp_x); | 51 GET_HIGH_WORD(hx, exp_x); |
| 53 » *expt = (hx >> 20) - (0x3ff + 1023) + k; | 52 *expt = (hx >> 20) - (0x3ff + 1023) + k; |
| 54 » SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); | 53 SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); |
| 55 » return exp_x; | 54 return exp_x; |
| 56 } | 55 } |
| 57 | 56 |
| 58 /* | 57 /* |
| 59 * __ldexp_cexp(x, expt) compute exp(x) * 2**expt. | 58 * __ldexp_cexp(x, expt) compute exp(x) * 2**expt. |
| 60 * It is intended for large arguments (real part >= ln(DBL_MAX)) | 59 * It is intended for large arguments (real part >= ln(DBL_MAX)) |
| 61 * where care is needed to avoid overflow. | 60 * where care is needed to avoid overflow. |
| 62 * | 61 * |
| 63 * The present implementation is narrowly tailored for our hyperbolic and | 62 * The present implementation is narrowly tailored for our hyperbolic and |
| 64 * exponential functions. We assume expt is small (0 or -1), and the caller | 63 * exponential functions. We assume expt is small (0 or -1), and the caller |
| 65 * has filtered out very large x, for which overflow would be inevitable. | 64 * has filtered out very large x, for which overflow would be inevitable. |
| 66 */ | 65 */ |
| 67 double complex __ldexp_cexp(double complex z, int expt) | 66 double complex __ldexp_cexp(double complex z, int expt) { |
| 68 { | 67 double x, y, exp_x, scale1, scale2; |
| 69 » double x, y, exp_x, scale1, scale2; | 68 int ex_expt, half_expt; |
| 70 » int ex_expt, half_expt; | |
| 71 | 69 |
| 72 » x = creal(z); | 70 x = creal(z); |
| 73 » y = cimag(z); | 71 y = cimag(z); |
| 74 » exp_x = __frexp_exp(x, &ex_expt); | 72 exp_x = __frexp_exp(x, &ex_expt); |
| 75 » expt += ex_expt; | 73 expt += ex_expt; |
| 76 | 74 |
| 77 » /* | 75 /* |
| 78 » * Arrange so that scale1 * scale2 == 2**expt. We use this to | 76 * Arrange so that scale1 * scale2 == 2**expt. We use this to |
| 79 » * compensate for scalbn being horrendously slow. | 77 * compensate for scalbn being horrendously slow. |
| 80 » */ | 78 */ |
| 81 » half_expt = expt / 2; | 79 half_expt = expt / 2; |
| 82 » INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); | 80 INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); |
| 83 » half_expt = expt - half_expt; | 81 half_expt = expt - half_expt; |
| 84 » INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); | 82 INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); |
| 85 | 83 |
| 86 » return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 *
scale2); | 84 return CMPLX(cos(y) * exp_x * scale1 * scale2, |
| 85 sin(y) * exp_x * scale1 * scale2); |
| 87 } | 86 } |
| OLD | NEW |