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 |