@@ -124,7 +124,7 @@ type-ldouble-yes := ldouble
# double support
type-double-suffix :=
type-double-routines := branred doasin dosincos mpa mpatan2 \
- mpatan mpexp mplog mpsqrt mptan sincos32 \
+ mpatan mpsqrt mptan sincos32 \
sincostab k_rem_pio2
# float support
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
@@ -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 */
@@ -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;
-}
deleted file mode 100644
@@ -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;
-}
deleted file mode 100644
@@ -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);
- }
-}
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
deleted file mode 100644
@@ -1 +0,0 @@
-/* Not needed. */
@@ -10,9 +10,8 @@ libm-sysdep_routines += s_ceil-sse4_1 s_ceilf-sse4_1 s_floor-sse4_1 \
libm-sysdep_routines += e_exp-fma e_log-fma e_pow-fma s_atan-fma \
e_asin-fma e_atan2-fma s_sin-fma s_tan-fma \
- mplog-fma mpa-fma \
+ mpa-fma \
sincos32-fma doasin-fma dosincos-fma \
- mpexp-fma \
mpatan2-fma mpatan-fma mpsqrt-fma mptan-fma
CFLAGS-doasin-fma.c = -mfma -mavx2
@@ -25,8 +24,6 @@ CFLAGS-e_pow-fma.c = -mfma -mavx2 $(config-cflags-nofma)
CFLAGS-mpa-fma.c = -mfma -mavx2
CFLAGS-mpatan-fma.c = -mfma -mavx2
CFLAGS-mpatan2-fma.c = -mfma -mavx2
-CFLAGS-mpexp-fma.c = -mfma -mavx2
-CFLAGS-mplog-fma.c = -mfma -mavx2
CFLAGS-mpsqrt-fma.c = -mfma -mavx2
CFLAGS-mptan-fma.c = -mfma -mavx2
CFLAGS-s_atan-fma.c = -mfma -mavx2
@@ -50,9 +47,8 @@ CFLAGS-s_sincosf-fma.c = -mfma -mavx2
libm-sysdep_routines += e_exp-fma4 e_log-fma4 e_pow-fma4 s_atan-fma4 \
e_asin-fma4 e_atan2-fma4 s_sin-fma4 s_tan-fma4 \
- mplog-fma4 mpa-fma4 \
+ mpa-fma4 \
sincos32-fma4 doasin-fma4 dosincos-fma4 \
- mpexp-fma4 \
mpatan2-fma4 mpatan-fma4 mpsqrt-fma4 mptan-fma4
CFLAGS-doasin-fma4.c = -mfma4
@@ -65,8 +61,6 @@ CFLAGS-e_pow-fma4.c = -mfma4 $(config-cflags-nofma)
CFLAGS-mpa-fma4.c = -mfma4
CFLAGS-mpatan-fma4.c = -mfma4
CFLAGS-mpatan2-fma4.c = -mfma4
-CFLAGS-mpexp-fma4.c = -mfma4
-CFLAGS-mplog-fma4.c = -mfma4
CFLAGS-mpsqrt-fma4.c = -mfma4
CFLAGS-mptan-fma4.c = -mfma4
CFLAGS-s_atan-fma4.c = -mfma4
@@ -76,15 +70,12 @@ CFLAGS-s_tan-fma4.c = -mfma4
libm-sysdep_routines += e_exp-avx e_log-avx s_atan-avx \
e_atan2-avx s_sin-avx s_tan-avx \
- mplog-avx mpa-avx \
- mpexp-avx
+ mpa-avx
CFLAGS-e_atan2-avx.c = -msse2avx -DSSE2AVX
CFLAGS-e_exp-avx.c = -msse2avx -DSSE2AVX
CFLAGS-e_log-avx.c = -msse2avx -DSSE2AVX
CFLAGS-mpa-avx.c = -msse2avx -DSSE2AVX
-CFLAGS-mpexp-avx.c = -msse2avx -DSSE2AVX
-CFLAGS-mplog-avx.c = -msse2avx -DSSE2AVX
CFLAGS-s_atan-avx.c = -msse2avx -DSSE2AVX
CFLAGS-s_sin-avx.c = -msse2avx -DSSE2AVX
CFLAGS-s_tan-avx.c = -msse2avx -DSSE2AVX
@@ -1,8 +1,4 @@
#define __ieee754_log __ieee754_log_avx
-#define __mplog __mplog_avx
-#define __add __add_avx
-#define __dbl_mp __dbl_mp_avx
-#define __sub __sub_avx
#define SECTION __attribute__ ((section (".text.avx")))
#include <sysdeps/ieee754/dbl-64/e_log.c>
@@ -1,8 +1,4 @@
#define __ieee754_log __ieee754_log_fma
-#define __mplog __mplog_fma
-#define __add __add_fma
-#define __dbl_mp __dbl_mp_fma
-#define __sub __sub_fma
#define SECTION __attribute__ ((section (".text.fma")))
#include <sysdeps/ieee754/dbl-64/e_log.c>
@@ -1,8 +1,4 @@
#define __ieee754_log __ieee754_log_fma4
-#define __mplog __mplog_fma4
-#define __add __add_fma4
-#define __dbl_mp __dbl_mp_fma4
-#define __sub __sub_fma4
#define SECTION __attribute__ ((section (".text.fma4")))
#include <sysdeps/ieee754/dbl-64/e_log.c>
deleted file mode 100644
@@ -1,9 +0,0 @@
-#define __mpexp __mpexp_avx
-#define __add __add_avx
-#define __dbl_mp __dbl_mp_avx
-#define __dvd __dvd_avx
-#define __mul __mul_avx
-#define AVOID_MPEXP_H 1
-#define SECTION __attribute__ ((section (".text.avx")))
-
-#include <sysdeps/ieee754/dbl-64/mpexp.c>
deleted file mode 100644
@@ -1,9 +0,0 @@
-#define __mpexp __mpexp_fma
-#define __add __add_fma
-#define __dbl_mp __dbl_mp_fma
-#define __dvd __dvd_fma
-#define __mul __mul_fma
-#define AVOID_MPEXP_H 1
-#define SECTION __attribute__ ((section (".text.fma")))
-
-#include <sysdeps/ieee754/dbl-64/mpexp.c>
deleted file mode 100644
@@ -1,9 +0,0 @@
-#define __mpexp __mpexp_fma4
-#define __add __add_fma4
-#define __dbl_mp __dbl_mp_fma4
-#define __dvd __dvd_fma4
-#define __mul __mul_fma4
-#define AVOID_MPEXP_H 1
-#define SECTION __attribute__ ((section (".text.fma4")))
-
-#include <sysdeps/ieee754/dbl-64/mpexp.c>
deleted file mode 100644
@@ -1,8 +0,0 @@
-#define __mplog __mplog_avx
-#define __add __add_avx
-#define __mpexp __mpexp_avx
-#define __mul __mul_avx
-#define __sub __sub_avx
-#define SECTION __attribute__ ((section (".text.avx")))
-
-#include <sysdeps/ieee754/dbl-64/mplog.c>
deleted file mode 100644
@@ -1,8 +0,0 @@
-#define __mplog __mplog_fma
-#define __add __add_fma
-#define __mpexp __mpexp_fma
-#define __mul __mul_fma
-#define __sub __sub_fma
-#define SECTION __attribute__ ((section (".text.fma")))
-
-#include <sysdeps/ieee754/dbl-64/mplog.c>
deleted file mode 100644
@@ -1,8 +0,0 @@
-#define __mplog __mplog_fma4
-#define __add __add_fma4
-#define __mpexp __mpexp_fma4
-#define __mul __mul_fma4
-#define __sub __sub_fma4
-#define SECTION __attribute__ ((section (".text.fma4")))
-
-#include <sysdeps/ieee754/dbl-64/mplog.c>