Update.

1999-08-06  Jakub Jelinek  <jj@ultra.linux.cz>

	* sysdeps/ieee754/ldbl-128/e_expl.c: New file.
	* sysdeps/ieee754/ldbl-128/t_expl.h: New file.
	* sysdeps/ieee754/ldbl-128/Dist: Add t_expl.h.

1 files changed, 248 insertions, 0 deletions

diff --git a/sysdeps/ieee754/ldbl-128/e_expl.c b/sysdeps/ieee754/ldbl-128/e_expl.c
new file mode 100644
index 0000000000..959d958736
--- /dev/null
+++ b/sysdeps/ieee754/ldbl-128/e_expl.c
@@ -0,0 +1,248 @@
+/* Quad-precision floating point e^x.
+   Copyright (C) 1999 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  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 "t_expl.h"
+
+static const long double C[] = {
+/* Smallest integer x for which e^x overflows.  */
+#define himark C[0]
+ 11356.523406294143949491931077970765L,
+ 
+/* Largest integer x for which e^x underflows.  */
+#define lomark C[1]
+-11433.4627433362978788372438434526231L,
+
+/* 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-4900L,
+
+/* 2^16383 */
+#define TWO16383 C[9]
+ 5.94865747678615882542879663314003565E+4931L,
+
+/* 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;
+      long double x22, n, t, result, xl;
+      union ieee854_long_double ex2_u, scale_u;
+      fenv_t oldenv;
+
+      feholdexcept (&oldenv);
+#ifdef FE_TONEAREST
+      fesetround (FE_TONEAREST);
+#endif
+
+      /* Calculate n.  */
+      n = x * M_1_LN2 + THREEp111;
+      n -= THREEp111;
+      x = x - n * M_LN2_0;
+      xl = n * M_LN2_1;
+
+      /* Calculate t/256.  */
+      t = x + THREEp103;
+      t -= THREEp103;
+
+      /* Compute tval1 = t.  */
+      tval1 = (int) (t * TWO8);
+
+      x -= __expl_table[T_EXPL_ARG1+2*tval1];
+      xl -= __expl_table[T_EXPL_ARG1+2*tval1+1];
+
+      /* Calculate t/32768.  */
+      t = x + THREEp96;
+      t -= THREEp96;
+
+      /* Compute tval2 = t.  */
+      tval2 = (int) (t * 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 = abs(n_i) >= -LDBL_MIN_EXP - 1;
+      ex2_u.ieee.exponent += n_i >> unsafe;
+
+      /* 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 TWO16383*x;
+}