|
16620
|
1 /* Exponential base 2 function. |
|
|
2 Copyright (C) 2011-2012 Free Software Foundation, Inc. |
|
|
3 |
|
|
4 This program is free software: you can redistribute it and/or modify |
|
|
5 it under the terms of the GNU General Public License as published by |
|
|
6 the Free Software Foundation; either version 3 of the License, or |
|
|
7 (at your option) any later version. |
|
|
8 |
|
|
9 This program is distributed in the hope that it will be useful, |
|
|
10 but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
12 GNU General Public License for more details. |
|
|
13 |
|
|
14 You should have received a copy of the GNU General Public License |
|
|
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */ |
|
|
16 |
|
|
17 #include <config.h> |
|
|
18 |
|
|
19 /* Specification. */ |
|
|
20 #include <math.h> |
|
|
21 |
|
|
22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE |
|
|
23 |
|
|
24 long double |
|
|
25 exp2l (long double x) |
|
|
26 { |
|
|
27 return exp2 (x); |
|
|
28 } |
|
|
29 |
|
|
30 #else |
|
|
31 |
|
|
32 # include <float.h> |
|
|
33 |
|
|
34 /* gl_expl_table[i] = exp((i - 128) * log(2)/256). */ |
|
|
35 extern const long double gl_expl_table[257]; |
|
|
36 |
|
|
37 /* Best possible approximation of log(2) as a 'long double'. */ |
|
|
38 #define LOG2 0.693147180559945309417232121458176568075L |
|
|
39 |
|
|
40 /* Best possible approximation of 1/log(2) as a 'long double'. */ |
|
|
41 #define LOG2_INVERSE 1.44269504088896340735992468100189213743L |
|
|
42 |
|
|
43 /* Best possible approximation of log(2)/256 as a 'long double'. */ |
|
|
44 #define LOG2_BY_256 0.00270760617406228636491106297444600221904L |
|
|
45 |
|
|
46 /* Best possible approximation of 256/log(2) as a 'long double'. */ |
|
|
47 #define LOG2_BY_256_INVERSE 369.329930467574632284140718336484387181L |
|
|
48 |
|
|
49 long double |
|
|
50 exp2l (long double x) |
|
|
51 { |
|
|
52 /* exp2(x) = exp(x*log(2)). |
|
|
53 If we would compute it like this, there would be rounding errors for |
|
|
54 integer or near-integer values of x. To avoid these, we inline the |
|
|
55 algorithm for exp(), and the multiplication with log(2) cancels a |
|
|
56 division by log(2). */ |
|
|
57 |
|
|
58 if (isnanl (x)) |
|
|
59 return x; |
|
|
60 |
|
|
61 if (x > (long double) LDBL_MAX_EXP) |
|
|
62 /* x > LDBL_MAX_EXP |
|
|
63 hence exp2(x) > 2^LDBL_MAX_EXP, overflows to Infinity. */ |
|
|
64 return HUGE_VALL; |
|
|
65 |
|
|
66 if (x < (long double) (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG)) |
|
|
67 /* x < (LDBL_MIN_EXP - 1 - LDBL_MANT_DIG) |
|
|
68 hence exp2(x) < 2^(LDBL_MIN_EXP-1-LDBL_MANT_DIG), |
|
|
69 underflows to zero. */ |
|
|
70 return 0.0L; |
|
|
71 |
|
|
72 /* Decompose x into |
|
|
73 x = n + m/256 + y/log(2) |
|
|
74 where |
|
|
75 n is an integer, |
|
|
76 m is an integer, -128 <= m <= 128, |
|
|
77 y is a number, |y| <= log(2)/512 + epsilon = 0.00135... |
|
|
78 Then |
|
|
79 exp2(x) = 2^n * exp(m * log(2)/256) * exp(y) |
|
|
80 The first factor is an ldexpl() call. |
|
|
81 The second factor is a table lookup. |
|
|
82 The third factor is computed |
|
|
83 - either as sinh(y) + cosh(y) |
|
|
84 where sinh(y) is computed through the power series: |
|
|
85 sinh(y) = y + y^3/3! + y^5/5! + ... |
|
|
86 and cosh(y) is computed as hypot(1, sinh(y)), |
|
|
87 - or as exp(2*z) = (1 + tanh(z)) / (1 - tanh(z)) |
|
|
88 where z = y/2 |
|
|
89 and tanh(z) is computed through its power series: |
|
|
90 tanh(z) = z |
|
|
91 - 1/3 * z^3 |
|
|
92 + 2/15 * z^5 |
|
|
93 - 17/315 * z^7 |
|
|
94 + 62/2835 * z^9 |
|
|
95 - 1382/155925 * z^11 |
|
|
96 + 21844/6081075 * z^13 |
|
|
97 - 929569/638512875 * z^15 |
|
|
98 + ... |
|
16651
|
99 Since |z| <= log(2)/1024 < 0.0007, the relative contribution of the |
|
|
100 z^13 term is < 0.0007^12 < 2^-120 <= 2^-LDBL_MANT_DIG, therefore we |
|
|
101 can truncate the series after the z^11 term. */ |
|
16620
|
102 |
|
|
103 { |
|
|
104 long double nm = roundl (x * 256.0L); /* = 256 * n + m */ |
|
|
105 long double z = (x * 256.0L - nm) * (LOG2_BY_256 * 0.5L); |
|
|
106 |
|
|
107 /* Coefficients of the power series for tanh(z). */ |
|
|
108 #define TANH_COEFF_1 1.0L |
|
|
109 #define TANH_COEFF_3 -0.333333333333333333333333333333333333334L |
|
|
110 #define TANH_COEFF_5 0.133333333333333333333333333333333333334L |
|
|
111 #define TANH_COEFF_7 -0.053968253968253968253968253968253968254L |
|
|
112 #define TANH_COEFF_9 0.0218694885361552028218694885361552028218L |
|
|
113 #define TANH_COEFF_11 -0.00886323552990219656886323552990219656886L |
|
|
114 #define TANH_COEFF_13 0.00359212803657248101692546136990581435026L |
|
|
115 #define TANH_COEFF_15 -0.00145583438705131826824948518070211191904L |
|
|
116 |
|
|
117 long double z2 = z * z; |
|
|
118 long double tanh_z = |
|
|
119 (((((TANH_COEFF_11 |
|
|
120 * z2 + TANH_COEFF_9) |
|
|
121 * z2 + TANH_COEFF_7) |
|
|
122 * z2 + TANH_COEFF_5) |
|
|
123 * z2 + TANH_COEFF_3) |
|
|
124 * z2 + TANH_COEFF_1) |
|
|
125 * z; |
|
|
126 |
|
|
127 long double exp_y = (1.0L + tanh_z) / (1.0L - tanh_z); |
|
|
128 |
|
|
129 int n = (int) roundl (nm * (1.0L / 256.0L)); |
|
|
130 int m = (int) nm - 256 * n; |
|
|
131 |
|
|
132 return ldexpl (gl_expl_table[128 + m] * exp_y, n); |
|
|
133 } |
|
|
134 } |
|
|
135 |
|
|
136 #endif |
