2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2017 Free Software Foundation, Inc.
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.
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.
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/>.
19 /*************************************************************************/
20 /* MODULE_NAME:mpexp.c */
22 /* FUNCTIONS: mpexp */
24 /* FILES NEEDED: mpa.h endian.h mpexp.h */
27 /* Multi-Precision exponential function subroutine */
28 /* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */
29 /*************************************************************************/
39 /* Multi-Precision exponential function subroutine (for p >= 4,
40 2**(-55) <= abs(x) <= 1024). */
43 __mpexp (mp_no
*x
, mp_no
*y
, int p
)
45 int i
, j
, k
, m
, m1
, m2
, n
;
47 static const int np
[33] =
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
53 static const int m1p
[33] =
65 static const int m1np
[7][18] =
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}
75 mp_no mps
, mpk
, mpt1
, mpt2
;
77 /* Choose m,n and compute a=2**(-m). */
82 for (; b
< HALFRAD
; m2
--)
86 for (i
= 2; i
<= p
; i
++)
96 if (__glibc_unlikely (m
<= 0))
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
105 for (i
= n
- 1; i
> 0; i
--, n
--)
106 if (m1np
[i
][p
] + m2
> 0)
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
113 __pow_mp (-m
, &mpt1
, p
);
114 __mul (x
, &mpt1
, &mps
, p
);
116 /* Compute the Taylor series for e^s:
118 1 + x/1! + x^2/2! + x^3/3! ...
120 for N iterations. We compute this as:
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!
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
);
131 /* Evaluate the rest. The result will be in mpt2. */
132 for (k
= n
- 1; k
> 0; k
--)
134 /* n! / k! = n * (n - 1) ... * (n - k + 1) */
137 __dbl_mp (kf
, &mpk
, p
);
138 __add (&mpt2
, &mpk
, &mpt1
, p
);
139 __mul (&mps
, &mpt1
, &mpt2
, p
);
141 __dbl_mp (kf
, &mpk
, p
);
142 __dvd (&mpt2
, &mpk
, &mpt1
, p
);
143 __add (&__mpone
, &mpt1
, &mpt2
, p
);
145 /* Raise polynomial value to the power of 2**m. Put result in y. */
146 for (k
= 0, j
= 0; k
< m
;)
148 __sqr (&mpt2
, &mpt1
, p
);
155 __sqr (&mpt1
, &mpt2
, p
);