diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_exp.c | 1 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpa.h | 31 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mpexp.c | 163 | ||||
| -rw-r--r-- | sysdeps/ieee754/dbl-64/mplog.c | 65 |
4 files changed, 0 insertions, 260 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index 1bcb382c77..62035a890b 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -23,7 +23,6 @@ /* exp1 */ /* */ /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */ -/* mpa.c mpexp.x */ /* */ /* An ultimate exp routine. Given an IEEE double machine number x */ /* it computes an almost correctly rounded (to nearest) value of e^x */ diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index f21b66d0f3..1e188de4d1 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -119,36 +119,5 @@ void __dvd (const mp_no *, const mp_no *, mp_no *, int); extern void __mpatan (mp_no *, mp_no *, int); extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int); extern void __mpsqrt (mp_no *, mp_no *, int); -extern void __mpexp (mp_no *, mp_no *, int); extern void __c32 (mp_no *, mp_no *, mp_no *, int); extern int __mpranred (double, mp_no *, int); - -/* Given a power POW, build a multiprecision number 2^POW. */ -static inline void -__pow_mp (int pow, mp_no *y, int p) -{ - int i, rem; - - /* The exponent is E such that E is a factor of 2^24. The remainder (of the - form 2^x) goes entirely into the first digit of the mantissa as it is - always less than 2^24. */ - EY = pow / 24; - rem = pow - EY * 24; - EY++; - - /* If the remainder is negative, it means that POW was negative since - |EY * 24| <= |pow|. Adjust so that REM is positive and still less than - 24 because of which, the mantissa digit is less than 2^24. */ - if (rem < 0) - { - EY--; - rem += 24; - } - /* The sign of any 2^x is always positive. */ - Y[0] = 1; - Y[1] = 1 << rem; - - /* Everything else is 0. */ - for (i = 2; i <= p; i++) - Y[i] = 0; -} diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c deleted file mode 100644 index e5228e55b3..0000000000 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ /dev/null @@ -1,163 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2018 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; -} diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c deleted file mode 100644 index f6cd5f4aa1..0000000000 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ /dev/null @@ -1,65 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2018 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:mplog.c */ -/* */ -/* FUNCTIONS: mplog */ -/* */ -/* FILES NEEDED: endian.h mpa.h mplog.h */ -/* mpexp.c */ -/* */ -/* Multi-Precision logarithm function subroutine (for precision p >= 4, */ -/* 2**(-1024) < x < 2**1024) and x is outside of the interval */ -/* [1-2**(-54),1+2**(-54)]. Upon entry, x should be set to the */ -/* multi-precision value of the input and y should be set into a multi- */ -/* precision value of an approximation of log(x) with relative error */ -/* bound of at most 2**(-52). The routine improves the accuracy of y. */ -/* */ -/************************************************************************/ -#include "endian.h" -#include "mpa.h" - -void -__mplog (mp_no *x, mp_no *y, int p) -{ - int i, m; - static const int mp[33] = - { - 0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 - }; - mp_no mpt1, mpt2; - - /* Choose m. */ - m = mp[p]; - - /* Perform m newton iterations to solve for y: exp(y) - x = 0. The - iterations formula is: y(n + 1) = y(n) + (x * exp(-y(n)) - 1). */ - __cpy (y, &mpt1, p); - for (i = 0; i < m; i++) - { - mpt1.d[0] = -mpt1.d[0]; - __mpexp (&mpt1, &mpt2, p); - __mul (x, &mpt2, &mpt1, p); - __sub (&mpt1, &__mpone, &mpt2, p); - __add (y, &mpt2, &mpt1, p); - __cpy (&mpt1, y, p); - } -} |