OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/k_expf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/k_expf.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 = 235; /* constant for reduction */ | 30 static const uint32_t k = 235; /* constant for reduction */ |
31 static const float kln2 = 162.88958740F; /* k * ln2 */ | 31 static const float kln2 = 162.88958740F; /* k * ln2 */ |
32 | 32 |
33 /* | 33 /* |
34 * See __cexp.c for details. | 34 * See __cexp.c for details. |
35 * | 35 * |
36 * Input: ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7 | 36 * Input: ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7 |
37 * Output: 2**127 <= y < 2**128 | 37 * Output: 2**127 <= y < 2**128 |
38 */ | 38 */ |
39 static float __frexp_expf(float x, int *expt) | 39 static float __frexp_expf(float x, int* expt) { |
40 { | 40 float exp_x; |
41 » float exp_x; | 41 uint32_t hx; |
42 » uint32_t hx; | |
43 | 42 |
44 » exp_x = expf(x - kln2); | 43 exp_x = expf(x - kln2); |
45 » GET_FLOAT_WORD(hx, exp_x); | 44 GET_FLOAT_WORD(hx, exp_x); |
46 » *expt = (hx >> 23) - (0x7f + 127) + k; | 45 *expt = (hx >> 23) - (0x7f + 127) + k; |
47 » SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23)); | 46 SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23)); |
48 » return exp_x; | 47 return exp_x; |
49 } | 48 } |
50 | 49 |
51 float complex __ldexp_cexpf(float complex z, int expt) | 50 float complex __ldexp_cexpf(float complex z, int expt) { |
52 { | 51 float x, y, exp_x, scale1, scale2; |
53 » float x, y, exp_x, scale1, scale2; | 52 int ex_expt, half_expt; |
54 » int ex_expt, half_expt; | |
55 | 53 |
56 » x = crealf(z); | 54 x = crealf(z); |
57 » y = cimagf(z); | 55 y = cimagf(z); |
58 » exp_x = __frexp_expf(x, &ex_expt); | 56 exp_x = __frexp_expf(x, &ex_expt); |
59 » expt += ex_expt; | 57 expt += ex_expt; |
60 | 58 |
61 » half_expt = expt / 2; | 59 half_expt = expt / 2; |
62 » SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23); | 60 SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23); |
63 » half_expt = expt - half_expt; | 61 half_expt = expt - half_expt; |
64 » SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); | 62 SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); |
65 | 63 |
66 » return CMPLXF(cosf(y) * exp_x * scale1 * scale2, | 64 return CMPLXF(cosf(y) * exp_x * scale1 * scale2, |
67 » sinf(y) * exp_x * scale1 * scale2); | 65 sinf(y) * exp_x * scale1 * scale2); |
68 } | 66 } |
OLD | NEW |