1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 Free Software Foundation, Inc. |
5 | * |
6 | * This program is free software; you can redistribute it and/or modify |
7 | * it under the terms of the GNU Lesser General Public License as published by |
8 | * the Free Software Foundation; either version 2.1 of the License, or |
9 | * (at your option) any later version. |
10 | * |
11 | * This program is distributed in the hope that it will be useful, |
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
14 | * GNU Lesser General Public License for more details. |
15 | * |
16 | * You should have received a copy of the GNU Lesser General Public License |
17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
18 | */ |
19 | /*************************************************************************/ |
20 | /* MODULE_NAME:mpexp.c */ |
21 | /* */ |
22 | /* FUNCTIONS: mpexp */ |
23 | /* */ |
24 | /* FILES NEEDED: mpa.h endian.h mpexp.h */ |
25 | /* mpa.c */ |
26 | /* */ |
27 | /* Multi-Precision exponential function subroutine */ |
28 | /* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */ |
29 | /*************************************************************************/ |
30 | |
31 | #include "endian.h" |
32 | #include "mpa.h" |
33 | #include <assert.h> |
34 | |
35 | #ifndef SECTION |
36 | # define SECTION |
37 | #endif |
38 | |
39 | /* Multi-Precision exponential function subroutine (for p >= 4, |
40 | 2**(-55) <= abs(x) <= 1024). */ |
41 | void |
42 | SECTION |
43 | __mpexp (mp_no *x, mp_no *y, int p) |
44 | { |
45 | int i, j, k, m, m1, m2, n; |
46 | mantissa_t b; |
47 | static const int np[33] = |
48 | { |
49 | 0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, |
50 | 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8 |
51 | }; |
52 | |
53 | static const int m1p[33] = |
54 | { |
55 | 0, 0, 0, 0, |
56 | 17, 23, 23, 28, |
57 | 27, 38, 42, 39, |
58 | 43, 47, 43, 47, |
59 | 50, 54, 57, 60, |
60 | 64, 67, 71, 74, |
61 | 68, 71, 74, 77, |
62 | 70, 73, 76, 78, |
63 | 81 |
64 | }; |
65 | static const int m1np[7][18] = |
66 | { |
67 | {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, |
68 | {0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, |
69 | {0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0}, |
70 | {0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0}, |
71 | {0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0}, |
72 | {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63}, |
73 | {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54} |
74 | }; |
75 | mp_no mps, mpk, mpt1, mpt2; |
76 | |
77 | /* Choose m,n and compute a=2**(-m). */ |
78 | n = np[p]; |
79 | m1 = m1p[p]; |
80 | b = X[1]; |
81 | m2 = 24 * EX; |
82 | for (; b < HALFRAD; m2--) |
83 | b *= 2; |
84 | if (b == HALFRAD) |
85 | { |
86 | for (i = 2; i <= p; i++) |
87 | { |
88 | if (X[i] != 0) |
89 | break; |
90 | } |
91 | if (i == p + 1) |
92 | m2--; |
93 | } |
94 | |
95 | m = m1 + m2; |
96 | if (__glibc_unlikely (m <= 0)) |
97 | { |
98 | /* The m1np array which is used to determine if we can reduce the |
99 | polynomial expansion iterations, has only 18 elements. Besides, |
100 | numbers smaller than those required by p >= 18 should not come here |
101 | at all since the fast phase of exp returns 1.0 for anything less |
102 | than 2^-55. */ |
103 | assert (p < 18); |
104 | m = 0; |
105 | for (i = n - 1; i > 0; i--, n--) |
106 | if (m1np[i][p] + m2 > 0) |
107 | break; |
108 | } |
109 | |
110 | /* Compute s=x*2**(-m). Put result in mps. This is the range-reduced input |
111 | that we will use to compute e^s. For the final result, simply raise it |
112 | to 2^m. */ |
113 | __pow_mp (-m, &mpt1, p); |
114 | __mul (x, &mpt1, &mps, p); |
115 | |
116 | /* Compute the Taylor series for e^s: |
117 | |
118 | 1 + x/1! + x^2/2! + x^3/3! ... |
119 | |
120 | for N iterations. We compute this as: |
121 | |
122 | e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n! |
123 | = 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n! |
124 | |
125 | k! is computed on the fly as KF and at the end of the polynomial loop, KF |
126 | is n!, which can be used directly. */ |
127 | __cpy (&mps, &mpt2, p); |
128 | |
129 | double kf = 1.0; |
130 | |
131 | /* Evaluate the rest. The result will be in mpt2. */ |
132 | for (k = n - 1; k > 0; k--) |
133 | { |
134 | /* n! / k! = n * (n - 1) ... * (n - k + 1) */ |
135 | kf *= k + 1; |
136 | |
137 | __dbl_mp (kf, &mpk, p); |
138 | __add (&mpt2, &mpk, &mpt1, p); |
139 | __mul (&mps, &mpt1, &mpt2, p); |
140 | } |
141 | __dbl_mp (kf, &mpk, p); |
142 | __dvd (&mpt2, &mpk, &mpt1, p); |
143 | __add (&__mpone, &mpt1, &mpt2, p); |
144 | |
145 | /* Raise polynomial value to the power of 2**m. Put result in y. */ |
146 | for (k = 0, j = 0; k < m;) |
147 | { |
148 | __sqr (&mpt2, &mpt1, p); |
149 | k++; |
150 | if (k == m) |
151 | { |
152 | j = 1; |
153 | break; |
154 | } |
155 | __sqr (&mpt1, &mpt2, p); |
156 | k++; |
157 | } |
158 | if (j) |
159 | __cpy (&mpt1, y, p); |
160 | else |
161 | __cpy (&mpt2, y, p); |
162 | return; |
163 | } |
164 | |