1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2017 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). */
41void
42SECTION
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