diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/dbl-64/mpexp.c')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/dbl-64/mpexp.c | 163 |
1 files changed, 163 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/dbl-64/mpexp.c b/REORG.TODO/sysdeps/ieee754/dbl-64/mpexp.c new file mode 100644 index 0000000000..e08f424133 --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/dbl-64/mpexp.c @@ -0,0 +1,163 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001-2017 Free Software Foundation, Inc. + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, see <http://www.gnu.org/licenses/>. + */ +/*************************************************************************/ +/* MODULE_NAME:mpexp.c */ +/* */ +/* FUNCTIONS: mpexp */ +/* */ +/* FILES NEEDED: mpa.h endian.h mpexp.h */ +/* mpa.c */ +/* */ +/* Multi-Precision exponential function subroutine */ +/* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */ +/*************************************************************************/ + +#include "endian.h" +#include "mpa.h" +#include <assert.h> + +#ifndef SECTION +# define SECTION +#endif + +/* Multi-Precision exponential function subroutine (for p >= 4, + 2**(-55) <= abs(x) <= 1024). */ +void +SECTION +__mpexp (mp_no *x, mp_no *y, int p) +{ + int i, j, k, m, m1, m2, n; + mantissa_t b; + static const int np[33] = + { + 0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, + 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8 + }; + + static const int m1p[33] = + { + 0, 0, 0, 0, + 17, 23, 23, 28, + 27, 38, 42, 39, + 43, 47, 43, 47, + 50, 54, 57, 60, + 64, 67, 71, 74, + 68, 71, 74, 77, + 70, 73, 76, 78, + 81 + }; + static const int m1np[7][18] = + { + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0}, + {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54} + }; + mp_no mps, mpk, mpt1, mpt2; + + /* Choose m,n and compute a=2**(-m). */ + n = np[p]; + m1 = m1p[p]; + b = X[1]; + m2 = 24 * EX; + for (; b < HALFRAD; m2--) + b *= 2; + if (b == HALFRAD) + { + for (i = 2; i <= p; i++) + { + if (X[i] != 0) + break; + } + if (i == p + 1) + m2--; + } + + m = m1 + m2; + if (__glibc_unlikely (m <= 0)) + { + /* The m1np array which is used to determine if we can reduce the + polynomial expansion iterations, has only 18 elements. Besides, + numbers smaller than those required by p >= 18 should not come here + at all since the fast phase of exp returns 1.0 for anything less + than 2^-55. */ + assert (p < 18); + m = 0; + for (i = n - 1; i > 0; i--, n--) + if (m1np[i][p] + m2 > 0) + break; + } + + /* Compute s=x*2**(-m). Put result in mps. This is the range-reduced input + that we will use to compute e^s. For the final result, simply raise it + to 2^m. */ + __pow_mp (-m, &mpt1, p); + __mul (x, &mpt1, &mps, p); + + /* Compute the Taylor series for e^s: + + 1 + x/1! + x^2/2! + x^3/3! ... + + for N iterations. We compute this as: + + e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n! + = 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n! + + k! is computed on the fly as KF and at the end of the polynomial loop, KF + is n!, which can be used directly. */ + __cpy (&mps, &mpt2, p); + + double kf = 1.0; + + /* Evaluate the rest. The result will be in mpt2. */ + for (k = n - 1; k > 0; k--) + { + /* n! / k! = n * (n - 1) ... * (n - k + 1) */ + kf *= k + 1; + + __dbl_mp (kf, &mpk, p); + __add (&mpt2, &mpk, &mpt1, p); + __mul (&mps, &mpt1, &mpt2, p); + } + __dbl_mp (kf, &mpk, p); + __dvd (&mpt2, &mpk, &mpt1, p); + __add (&__mpone, &mpt1, &mpt2, p); + + /* Raise polynomial value to the power of 2**m. Put result in y. */ + for (k = 0, j = 0; k < m;) + { + __sqr (&mpt2, &mpt1, p); + k++; + if (k == m) + { + j = 1; + break; + } + __sqr (&mpt1, &mpt2, p); + k++; + } + if (j) + __cpy (&mpt1, y, p); + else + __cpy (&mpt2, y, p); + return; +} |