about summary refs log tree commit diff
path: root/sysdeps/ieee754/ldbl-128ibm/e_expl.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/e_expl.c')
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/e_expl.c257
1 files changed, 0 insertions, 257 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_expl.c b/sysdeps/ieee754/ldbl-128ibm/e_expl.c
deleted file mode 100644
index 3c4088f75f..0000000000
--- a/sysdeps/ieee754/ldbl-128ibm/e_expl.c
+++ /dev/null
@@ -1,257 +0,0 @@
-/* Quad-precision floating point e^x.
-   Copyright (C) 1999,2004,2006 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Jakub Jelinek <jj@ultra.linux.cz>
-   Partly based on double-precision code
-   by Geoffrey Keating <geoffk@ozemail.com.au>
-
-   The GNU C Library 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.
-
-   The GNU C Library 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 the GNU C Library; if not, write to the Free
-   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
-   02111-1307 USA.  */
-
-/* The basic design here is from
-   Abraham Ziv, "Fast Evaluation of Elementary Mathematical Functions with
-   Correctly Rounded Last Bit", ACM Trans. Math. Soft., 17 (3), September 1991,
-   pp. 410-423.
-
-   We work with number pairs where the first number is the high part and
-   the second one is the low part. Arithmetic with the high part numbers must
-   be exact, without any roundoff errors.
-
-   The input value, X, is written as
-   X = n * ln(2)_0 + arg1[t1]_0 + arg2[t2]_0 + x
-       - n * ln(2)_1 + arg1[t1]_1 + arg2[t2]_1 + xl
-
-   where:
-   - n is an integer, 16384 >= n >= -16495;
-   - ln(2)_0 is the first 93 bits of ln(2), and |ln(2)_0-ln(2)-ln(2)_1| < 2^-205
-   - t1 is an integer, 89 >= t1 >= -89
-   - t2 is an integer, 65 >= t2 >= -65
-   - |arg1[t1]-t1/256.0| < 2^-53
-   - |arg2[t2]-t2/32768.0| < 2^-53
-   - x + xl is whatever is left, |x + xl| < 2^-16 + 2^-53
-
-   Then e^x is approximated as
-
-   e^x = 2^n_1 ( 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
-	       + 2^n_0 e^(arg1[t1]_0 + arg1[t1]_1) e^(arg2[t2]_0 + arg2[t2]_1)
-		 * p (x + xl + n * ln(2)_1))
-   where:
-   - p(x) is a polynomial approximating e(x)-1
-   - e^(arg1[t1]_0 + arg1[t1]_1) is obtained from a table
-   - e^(arg2[t2]_0 + arg2[t2]_1) likewise
-   - n_1 + n_0 = n, so that |n_0| < -LDBL_MIN_EXP-1.
-
-   If it happens that n_1 == 0 (this is the usual case), that multiplication
-   is omitted.
-   */
-
-#ifndef _GNU_SOURCE
-#define _GNU_SOURCE
-#endif
-#include <float.h>
-#include <ieee754.h>
-#include <math.h>
-#include <fenv.h>
-#include <inttypes.h>
-#include <math_private.h>
-#include <sysdeps/ieee754/ldbl-128/t_expl.h>
-
-static const long double C[] = {
-/* Smallest integer x for which e^x overflows.  */
-#define himark C[0]
- 709.08956571282405153382846025171462914L,
-
-/* Largest integer x for which e^x underflows.  */
-#define lomark C[1]
--709.08956571282405153382846025171462914L,
-
-/* 3x2^96 */
-#define THREEp96 C[2]
- 59421121885698253195157962752.0L,
-
-/* 3x2^103 */
-#define THREEp103 C[3]
- 30423614405477505635920876929024.0L,
-
-/* 3x2^111 */
-#define THREEp111 C[4]
- 7788445287802241442795744493830144.0L,
-
-/* 1/ln(2) */
-#define M_1_LN2 C[5]
- 1.44269504088896340735992468100189204L,
-
-/* first 93 bits of ln(2) */
-#define M_LN2_0 C[6]
- 0.693147180559945309417232121457981864L,
-
-/* ln2_0 - ln(2) */
-#define M_LN2_1 C[7]
--1.94704509238074995158795957333327386E-31L,
-
-/* very small number */
-#define TINY C[8]
- 1.0e-308L,
-
-/* 2^16383 */
-#define TWO1023 C[9]
- 8.988465674311579538646525953945123668E+307L,
-
-/* 256 */
-#define TWO8 C[10]
- 256.0L,
-
-/* 32768 */
-#define TWO15 C[11]
- 32768.0L,
-
-/* Chebyshev polynom coeficients for (exp(x)-1)/x */
-#define P1 C[12]
-#define P2 C[13]
-#define P3 C[14]
-#define P4 C[15]
-#define P5 C[16]
-#define P6 C[17]
- 0.5L,
- 1.66666666666666666666666666666666683E-01L,
- 4.16666666666666666666654902320001674E-02L,
- 8.33333333333333333333314659767198461E-03L,
- 1.38888888889899438565058018857254025E-03L,
- 1.98412698413981650382436541785404286E-04L,
-};
-
-long double
-__ieee754_expl (long double x)
-{
-  /* Check for usual case.  */
-  if (isless (x, himark) && isgreater (x, lomark))
-    {
-      int tval1, tval2, unsafe, n_i, exponent2;
-      long double x22, n, result, xl;
-      union ibm_extended_long_double ex2_u, scale_u;
-      fenv_t oldenv;
-
-      feholdexcept (&oldenv);
-#ifdef FE_TONEAREST
-      fesetround (FE_TONEAREST);
-#endif
-
-      n = roundl(x*M_1_LN2);
-      x = x-n*M_LN2_0;
-      xl = n*M_LN2_1;
-
-      tval1 = roundl(x*TWO8);
-      x -= __expl_table[T_EXPL_ARG1+2*tval1];
-      xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
-
-      tval2 = roundl(x*TWO15);
-      x -= __expl_table[T_EXPL_ARG2+2*tval2];
-      xl -= __expl_table[T_EXPL_ARG2+2*tval2+1];
-
-      x = x + xl;
-
-      /* Compute ex2 = 2^n_0 e^(argtable[tval1]) e^(argtable[tval2]).  */
-      ex2_u.d = __expl_table[T_EXPL_RES1 + tval1]
-		* __expl_table[T_EXPL_RES2 + tval2];
-      n_i = (int)n;
-      /* 'unsafe' is 1 iff n_1 != 0.  */
-      unsafe = fabsl(n_i) >= -LDBL_MIN_EXP - 1;
-      ex2_u.ieee.exponent += n_i >> unsafe;
-      /* Fortunately, there are no subnormal lowpart doubles in
-	 __expl_table, only normal values and zeros.
-	 But after scaling it can be subnormal.  */
-      exponent2 = ex2_u.ieee.exponent2 + (n_i >> unsafe);
-      if (ex2_u.ieee.exponent2 == 0)
-	/* assert ((ex2_u.ieee.mantissa2|ex2_u.ieee.mantissa3) == 0) */;
-      else if (exponent2 > 0)
-	ex2_u.ieee.exponent2 = exponent2;
-      else if (exponent2 <= -54)
-	{
-	  ex2_u.ieee.exponent2 = 0;
-	  ex2_u.ieee.mantissa2 = 0;
-	  ex2_u.ieee.mantissa3 = 0;
-	}
-      else
-	{
-	  static const double
-	    two54 = 1.80143985094819840000e+16, /* 4350000000000000 */
-	    twom54 = 5.55111512312578270212e-17; /* 3C90000000000000 */
-	  ex2_u.dd[1] *= two54;
-	  ex2_u.ieee.exponent2 += n_i >> unsafe;
-	  ex2_u.dd[1] *= twom54;
-	}
-
-      /* Compute scale = 2^n_1.  */
-      scale_u.d = 1.0L;
-      scale_u.ieee.exponent += n_i - (n_i >> unsafe);
-
-      /* Approximate e^x2 - 1, using a seventh-degree polynomial,
-	 with maximum error in [-2^-16-2^-53,2^-16+2^-53]
-	 less than 4.8e-39.  */
-      x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
-
-      /* Return result.  */
-      fesetenv (&oldenv);
-
-      result = x22 * ex2_u.d + ex2_u.d;
-
-      /* Now we can test whether the result is ultimate or if we are unsure.
-	 In the later case we should probably call a mpn based routine to give
-	 the ultimate result.
-	 Empirically, this routine is already ultimate in about 99.9986% of
-	 cases, the test below for the round to nearest case will be false
-	 in ~ 99.9963% of cases.
-	 Without proc2 routine maximum error which has been seen is
-	 0.5000262 ulp.
-
-	  union ieee854_long_double ex3_u;
-
-	  #ifdef FE_TONEAREST
-	    fesetround (FE_TONEAREST);
-	  #endif
-	  ex3_u.d = (result - ex2_u.d) - x22 * ex2_u.d;
-	  ex2_u.d = result;
-	  ex3_u.ieee.exponent += LDBL_MANT_DIG + 15 + IEEE854_LONG_DOUBLE_BIAS
-	  			 - ex2_u.ieee.exponent;
-	  n_i = abs (ex3_u.d);
-	  n_i = (n_i + 1) / 2;
-	  fesetenv (&oldenv);
-	  #ifdef FE_TONEAREST
-	  if (fegetround () == FE_TONEAREST)
-	    n_i -= 0x4000;
-	  #endif
-	  if (!n_i) {
-	    return __ieee754_expl_proc2 (origx);
-	  }
-       */
-      if (!unsafe)
-	return result;
-      else
-	return result * scale_u.d;
-    }
-  /* Exceptional cases:  */
-  else if (isless (x, himark))
-    {
-      if (__isinfl (x))
-	/* e^-inf == 0, with no error.  */
-	return 0;
-      else
-	/* Underflow */
-	return TINY * TINY;
-    }
-  else
-    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
-    return TWO1023*x;
-}