File math-remove-slow-path.patch of Package glibc

Overview Repositories Revisions Requests Users Attributes Meta

File math-remove-slow-path.patch of Package glibc (Revision e53058db315273b86de539470702f9df)

Currently displaying revision e53058db315273b86de539470702f9df , Show latest

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Cleanup ifdefs.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sin.c (__sincos): Refactor using the same
	logic as sin and cos.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (do_sin): Use TAYLOR_SIN for small
	inputs.  Return correct sign.
	(do_sincos): Remove small input check before do_sin, let do_sin set
	the sign.
	(__sin): Likewise.
	(__cos): Likewise.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (TAYLOR_SLOW): Remove.
	(do_cos_slow): Likewise.
	(do_sin_slow): Likewise.
	(reduce_and_compute): Likewise.
	(slow): Likewise.
	(slow1): Likewise.
	(slow2): Likewise.
	(sloww): Likewise.
	(sloww1): Likewise.
	(sloww2): Likewise.
	(bslow): Likewise.
	(bslow1): Likewise.
	(bslow2): Likewise.
	(cslow2): Likewise.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (TAYLOR_SIN): Remove cor parameter.
	(do_cos): Remove corp parameter and calculations.
	(do_sin): Likewise.
	(do_sincos): Remove cor variable.
	(__sin): Use do_sincos for huge inputs.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.
	(reduce_and_compute_sincos): Remove unused function.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_1): Rename to
	reduce_sincos, improve accuracy to 136 bits.
	(do_sincos_1): Rename to do_sincos, remove fallbacks to slow functions.
	(__sin): Use improved reduction and simplified do_sincos calculation.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_2): Remove function.
	(do_sincos_2): Likewise.
	(__sin): Remove middle range reduction case.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Remove middle range
	reduction case.

2018-04-03  Wilco Dijkstra  <wdijkstr@arm.com>

* sysdeps/aarch64/libm-test-ulps: Update ULP for sin, cos, sincos.
	* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Remove slow paths for small
	inputs.
	(__cos): Likewise.
	* sysdeps/x86_64/fpu/libm-test-ulps: Update ULP for sin, cos, sincos.

2018-02-12  Szabolcs Nagy  <szabolcs.nagy@arm.com>

* manual/probes.texi: Remove slowexp probes.
	* math/Makefile: Remove slowexp.
	* sysdeps/generic/math_private.h (__slowexp): Remove.
	* sysdeps/ieee754/dbl-64/e_exp.c (__ieee754_exp): Remove __slowexp and
	document error bounds.
	* sysdeps/i386/fpu/slowexp.c: Remove.
	* sysdeps/ia64/fpu/slowexp.c: Remove.
	* sysdeps/ieee754/dbl-64/slowexp.c: Remove.
	* sysdeps/ieee754/dbl-64/uexp.h (err_0): Remove.
	* sysdeps/m68k/m680x0/fpu/slowexp.c: Remove.
	* sysdeps/powerpc/power4/fpu/Makefile (CPPFLAGS-slowexp.c): Remove.
	* sysdeps/x86_64/fpu/multiarch/Makefile: Remove slowexp-fma.
	* sysdeps/x86_64/fpu/multiarch/e_exp-avx.c (__slowexp): Remove.
	* sysdeps/x86_64/fpu/multiarch/e_exp-fma.c (__slowexp): Remove.
	* sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c (__slowexp): Remove.
	* sysdeps/x86_64/fpu/multiarch/slowexp-avx.c: Remove.
	* sysdeps/x86_64/fpu/multiarch/slowexp-fma.c: Remove.
	* sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: Remove.

2018-02-12  Wilco Dijkstra  <wdijkstr@arm.com>

[BZ #13932]
	* sysdeps/ieee754/dbl-64/uexp.h (err_1): Remove.
	* benchtests/pow-inputs: Update comment for slow path cases.
	* manual/probes.texi (slowpow_p10): Delete removed probe.
	(slowpow_p10): Likewise.
	* math/Makefile: Remove halfulp.c and slowpow.c.
	* sysdeps/aarch64/libm-test-ulps: Set ULP of pow to 1.
	* sysdeps/generic/math_private.h (__exp1): Remove error argument.
	(__halfulp): Remove.
	(__slowpow): Remove.
	* sysdeps/i386/fpu/halfulp.c: Delete file.
	* sysdeps/i386/fpu/slowpow.c: Likewise.
	* sysdeps/ia64/fpu/halfulp.c: Likewise.
	* sysdeps/ia64/fpu/slowpow.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_exp.c (__exp1): Remove error argument,
	improve comments and add error analysis.
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Add error analysis.
	(power1): Remove function:
	(log1): Remove error argument, add error analysis.
	(my_log2): Remove function.
	* sysdeps/ieee754/dbl-64/halfulp.c: Delete file.
	* sysdeps/ieee754/dbl-64/slowpow.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/halfulp.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/slowpow.c: Likewise.
	* sysdeps/powerpc/power4/fpu/Makefile: Remove CPPFLAGS-slowpow.c.
	* sysdeps/x86_64/fpu/libm-test-ulps: Set ULP of pow to 1.
	* sysdeps/x86_64/fpu/multiarch/Makefile: Remove slowpow-fma.c,
	slowpow-fma4.c, halfulp-fma.c, halfulp-fma4.c.
	* sysdeps/x86_64/fpu/multiarch/e_pow-fma.c (__slowpow): Remove define.
	* sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c (__slowpow): Likewise.
	* sysdeps/x86_64/fpu/multiarch/halfulp-fma.c: Delete file.
	* sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowpow-fma.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c: Likewise.

2018-02-07  Wilco Dijkstra  <wdijkstr@arm.com>

* manual/probes.texi (slowlog): Delete documentation of removed probe.
	(slowlog_inexact): Likewise
	* sysdeps/ieee754/dbl-64/e_log.c (__ieee754_log): Remove slow paths.
	* sysdeps/ieee754/dbl-64/ulog.h: Remove unused declarations.

Index: glibc-2.26/benchtests/pow-inputs
===================================================================
--- glibc-2.26.orig/benchtests/pow-inputs
+++ glibc-2.26/benchtests/pow-inputs
@@ -302,8 +302,7 @@
 0x1.c004d2256a5b8p402, -0x1.a01df480fdcb7p98
 0x1.52b9d41aaa1e9p-589, -0x1.292cb15f1459dp46
 -0x1.ea9ca6fa0919ep-279, -0x1.601e44b6a588cp40
-# pow slow path at 240 bits
-# Implemented in sysdeps/ieee754/dbl-64/slowpow.c
+# old pow slow path at 240 bits
 ## name: 240bits
 0x1.01fcd33493ea3p596, -0x1.724bd4e887783p-14
 0x1.032ff59ab34fdp-540, -0x1.61e3632080b87p-24
@@ -405,8 +404,7 @@
 0x1.fae913d4f952ep-809, -0x1.4b649402fce63p-6
 0x1.fe6d725408f24p484, -0x1.25f4f6441d2e4p-12
 0x1.ff6393f9150ccp-718, 0x1.a0cb50a9bf2f3p-31
-# pow slowest path at 768 bits
-# Implemented in sysdeps/ieee754/dbl-64/slowpow.c
+# old pow slowest path at 768 bits
 ## name: 768bits
 1.0000000000000020, 1.5
 0x1.006777b4b61dep843, -0x1.67e3145491872p-1
Index: glibc-2.26/manual/probes.texi
===================================================================
--- glibc-2.26.orig/manual/probes.texi
+++ glibc-2.26/manual/probes.texi
@@ -265,53 +265,6 @@ Unless explicitly mentioned otherwise, a
 precision in the mantissa of the multiple precision number.  Hence, a precision
 level of 32 implies 768 bits of precision in the mantissa.
 
-@deftp Probe slowexp_p6 (double @var{$arg1}, double @var{$arg2})
-This probe is triggered when the @code{exp} function is called with an
-input that results in multiple precision computation with precision
-6.  Argument @var{$arg1} is the input value and @var{$arg2} is the
-computed output.
-@end deftp
-
-@deftp Probe slowexp_p32 (double @var{$arg1}, double @var{$arg2})
-This probe is triggered when the @code{exp} function is called with an
-input that results in multiple precision computation with precision
-32.  Argument @var{$arg1} is the input value and @var{$arg2} is the
-computed output.
-@end deftp
-
-@deftp Probe slowpow_p10 (double @var{$arg1}, double @var{$arg2}, double @var{$arg3}, double @var{$arg4})
-This probe is triggered when the @code{pow} function is called with
-inputs that result in multiple precision computation with precision
-10.  Arguments @var{$arg1} and @var{$arg2} are the input values,
-@code{$arg3} is the value computed in the fast phase of the algorithm
-and @code{$arg4} is the final accurate value.
-@end deftp
-
-@deftp Probe slowpow_p32 (double @var{$arg1}, double @var{$arg2}, double @var{$arg3}, double @var{$arg4})
-This probe is triggered when the @code{pow} function is called with an
-input that results in multiple precision computation with precision
-32.  Arguments @var{$arg1} and @var{$arg2} are the input values,
-@code{$arg3} is the value computed in the fast phase of the algorithm
-and @code{$arg4} is the final accurate value.
-@end deftp
-
-@deftp Probe slowlog (int @var{$arg1}, double @var{$arg2}, double @var{$arg3})
-This probe is triggered when the @code{log} function is called with an
-input that results in multiple precision computation.  Argument
-@var{$arg1} is the precision with which the computation succeeded.
-Argument @var{$arg2} is the input and @var{$arg3} is the computed
-output.
-@end deftp
-
-@deftp Probe slowlog_inexact (int @var{$arg1}, double @var{$arg2}, double @var{$arg3})
-This probe is triggered when the @code{log} function is called with an
-input that results in multiple precision computation and none of the
-multiple precision computations result in an accurate result.
-Argument @var{$arg1} is the maximum precision with which computations
-were performed.  Argument @var{$arg2} is the input and @var{$arg3} is
-the computed output.
-@end deftp
-
 @deftp Probe slowatan2 (int @var{$arg1}, double @var{$arg2}, double @var{$arg3}, double @var{$arg4})
 This probe is triggered when the @code{atan2} function is called with
 an input that results in multiple precision computation.  Argument
Index: glibc-2.26/math/Makefile
===================================================================
--- glibc-2.26.orig/math/Makefile
+++ glibc-2.26/math/Makefile
@@ -110,9 +110,9 @@ type-ldouble-yes := ldouble
 
 # double support
 type-double-suffix :=
-type-double-routines := branred doasin dosincos halfulp mpa mpatan2	\
-		       mpatan mpexp mplog mpsqrt mptan sincos32 slowexp	\
-		       slowpow sincostab k_rem_pio2
+type-double-routines := branred doasin dosincos mpa mpatan2	\
+		       mpatan mpexp mplog mpsqrt mptan sincos32	\
+		       sincostab k_rem_pio2
 
 # float support
 type-float-suffix := f
Index: glibc-2.26/sysdeps/aarch64/libm-test-ulps
===================================================================
--- glibc-2.26.orig/sysdeps/aarch64/libm-test-ulps
+++ glibc-2.26/sysdeps/aarch64/libm-test-ulps
@@ -1012,7 +1012,9 @@ ildouble: 2
 ldouble: 2
 
 Function: "cos":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
@@ -1932,7 +1934,9 @@ ildouble: 1
 ldouble: 1
 
 Function: "pow":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 2
 ldouble: 2
@@ -1992,7 +1996,9 @@ ildouble: 2
 ldouble: 2
 
 Function: "sin":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
@@ -2022,7 +2028,9 @@ ildouble: 3
 ldouble: 3
 
 Function: "sincos":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1
Index: glibc-2.26/sysdeps/generic/math_private.h
===================================================================
--- glibc-2.26.orig/sysdeps/generic/math_private.h
+++ glibc-2.26/sysdeps/generic/math_private.h
@@ -255,20 +255,17 @@ extern float __kernel_standard_f (float,
 extern long double __kernel_standard_l (long double,long double,int);
 
 /* Prototypes for functions of the IBM Accurate Mathematical Library.  */
-extern double __exp1 (double __x, double __xx, double __error);
+extern double __exp1 (double __x, double __xx);
 extern double __sin (double __x);
 extern double __cos (double __x);
 extern int __branred (double __x, double *__a, double *__aa);
 extern void __doasin (double __x, double __dx, double __v[]);
 extern void __dubsin (double __x, double __dx, double __v[]);
 extern void __dubcos (double __x, double __dx, double __v[]);
-extern double __halfulp (double __x, double __y);
 extern double __sin32 (double __x, double __res, double __res1);
 extern double __cos32 (double __x, double __res, double __res1);
 extern double __mpsin (double __x, double __dx, bool __range_reduce);
 extern double __mpcos (double __x, double __dx, bool __range_reduce);
-extern double __slowexp (double __x);
-extern double __slowpow (double __x, double __y, double __z);
 extern void __docos (double __x, double __dx, double __v[]);
 
 #ifndef math_opt_barrier
Index: glibc-2.26/sysdeps/i386/fpu/halfulp.c
===================================================================
--- glibc-2.26.orig/sysdeps/i386/fpu/halfulp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/i386/fpu/slowexp.c
===================================================================
--- glibc-2.26.orig/sysdeps/i386/fpu/slowexp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/i386/fpu/slowpow.c
===================================================================
--- glibc-2.26.orig/sysdeps/i386/fpu/slowpow.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/ia64/fpu/halfulp.c
===================================================================
--- glibc-2.26.orig/sysdeps/ia64/fpu/halfulp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/ia64/fpu/slowexp.c
===================================================================
--- glibc-2.26.orig/sysdeps/ia64/fpu/slowexp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/ia64/fpu/slowpow.c
===================================================================
--- glibc-2.26.orig/sysdeps/ia64/fpu/slowpow.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/ieee754/dbl-64/e_exp.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/e_exp.c
+++ glibc-2.26/sysdeps/ieee754/dbl-64/e_exp.c
@@ -23,10 +23,10 @@
 /*           exp1                                                          */
 /*                                                                         */
 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
-/*              mpa.c mpexp.x slowexp.c                                    */
+/*              mpa.c mpexp.x                                              */
 /*                                                                         */
 /* An ultimate exp routine. Given an IEEE double machine number x          */
-/* it computes the correctly rounded (to nearest) value of e^x             */
+/* it computes an almost correctly rounded (to nearest) value of e^x       */
 /* Assumption: Machine arithmetic operations are performed in              */
 /* round to nearest mode of IEEE 754 standard.                             */
 /*                                                                         */
@@ -46,10 +46,6 @@
 # define SECTION
 #endif
 
-double __slowexp (double);
-
-/* An ultimate exp routine. Given an IEEE double machine number x it computes
-   the correctly rounded (to nearest) value of e^x.  */
 double
 SECTION
 __ieee754_exp (double x)
@@ -93,17 +89,10 @@ __ieee754_exp (double x)
 
 	rem = (bet + bet * eps) + al * eps;
 	res = al + rem;
-	cor = (al - res) + rem;
-	if (res == (res + cor * err_0))
-	  {
-	    retval = res * binexp.x;
-	    goto ret;
-	  }
-	else
-	  {
-	    retval = __slowexp (x);
-	    goto ret;
-	  }			/*if error is over bound */
+	/* Maximum relative error is 7.8e-22 (70.1 bits).
+	   Maximum ULP error is 0.500007.  */
+	retval = res * binexp.x;
+	goto ret;
       }
 
     if (n <= smallint)
@@ -166,38 +155,22 @@ __ieee754_exp (double x)
 	if (ex >= -1022)
 	  {
 	    binexp.i[HIGH_HALF] = (1023 + ex) << 20;
-	    if (res == (res + cor * err_0))
-	      {
-		retval = res * binexp.x;
-		goto ret;
-	      }
-	    else
-	      {
-		retval = __slowexp (x);
-		goto check_uflow_ret;
-	      }			/*if error is over bound */
+	    /* Does not underflow: res >= 1.0, binexp >= 0x1p-1022
+	       Maximum relative error is 7.8e-22 (70.1 bits).
+	       Maximum ULP error is 0.500007.  */
+	    retval = res * binexp.x;
+	    goto ret;
 	  }
 	ex = -(1022 + ex);
 	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
 	res *= binexp.x;
 	cor *= binexp.x;
-	eps = 1.0000000001 + err_0 * binexp.x;
 	t = 1.0 + res;
 	y = ((1.0 - t) + res) + cor;
 	res = t + y;
-	cor = (t - res) + y;
-	if (res == (res + eps * cor))
-	  {
-	    binexp.i[HIGH_HALF] = 0x00100000;
-	    retval = (res - 1.0) * binexp.x;
-	    goto check_uflow_ret;
-	  }
-	else
-	  {
-	    retval = __slowexp (x);
-	    goto check_uflow_ret;
-	  }			/*   if error is over bound    */
-      check_uflow_ret:
+	/* Maximum ULP error is 0.5000035.  */
+	binexp.i[HIGH_HALF] = 0x00100000;
+	retval = (res - 1.0) * binexp.x;
 	if (retval < DBL_MIN)
 	  {
 	    double force_underflow = tiny * tiny;
@@ -210,10 +183,9 @@ __ieee754_exp (double x)
     else
       {
 	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
-	if (res == (res + cor * err_0))
-	  retval = res * binexp.x * t256.x;
-	else
-	  retval = __slowexp (x);
+	/* Maximum relative error is 7.8e-22 (70.1 bits).
+	   Maximum ULP error is 0.500007.  */
+	retval = res * binexp.x * t256.x;
 	if (isinf (retval))
 	  goto ret_huge;
 	else
@@ -233,13 +205,10 @@ ret:
 strong_alias (__ieee754_exp, __exp_finite)
 #endif
 
-/* Compute e^(x+xx).  The routine also receives bound of error of previous
-   calculation.  If after computing exp the error exceeds the allowed bounds,
-   the routine returns a non-positive number.  Otherwise it returns the
-   computed result, which is always positive.  */
+/* Compute e^(x+xx).  */
 double
 SECTION
-__exp1 (double x, double xx, double error)
+__exp1 (double x, double xx)
 {
   double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
   mynumber junk1, junk2, binexp = {{0, 0}};
@@ -249,6 +218,7 @@ __exp1 (double x, double xx, double erro
   m = junk1.i[HIGH_HALF];
   n = m & hugeint;		/* no sign */
 
+  /* fabs (x) > 5.551112e-17 and fabs (x) < 7.080010e+02.  */
   if (n > smallint && n < bigint)
     {
       y = x * log2e.x + three51.x;
@@ -276,11 +246,9 @@ __exp1 (double x, double xx, double erro
 
       rem = (bet + bet * eps) + al * eps;
       res = al + rem;
-      cor = (al - res) + rem;
-      if (res == (res + cor * (1.0 + error + err_1)))
-	return res * binexp.x;
-      else
-	return -10.0;
+      /* Maximum relative error before rounding is 8.8e-22 (69.9 bits).
+	 Maximum ULP error is 0.500008.  */
+      return res * binexp.x;
     }
 
   if (n <= smallint)
@@ -318,6 +286,7 @@ __exp1 (double x, double xx, double erro
   cor = (al - res) + rem;
   if (m >> 31)
     {
+      /* x < 0.  */
       ex = junk1.i[LOW_HALF];
       if (res < 1.0)
 	{
@@ -328,34 +297,25 @@ __exp1 (double x, double xx, double erro
       if (ex >= -1022)
 	{
 	  binexp.i[HIGH_HALF] = (1023 + ex) << 20;
-	  if (res == (res + cor * (1.0 + error + err_1)))
-	    return res * binexp.x;
-	  else
-	    return -10.0;
+	  /* Maximum ULP error is 0.500008.  */
+	  return res * binexp.x;
 	}
+      /* Denormal case - ex < -1022.  */
       ex = -(1022 + ex);
       binexp.i[HIGH_HALF] = (1023 - ex) << 20;
       res *= binexp.x;
       cor *= binexp.x;
-      eps = 1.00000000001 + (error + err_1) * binexp.x;
       t = 1.0 + res;
       y = ((1.0 - t) + res) + cor;
       res = t + y;
-      cor = (t - res) + y;
-      if (res == (res + eps * cor))
-	{
-	  binexp.i[HIGH_HALF] = 0x00100000;
-	  return (res - 1.0) * binexp.x;
-	}
-      else
-	return -10.0;
+      binexp.i[HIGH_HALF] = 0x00100000;
+      /* Maximum ULP error is 0.500004.  */
+      return (res - 1.0) * binexp.x;
     }
   else
     {
       binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
-      if (res == (res + cor * (1.0 + error + err_1)))
-	return res * binexp.x * t256.x;
-      else
-	return -10.0;
+      /* Maximum ULP error is 0.500008.  */
+      return res * binexp.x * t256.x;
     }
 }
Index: glibc-2.26/sysdeps/ieee754/dbl-64/e_log.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/e_log.c
+++ glibc-2.26/sysdeps/ieee754/dbl-64/e_log.c
@@ -23,11 +23,10 @@
 /*      FUNCTION:ulog                                                */
 /*                                                                   */
 /*      FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h           */
-/*                    mpexp.c mplog.c mpa.c                          */
 /*                    ulog.tbl                                       */
 /*                                                                   */
 /* An ultimate log routine. Given an IEEE double machine number x    */
-/* it computes the correctly rounded (to nearest) value of log(x).   */
+/* it computes the rounded (to nearest) value of log(x).	     */
 /* Assumption: Machine arithmetic operations are performed in        */
 /* round to nearest mode of IEEE 754 standard.                       */
 /*                                                                   */
@@ -40,34 +39,26 @@
 #include "MathLib.h"
 #include <math.h>
 #include <math_private.h>
-#include <stap-probe.h>
 
 #ifndef SECTION
 # define SECTION
 #endif
 
-void __mplog (mp_no *, mp_no *, int);
-
 /*********************************************************************/
-/* An ultimate log routine. Given an IEEE double machine number x     */
-/* it computes the correctly rounded (to nearest) value of log(x).   */
+/* An ultimate log routine. Given an IEEE double machine number x    */
+/* it computes the rounded (to nearest) value of log(x).	     */
 /*********************************************************************/
 double
 SECTION
 __ieee754_log (double x)
 {
-#define M 4
-  static const int pr[M] = { 8, 10, 18, 32 };
-  int i, j, n, ux, dx, p;
+  int i, j, n, ux, dx;
   double dbl_n, u, p0, q, r0, w, nln2a, luai, lubi, lvaj, lvbj,
-	 sij, ssij, ttij, A, B, B0, y, y1, y2, polI, polII, sa, sb,
-	 t1, t2, t7, t8, t, ra, rb, ww,
-	 a0, aa0, s1, s2, ss2, s3, ss3, a1, aa1, a, aa, b, bb, c;
+	 sij, ssij, ttij, A, B, B0, polI, polII, t8, a, aa, b, bb, c;
 #ifndef DLA_FMS
-  double t3, t4, t5, t6;
+  double t1, t2, t3, t4, t5;
 #endif
   number num;
-  mp_no mpx, mpy, mpy1, mpy2, mperr;
 
 #include "ulog.tbl"
 #include "ulog.h"
@@ -101,7 +92,7 @@ __ieee754_log (double x)
   if (w == 0.0)
     return 0.0;
 
-  /*--- Stage I, the case abs(x-1) < 0.03 */
+  /*--- The case abs(x-1) < 0.03 */
 
   t8 = MHALF * w;
   EMULV (t8, w, a, aa, t1, t2, t3, t4, t5);
@@ -118,50 +109,12 @@ __ieee754_log (double x)
   polII *= w * w * w;
   c = (aa + bb) + polII;
 
-  /* End stage I, case abs(x-1) < 0.03 */
-  if ((y = b + (c + b * E2)) == b + (c - b * E2))
-    return y;
-
-  /*--- Stage II, the case abs(x-1) < 0.03 */
-
-  a = d19.d + w * d20.d;
-  a = d18.d + w * a;
-  a = d17.d + w * a;
-  a = d16.d + w * a;
-  a = d15.d + w * a;
-  a = d14.d + w * a;
-  a = d13.d + w * a;
-  a = d12.d + w * a;
-  a = d11.d + w * a;
-
-  EMULV (w, a, s2, ss2, t1, t2, t3, t4, t5);
-  ADD2 (d10.d, dd10.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d9.d, dd9.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d8.d, dd8.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d7.d, dd7.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d6.d, dd6.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d5.d, dd5.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d4.d, dd4.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d3.d, dd3.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (d2.d, dd2.d, s2, ss2, s3, ss3, t1, t2);
-  MUL2 (w, 0, s3, ss3, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  MUL2 (w, 0, s2, ss2, s3, ss3, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (w, 0, s3, ss3, b, bb, t1, t2);
-
-  /* End stage II, case abs(x-1) < 0.03 */
-  if ((y = b + (bb + b * E4)) == b + (bb - b * E4))
-    return y;
-  goto stage_n;
+  /* Here b contains the high part of the result, and c the low part.
+     Maximum error is b * 2.334e-19, so accuracy is >61 bits.
+     Therefore max ULP error of b + c is ~0.502.  */
+  return b + c;
 
-  /*--- Stage I, the case abs(x-1) > 0.03 */
+  /*--- The case abs(x-1) > 0.03 */
 case_03:
 
   /* Find n,u such that x = u*2**n,   1/sqrt(2) < u < sqrt(2)  */
@@ -203,58 +156,10 @@ case_03:
   B0 = (((lubi + lvbj) + ssij) + ttij) + dbl_n * LN2B;
   B = polI + B0;
 
-  /* End stage I, case abs(x-1) >= 0.03 */
-  if ((y = A + (B + E1)) == A + (B - E1))
-    return y;
-
-
-  /*--- Stage II, the case abs(x-1) > 0.03 */
-
-  /* Improve the accuracy of r0 */
-  EMULV (p0, r0, sa, sb, t1, t2, t3, t4, t5);
-  t = r0 * ((1 - sa) - sb);
-  EADD (r0, t, ra, rb);
-
-  /* Compute w */
-  MUL2 (q, 0, ra, rb, w, ww, t1, t2, t3, t4, t5, t6, t7, t8);
-
-  EADD (A, B0, a0, aa0);
-
-  /* Evaluate polynomial III */
-  s1 = (c3.d + (c4.d + c5.d * w) * w) * w;
-  EADD (c2.d, s1, s2, ss2);
-  MUL2 (s2, ss2, w, ww, s3, ss3, t1, t2, t3, t4, t5, t6, t7, t8);
-  MUL2 (s3, ss3, w, ww, s2, ss2, t1, t2, t3, t4, t5, t6, t7, t8);
-  ADD2 (s2, ss2, w, ww, s3, ss3, t1, t2);
-  ADD2 (s3, ss3, a0, aa0, a1, aa1, t1, t2);
-
-  /* End stage II, case abs(x-1) >= 0.03 */
-  if ((y = a1 + (aa1 + E3)) == a1 + (aa1 - E3))
-    return y;
-
-
-  /* Final stages. Use multi-precision arithmetic. */
-stage_n:
-
-  for (i = 0; i < M; i++)
-    {
-      p = pr[i];
-      __dbl_mp (x, &mpx, p);
-      __dbl_mp (y, &mpy, p);
-      __mplog (&mpx, &mpy, p);
-      __dbl_mp (e[i].d, &mperr, p);
-      __add (&mpy, &mperr, &mpy1, p);
-      __sub (&mpy, &mperr, &mpy2, p);
-      __mp_dbl (&mpy1, &y1, p);
-      __mp_dbl (&mpy2, &y2, p);
-      if (y1 == y2)
-	{
-	  LIBC_PROBE (slowlog, 3, &p, &x, &y1);
-	  return y1;
-	}
-    }
-  LIBC_PROBE (slowlog_inexact, 3, &p, &x, &y1);
-  return y1;
+  /* Here A contains the high part of the result, and B the low part.
+     Maximum abs error is 6.095e-21 and min log (x) is 0.0295 since x > 1.03.
+     Therefore max ULP error of A + B is ~0.502.  */
+  return A + B;
 }
 
 #ifndef __ieee754_log
Index: glibc-2.26/sysdeps/ieee754/dbl-64/e_pow.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/e_pow.c
+++ glibc-2.26/sysdeps/ieee754/dbl-64/e_pow.c
@@ -20,13 +20,9 @@
 /*  MODULE_NAME: upow.c                                                    */
 /*                                                                         */
 /*  FUNCTIONS: upow                                                        */
-/*             power1                                                      */
-/*             my_log2                                                     */
 /*             log1                                                        */
 /*             checkint                                                    */
 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
-/*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
-/*                          uexp.c  upow.c				   */
 /*               root.tbl uexp.tbl upow.tbl                                */
 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
 /* it computes the correctly rounded (to nearest) value of x^y.            */
@@ -50,11 +46,8 @@
 
 static const double huge = 1.0e300, tiny = 1.0e-300;
 
-double __exp1 (double x, double xx, double error);
-static double log1 (double x, double *delta, double *error);
-static double my_log2 (double x, double *delta, double *error);
-double __slowpow (double x, double y, double z);
-static double power1 (double x, double y);
+double __exp1 (double x, double xx);
+static double log1 (double x, double *delta);
 static int checkint (double x);
 
 /* An ultimate power routine. Given two IEEE double machine numbers y, x it
@@ -63,7 +56,7 @@ double
 SECTION
 __ieee754_pow (double x, double y)
 {
-  double z, a, aa, error, t, a1, a2, y1, y2;
+  double z, a, aa, t, a1, a2, y1, y2;
   mynumber u, v;
   int k;
   int4 qx, qy;
@@ -100,7 +93,7 @@ __ieee754_pow (double x, double y)
 	   not matter if |y| <= 2**-64.  */
 	if (fabs (y) < 0x1p-64)
 	  y = y < 0 ? -0x1p-64 : 0x1p-64;
-	z = log1 (x, &aa, &error);	/* x^y  =e^(y log (X)) */
+	z = log1 (x, &aa);	/* x^y  =e^(y log (X)) */
 	t = y * CN;
 	y1 = t - (t - y);
 	y2 = y - y1;
@@ -111,9 +104,16 @@ __ieee754_pow (double x, double y)
 	aa = y2 * a1 + y * a2;
 	a1 = a + aa;
 	a2 = (a - a1) + aa;
-	error = error * fabs (y);
-	t = __exp1 (a1, a2, 1.9e16 * error);	/* return -10 or 0 if wasn't computed exactly */
-	retval = (t > 0) ? t : power1 (x, y);
+
+	/* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits).
+	   Maximum relative error REexp of __exp1 is 8.8e-22 (69.9 bits).
+	   We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp).
+	   Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX),
+	   this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp).
+	   So the relative error is 710 * 1.0e-21 + 8.8e-22 = 7.1e-19
+	   (60.2 bits).  The worst-case ULP error is 0.5064.  */
+
+	retval = __exp1 (a1, a2);
       }
 
       if (isinf (retval))
@@ -218,33 +218,11 @@ __ieee754_pow (double x, double y)
 strong_alias (__ieee754_pow, __pow_finite)
 #endif
 
-/* Compute x^y using more accurate but more slow log routine.  */
-static double
-SECTION
-power1 (double x, double y)
-{
-  double z, a, aa, error, t, a1, a2, y1, y2;
-  z = my_log2 (x, &aa, &error);
-  t = y * CN;
-  y1 = t - (t - y);
-  y2 = y - y1;
-  t = z * CN;
-  a1 = t - (t - z);
-  a2 = z - a1;
-  a = y * z;
-  aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y;
-  a1 = a + aa;
-  a2 = (a - a1) + aa;
-  error = error * fabs (y);
-  t = __exp1 (a1, a2, 1.9e16 * error);
-  return (t >= 0) ? t : __slowpow (x, y, z);
-}
-
 /* Compute log(x) (x is left argument). The result is the returned double + the
-   parameter DELTA.  The result is bounded by ERROR.  */
+   parameter DELTA.  */
 static double
 SECTION
-log1 (double x, double *delta, double *error)
+log1 (double x, double *delta)
 {
   unsigned int i, j;
   int m;
@@ -260,9 +238,7 @@ log1 (double x, double *delta, double *e
 
   u.x = x;
   m = u.i[HIGH_HALF];
-  *error = 0;
-  *delta = 0;
-  if (m < 0x00100000)		/*  1<x<2^-1007 */
+  if (m < 0x00100000)		/* Handle denormal x.  */
     {
       x = x * t52.x;
       add = -52.0;
@@ -284,7 +260,7 @@ log1 (double x, double *delta, double *e
   v.x = u.x + bigu.x;
   uu = v.x - bigu.x;
   i = (v.i[LOW_HALF] & 0x000003ff) << 2;
-  if (two52.i[LOW_HALF] == 1023)	/* nx = 0              */
+  if (two52.i[LOW_HALF] == 1023)	/* Exponent of x is 0.  */
     {
       if (i > 1192 && i < 1208)	/* |x-1| < 1.5*2**-10  */
 	{
@@ -296,8 +272,8 @@ log1 (double x, double *delta, double *e
 							   * (r7 + t * r8)))))
 		- 0.5 * t2 * (t + t1));
 	  res = e1 + e2;
-	  *error = 1.0e-21 * fabs (t);
 	  *delta = (e1 - res) + e2;
+	  /* Max relative error is 1.464844e-24, so accurate to 79.1 bits.  */
 	  return res;
 	}			/* |x-1| < 1.5*2**-10  */
       else
@@ -316,12 +292,12 @@ log1 (double x, double *delta, double *e
 	  t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e
 		* (p2 + e * (p3 + e * p4)));
 	  res = t1 + t2;
-	  *error = 1.0e-24;
 	  *delta = (t1 - res) + t2;
+	  /* Max relative error is 1.0e-24, so accurate to 79.7 bits.  */
 	  return res;
 	}
-    }				/* nx = 0 */
-  else				/* nx != 0   */
+    }
+  else				/* Exponent of x != 0.  */
     {
       eps = u.x - uu;
       nx = (two52.x - two52e.x) + add;
@@ -334,113 +310,13 @@ log1 (double x, double *delta, double *e
       t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e
 	    * (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6)))));
       res = t1 + t2;
-      *error = 1.0e-21;
-      *delta = (t1 - res) + t2;
-      return res;
-    }				/* nx != 0   */
-}
-
-/* Slower but more accurate routine of log.  The returned result is double +
-   DELTA.  The result is bounded by ERROR.  */
-static double
-SECTION
-my_log2 (double x, double *delta, double *error)
-{
-  unsigned int i, j;
-  int m;
-  double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
-  double ou1, ou2, lu1, lu2, ov, lv1, lv2, a, a1, a2;
-  double y, yy, z, zz, j1, j2, j7, j8;
-#ifndef DLA_FMS
-  double j3, j4, j5, j6;
-#endif
-  mynumber u, v;
-#ifdef BIG_ENDI
-  mynumber /**/ two52 = {{0x43300000, 0x00000000}};	/* 2**52  */
-#else
-# ifdef LITTLE_ENDI
-  mynumber /**/ two52 = {{0x00000000, 0x43300000}};	/* 2**52  */
-# endif
-#endif
-
-  u.x = x;
-  m = u.i[HIGH_HALF];
-  *error = 0;
-  *delta = 0;
-  add = 0;
-  if (m < 0x00100000)
-    {				/* x < 2^-1022 */
-      x = x * t52.x;
-      add = -52.0;
-      u.x = x;
-      m = u.i[HIGH_HALF];
-    }
-
-  if ((m & 0x000fffff) < 0x0006a09e)
-    {
-      u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
-      two52.i[LOW_HALF] = (m >> 20);
-    }
-  else
-    {
-      u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
-      two52.i[LOW_HALF] = (m >> 20) + 1;
-    }
-
-  v.x = u.x + bigu.x;
-  uu = v.x - bigu.x;
-  i = (v.i[LOW_HALF] & 0x000003ff) << 2;
-  /*------------------------------------- |x-1| < 2**-11-------------------------------  */
-  if ((two52.i[LOW_HALF] == 1023) && (i == 1200))
-    {
-      t = x - 1.0;
-      EMULV (t, s3, y, yy, j1, j2, j3, j4, j5);
-      ADD2 (-0.5, 0, y, yy, z, zz, j1, j2);
-      MUL2 (t, 0, z, zz, y, yy, j1, j2, j3, j4, j5, j6, j7, j8);
-      MUL2 (t, 0, y, yy, z, zz, j1, j2, j3, j4, j5, j6, j7, j8);
-
-      e1 = t + z;
-      e2 = ((((t - e1) + z) + zz) + t * t * t
-	    * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8))))));
-      res = e1 + e2;
-      *error = 1.0e-25 * fabs (t);
-      *delta = (e1 - res) + e2;
-      return res;
-    }
-  /*----------------------------- |x-1| > 2**-11  --------------------------  */
-  else
-    {				/*Computing log(x) according to log table                        */
-      nx = (two52.x - two52e.x) + add;
-      ou1 = ui.x[i];
-      ou2 = ui.x[i + 1];
-      lu1 = ui.x[i + 2];
-      lu2 = ui.x[i + 3];
-      v.x = u.x * (ou1 + ou2) + bigv.x;
-      vv = v.x - bigv.x;
-      j = v.i[LOW_HALF] & 0x0007ffff;
-      j = j + j + j;
-      eps = u.x - uu * vv;
-      ov = vj.x[j];
-      lv1 = vj.x[j + 1];
-      lv2 = vj.x[j + 2];
-      a = (ou1 + ou2) * (1.0 + ov);
-      a1 = (a + 1.0e10) - 1.0e10;
-      a2 = a * (1.0 - a1 * uu * vv);
-      e1 = eps * a1;
-      e2 = eps * a2;
-      e = e1 + e2;
-      e2 = (e1 - e) + e2;
-      t = nx * ln2a.x + lu1 + lv1;
-      t1 = t + e;
-      t2 = ((((t - t1) + e) + (lu2 + lv2 + nx * ln2b.x + e2)) + e * e
-	    * (p2 + e * (p3 + e * p4)));
-      res = t1 + t2;
-      *error = 1.0e-27;
       *delta = (t1 - res) + t2;
+      /* Max relative error is 1.0e-21, so accurate to 69.7 bits.  */
       return res;
     }
 }
 
+
 /* This function receives a double x and checks if it is an integer.  If not,
    it returns 0, else it returns 1 if even or -1 if odd.  */
 static int
Index: glibc-2.26/sysdeps/ieee754/dbl-64/halfulp.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/halfulp.c
+++ /dev/null
@@ -1,152 +0,0 @@
-/*
- * 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:halfulp.c                                                */
-/*                                                                      */
-/*  FUNCTIONS:halfulp                                                   */
-/*  FILES NEEDED: mydefs.h dla.h endian.h                               */
-/*                uroot.c                                               */
-/*                                                                      */
-/*Routine halfulp(double x, double y) computes x^y where result does    */
-/*not need rounding. If the result is closer to 0 than can be           */
-/*represented it returns 0.                                             */
-/*     In the following cases the function does not compute anything    */
-/*and returns a negative number:                                        */
-/*1. if the result needs rounding,                                      */
-/*2. if y is outside the interval [0,  2^20-1],                         */
-/*3. if x can be represented by  x=2**n for some integer n.             */
-/************************************************************************/
-
-#include "endian.h"
-#include "mydefs.h"
-#include <dla.h>
-#include <math_private.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-static const int4 tab54[32] = {
-  262143, 11585,  1782, 511, 210, 107, 63, 42,
-  30,     22,     17,   14,  12,  10,  9, 7,
-  7,      6,      5,    5,   5,   4,   4, 4,
-  3,      3,      3,    3,   3,   3,   3, 3
-};
-
-
-double
-SECTION
-__halfulp (double x, double y)
-{
-  mynumber v;
-  double z, u, uu;
-#ifndef DLA_FMS
-  double j1, j2, j3, j4, j5;
-#endif
-  int4 k, l, m, n;
-  if (y <= 0)                 /*if power is negative or zero */
-    {
-      v.x = y;
-      if (v.i[LOW_HALF] != 0)
-	return -10.0;
-      v.x = x;
-      if (v.i[LOW_HALF] != 0)
-	return -10.0;
-      if ((v.i[HIGH_HALF] & 0x000fffff) != 0)
-	return -10;                                     /* if x =2 ^ n */
-      k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */
-      z = (double) k;
-      return (z * y == -1075.0) ? 0 : -10.0;
-    }
-  /* if y > 0  */
-  v.x = y;
-  if (v.i[LOW_HALF] != 0)
-    return -10.0;
-
-  v.x = x;
-  /*  case where x = 2**n for some integer n */
-  if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0)
-    {
-      k = (v.i[HIGH_HALF] >> 20) - 1023;
-      return (((double) k) * y == -1075.0) ? 0 : -10.0;
-    }
-
-  v.x = y;
-  k = v.i[HIGH_HALF];
-  m = k << 12;
-  l = 0;
-  while (m)
-    {
-      m = m << 1; l++;
-    }
-  n = (k & 0x000fffff) | 0x00100000;
-  n = n >> (20 - l);                       /*   n is the odd integer of y    */
-  k = ((k >> 20) - 1023) - l;               /*   y = n*2**k                   */
-  if (k > 5)
-    return -10.0;
-  if (k > 0)
-    for (; k > 0; k--)
-      n *= 2;
-  if (n > 34)
-    return -10.0;
-  k = -k;
-  if (k > 5)
-    return -10.0;
-
-  /*   now treat x        */
-  while (k > 0)
-    {
-      z = __ieee754_sqrt (x);
-      EMULV (z, z, u, uu, j1, j2, j3, j4, j5);
-      if (((u - x) + uu) != 0)
-	break;
-      x = z;
-      k--;
-    }
-  if (k)
-    return -10.0;
-
-  /* it is impossible that n == 2,  so the mantissa of x must be short  */
-
-  v.x = x;
-  if (v.i[LOW_HALF])
-    return -10.0;
-  k = v.i[HIGH_HALF];
-  m = k << 12;
-  l = 0;
-  while (m)
-    {
-      m = m << 1; l++;
-    }
-  m = (k & 0x000fffff) | 0x00100000;
-  m = m >> (20 - l);                       /*   m is the odd integer of x    */
-
-  /*   now check whether the length of m**n is at most 54 bits */
-
-  if (m > tab54[n - 3])
-    return -10.0;
-
-  /* yes, it is - now compute x**n by simple multiplications  */
-
-  u = x;
-  for (k = 1; k < n; k++)
-    u = u * x;
-  return u;
-}
Index: glibc-2.26/sysdeps/ieee754/dbl-64/s_sin.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/s_sin.c
+++ glibc-2.26/sysdeps/ieee754/dbl-64/s_sin.c
@@ -22,22 +22,11 @@
 /*                                                                          */
 /* FUNCTIONS: usin                                                          */
 /*            ucos                                                          */
-/*            slow                                                          */
-/*            slow1                                                         */
-/*            slow2                                                         */
-/*            sloww                                                         */
-/*            sloww1                                                        */
-/*            sloww2                                                        */
-/*            bsloww                                                        */
-/*            bsloww1                                                       */
-/*            bsloww2                                                       */
-/*            cslow2                                                        */
 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h  usncs.h                     */
-/*               branred.c sincos32.c dosincos.c mpa.c                      */
-/*               sincos.tbl                                                 */
+/*		 branred.c sincos.tbl					    */
 /*                                                                          */
-/* An ultimate sin and  routine. Given an IEEE double machine number x       */
-/* it computes the correctly rounded (to nearest) value of sin(x) or cos(x) */
+/* An ultimate sin and cos routine. Given an IEEE double machine number x   */
+/* it computes sin(x) or cos(x) with ~0.55 ULP.				    */
 /* Assumption: Machine arithmetic operations are performed in               */
 /* round to nearest mode of IEEE 754 standard.                              */
 /*                                                                          */
@@ -65,35 +54,11 @@
    a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
 
    The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
-   on.  The result is returned to LHS and correction in COR.  */
-#define TAYLOR_SIN(xx, a, da, cor) \
+   on.  The result is returned to LHS.  */
+#define TAYLOR_SIN(xx, a, da) \
 ({									      \
   double t = ((POLYNOMIAL (xx)  * (a) - 0.5 * (da))  * (xx) + (da));	      \
   double res = (a) + t;							      \
-  (cor) = ((a) - res) + t;						      \
-  res;									      \
-})
-
-/* This is again a variation of the Taylor series expansion with the term
-   x^3/3! expanded into the following for better accuracy:
-
-   bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3
-
-   The correction term is dx and bb + aa = -1/3!
-   */
-#define TAYLOR_SLOW(x0, dx, cor) \
-({									      \
-  static const double th2_36 = 206158430208.0;	/*    1.5*2**37   */	      \
-  double xx = (x0) * (x0);						      \
-  double x1 = ((x0) + th2_36) - th2_36;					      \
-  double y = aa * x1 * x1 * x1;						      \
-  double r = (x0) + y;							      \
-  double x2 = ((x0) - x1) + (dx);					      \
-  double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2)	      \
-	      * (x0)  + aa * x2 * x2 * x2 + (dx));			      \
-  t = (((x0) - r) + y) + t;						      \
-  double res = r + t;							      \
-  (cor) = (r - res) + t;						      \
   res;									      \
 })
 
@@ -123,31 +88,15 @@ static const double
   cs4 = -4.16666666666664434524222570944589E-02,
   cs6 = 1.38888874007937613028114285595617E-03;
 
-static const double t22 = 0x1.8p22;
-
-void __dubsin (double x, double dx, double w[]);
-void __docos (double x, double dx, double w[]);
-double __mpsin (double x, double dx, bool reduce_range);
-double __mpcos (double x, double dx, bool reduce_range);
-static double slow (double x);
-static double slow1 (double x);
-static double slow2 (double x);
-static double sloww (double x, double dx, double orig, bool shift_quadrant);
-static double sloww1 (double x, double dx, double orig, bool shift_quadrant);
-static double sloww2 (double x, double dx, double orig, int n);
-static double bsloww (double x, double dx, double orig, int n);
-static double bsloww1 (double x, double dx, double orig, int n);
-static double bsloww2 (double x, double dx, double orig, int n);
 int __branred (double x, double *a, double *aa);
-static double cslow2 (double x);
 
 /* Given a number partitioned into X and DX, this function computes the cosine
    of the number by combining the sin and cos of X (as computed by a variation
    of the Taylor series) with the values looked up from the sin/cos table to
-   get the result in RES and a correction value in COR.  */
+   get the result.  */
 static inline double
 __always_inline
-do_cos (double x, double dx, double *corp)
+do_cos (double x, double dx)
 {
   mynumber u;
 
@@ -157,60 +106,28 @@ do_cos (double x, double dx, double *cor
   u.x = big + fabs (x);
   x = fabs (x) - (u.x - big) + dx;
 
-  double xx, s, sn, ssn, c, cs, ccs, res, cor;
+  double xx, s, sn, ssn, c, cs, ccs, cor;
   xx = x * x;
   s = x + x * xx * (sn3 + xx * sn5);
   c = xx * (cs2 + xx * (cs4 + xx * cs6));
   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
   cor = (ccs - s * ssn - cs * c) - sn * s;
-  res = cs + cor;
-  cor = (cs - res) + cor;
-  *corp = cor;
-  return res;
-}
-
-/* A more precise variant of DO_COS.  EPS is the adjustment to the correction
-   COR.  */
-static inline double
-__always_inline
-do_cos_slow (double x, double dx, double eps, double *corp)
-{
-  mynumber u;
-
-  if (x <= 0)
-    dx = -dx;
-
-  u.x = big + fabs (x);
-  x = fabs (x) - (u.x - big);
-
-  double xx, y, x1, x2, e1, e2, res, cor;
-  double s, sn, ssn, c, cs, ccs;
-  xx = x * x;
-  s = x * xx * (sn3 + xx * sn5);
-  c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
-  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
-  x1 = (x + t22) - t22;
-  x2 = (x - x1) + dx;
-  e1 = (sn + t22) - t22;
-  e2 = (sn - e1) + ssn;
-  cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s;
-  y = cs - e1 * x1;
-  cor = cor + ((cs - y) - e1 * x1);
-  res = y + cor;
-  cor = (y - res) + cor;
-  cor = 1.0005 * cor + __copysign (eps, cor);
-  *corp = cor;
-  return res;
+  return cs + cor;
 }
 
 /* Given a number partitioned into X and DX, this function computes the sine of
    the number by combining the sin and cos of X (as computed by a variation of
    the Taylor series) with the values looked up from the sin/cos table to get
-   the result in RES and a correction value in COR.  */
+   the result.  */
 static inline double
 __always_inline
-do_sin (double x, double dx, double *corp)
+do_sin (double x, double dx)
 {
+  double xold = x;
+  /* Max ULP is 0.501 if |x| < 0.126, otherwise ULP is 0.518.  */
+  if (fabs (x) < 0.126)
+    return TAYLOR_SIN (x * x, x, dx);
+
   mynumber u;
 
   if (x <= 0)
@@ -218,85 +135,22 @@ do_sin (double x, double dx, double *cor
   u.x = big + fabs (x);
   x = fabs (x) - (u.x - big);
 
-  double xx, s, sn, ssn, c, cs, ccs, cor, res;
+  double xx, s, sn, ssn, c, cs, ccs, cor;
   xx = x * x;
   s = x + (dx + x * xx * (sn3 + xx * sn5));
   c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
   SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
   cor = (ssn + s * ccs - sn * c) + cs * s;
-  res = sn + cor;
-  cor = (sn - res) + cor;
-  *corp = cor;
-  return res;
-}
-
-/* A more precise variant of DO_SIN.  EPS is the adjustment to the correction
-   COR.  */
-static inline double
-__always_inline
-do_sin_slow (double x, double dx, double eps, double *corp)
-{
-  mynumber u;
-
-  if (x <= 0)
-    dx = -dx;
-  u.x = big + fabs (x);
-  x = fabs (x) - (u.x - big);
-
-  double xx, y, x1, x2, c1, c2, res, cor;
-  double s, sn, ssn, c, cs, ccs;
-  xx = x * x;
-  s = x * xx * (sn3 + xx * sn5);
-  c = xx * (cs2 + xx * (cs4 + xx * cs6));
-  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
-  x1 = (x + t22) - t22;
-  x2 = (x - x1) + dx;
-  c1 = (cs + t22) - t22;
-  c2 = (cs - c1) + ccs;
-  cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c;
-  y = sn + c1 * x1;
-  cor = cor + ((sn - y) + c1 * x1);
-  res = y + cor;
-  cor = (y - res) + cor;
-  cor = 1.0005 * cor + __copysign (eps, cor);
-  *corp = cor;
-  return res;
-}
-
-/* Reduce range of X and compute sin of a + da. When SHIFT_QUADRANT is true,
-   the routine returns the cosine of a + da by rotating the quadrant once and
-   computing the sine of the result.  */
-static inline double
-__always_inline
-reduce_and_compute (double x, bool shift_quadrant)
-{
-  double retval = 0, a, da;
-  unsigned int n = __branred (x, &a, &da);
-  int4 k = (n + shift_quadrant) % 4;
-  switch (k)
-    {
-    case 2:
-      a = -a;
-      da = -da;
-      /* Fall through.  */
-    case 0:
-      if (a * a < 0.01588)
-	retval = bsloww (a, da, x, n);
-      else
-	retval = bsloww1 (a, da, x, n);
-      break;
-
-    case 1:
-    case 3:
-      retval = bsloww2 (a, da, x, n);
-      break;
-    }
-  return retval;
+  return __copysign (sn + cor, xold);
 }
 
+/* Reduce range of x to within PI/2 with abs (x) < 105414350.  The high part
+   is written to *a, the low part to *da.  Range reduction is accurate to 136
+   bits so that when x is large and *a very close to zero, all 53 bits of *a
+   are correct.  */
 static inline int4
 __always_inline
-reduce_sincos_1 (double x, double *a, double *da)
+reduce_sincos (double x, double *a, double *da)
 {
   mynumber v;
 
@@ -305,156 +159,54 @@ reduce_sincos_1 (double x, double *a, do
   v.x = t;
   double y = (x - xn * mp1) - xn * mp2;
   int4 n = v.i[LOW_HALF] & 3;
-  double db = xn * mp3;
-  double b = y - db;
-  db = (y - b) - db;
 
-  *a = b;
-  *da = db;
-
-  return n;
-}
-
-/* Compute sin (A + DA).  cos can be computed by passing SHIFT_QUADRANT as
-   true, which results in shifting the quadrant N clockwise.  */
-static double
-__always_inline
-do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant)
-{
-  double xx, retval, res, cor;
-  double eps = fabs (x) * 1.2e-30;
-
-  int k1 = (n + shift_quadrant) & 3;
-  switch (k1)
-    {			/* quarter of unit circle */
-    case 2:
-      a = -a;
-      da = -da;
-      /* Fall through.  */
-    case 0:
-      xx = a * a;
-      if (xx < 0.01588)
-	{
-	  /* Taylor series.  */
-	  res = TAYLOR_SIN (xx, a, da, cor);
-	  cor = 1.02 * cor + __copysign (eps, cor);
-	  retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant);
-	}
-      else
-	{
-	  res = do_sin (a, da, &cor);
-	  cor = 1.035 * cor + __copysign (eps, cor);
-	  retval = ((res == res + cor) ? __copysign (res, a)
-		    : sloww1 (a, da, x, shift_quadrant));
-	}
-      break;
-
-    case 1:
-    case 3:
-      res = do_cos (a, da, &cor);
-      cor = 1.025 * cor + __copysign (eps, cor);
-      retval = ((res == res + cor) ? ((n & 2) ? -res : res)
-		: sloww2 (a, da, x, n));
-      break;
-    }
-
-  return retval;
-}
-
-static inline int4
-__always_inline
-reduce_sincos_2 (double x, double *a, double *da)
-{
-  mynumber v;
-
-  double t = (x * hpinv + toint);
-  double xn = t - toint;
-  v.x = t;
-  double xn1 = (xn + 8.0e22) - 8.0e22;
-  double xn2 = xn - xn1;
-  double y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
-  int4 n = v.i[LOW_HALF] & 3;
-  double db = xn1 * pp3;
-  t = y - db;
-  db = (y - t) - db;
-  db = (db - xn2 * pp3) - xn * pp4;
-  double b = t + db;
-  db = (t - b) + db;
+  double b, db, t1, t2;
+  t1 = xn * pp3;
+  t2 = y - t1;
+  db = (y - t2) - t1;
+
+  t1 = xn * pp4;
+  b = t2 - t1;
+  db += (t2 - b) - t1;
 
   *a = b;
   *da = db;
-
   return n;
 }
 
-/* Compute sin (A + DA).  cos can be computed by passing SHIFT_QUADRANT as
-   true, which results in shifting the quadrant N clockwise.  */
+/* Compute sin or cos (A + DA) for the given quadrant N.  */
 static double
 __always_inline
-do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant)
+do_sincos (double a, double da, int4 n)
 {
-  double res, retval, cor, xx;
+  double retval;
 
-  double eps = 1.0e-24;
-
-  int4 k = (n + shift_quadrant) & 3;
-
-  switch (k)
-    {
-    case 2:
-      a = -a;
-      da = -da;
-      /* Fall through.  */
-    case 0:
-      xx = a * a;
-      if (xx < 0.01588)
-	{
-	  /* Taylor series.  */
-	  res = TAYLOR_SIN (xx, a, da, cor);
-	  cor = 1.02 * cor + __copysign (eps, cor);
-	  retval = (res == res + cor) ? res : bsloww (a, da, x, n);
-	}
-      else
-	{
-	  res = do_sin (a, da, &cor);
-	  cor = 1.035 * cor + __copysign (eps, cor);
-	  retval = ((res == res + cor) ? __copysign (res, a)
-		    : bsloww1 (a, da, x, n));
-	}
-      break;
-
-    case 1:
-    case 3:
-      res = do_cos (a, da, &cor);
-      cor = 1.025 * cor + __copysign (eps, cor);
-      retval = ((res == res + cor) ? ((n & 2) ? -res : res)
-		: bsloww2 (a, da, x, n));
-      break;
-    }
+  if (n & 1)
+    /* Max ULP is 0.513.  */
+    retval = do_cos (a, da);
+  else
+    /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518.  */
+    retval = do_sin (a, da);
 
-  return retval;
+  return (n & 2) ? -retval : retval;
 }
 
+
 /*******************************************************************/
 /* An ultimate sin routine. Given an IEEE double machine number x   */
 /* it computes the correctly rounded (to nearest) value of sin(x)  */
 /*******************************************************************/
-#ifdef IN_SINCOS
-static double
-#else
+#ifndef IN_SINCOS
 double
 SECTION
-#endif
 __sin (double x)
 {
-  double xx, res, t, cor;
+  double t, a, da;
   mynumber u;
-  int4 k, m;
+  int4 k, m, n;
   double retval = 0;
 
-#ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
 
   u.x = x;
   m = u.i[HIGH_HALF];
@@ -464,56 +216,34 @@ __sin (double x)
       math_check_force_underflow (x);
       retval = x;
     }
- /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
-  else if (k < 0x3fd00000)
-    {
-      xx = x * x;
-      /* Taylor series.  */
-      t = POLYNOMIAL (xx) * (xx * x);
-      res = x + t;
-      cor = (x - res) + t;
-      retval = (res == res + 1.07 * cor) ? res : slow (x);
-    }				/*  else  if (k < 0x3fd00000)    */
-/*---------------------------- 0.25<|x|< 0.855469---------------------- */
+/*--------------------------- 2^-26<|x|< 0.855469---------------------- */
   else if (k < 0x3feb6000)
     {
-      res = do_sin (x, 0, &cor);
-      retval = (res == res + 1.096 * cor) ? res : slow1 (x);
-      retval = __copysign (retval, x);
+      /* Max ULP is 0.548.  */
+      retval = do_sin (x, 0);
     }				/*   else  if (k < 0x3feb6000)    */
 
 /*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
   else if (k < 0x400368fd)
     {
-
       t = hp0 - fabs (x);
-      res = do_cos (t, hp1, &cor);
-      retval = (res == res + 1.020 * cor) ? res : slow2 (x);
-      retval = __copysign (retval, x);
+      /* Max ULP is 0.51.  */
+      retval = __copysign (do_cos (t, hp1), x);
     }				/*   else  if (k < 0x400368fd)    */
 
-#ifndef IN_SINCOS
 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
   else if (k < 0x419921FB)
     {
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, false);
+      n = reduce_sincos (x, &a, &da);
+      retval = do_sincos (a, da, n);
     }				/*   else  if (k <  0x419921FB )    */
 
-/*---------------------105414350 <|x|< 281474976710656 --------------------*/
-  else if (k < 0x42F00000)
-    {
-      double a, da;
-
-      int4 n = reduce_sincos_2 (x, &a, &da);
-      retval = do_sincos_2 (a, da, x, n, false);
-    }				/*   else  if (k <  0x42F00000 )   */
-
-/* -----------------281474976710656 <|x| <2^1024----------------------------*/
+/* --------------------105414350 <|x| <2^1024------------------------------*/
   else if (k < 0x7ff00000)
-    retval = reduce_and_compute (x, false);
-
+    {
+      n = __branred (x, &a, &da);
+      retval = do_sincos (a, da, n);
+    }
 /*--------------------- |x| > 2^1024 ----------------------------------*/
   else
     {
@@ -521,7 +251,6 @@ __sin (double x)
 	__set_errno (EDOM);
       retval = x / x;
     }
-#endif
 
   return retval;
 }
@@ -532,23 +261,17 @@ __sin (double x)
 /* it computes the correctly rounded (to nearest) value of cos(x)  */
 /*******************************************************************/
 
-#ifdef IN_SINCOS
-static double
-#else
 double
 SECTION
-#endif
 __cos (double x)
 {
-  double y, xx, res, cor, a, da;
+  double y, a, da;
   mynumber u;
-  int4 k, m;
+  int4 k, m, n;
 
   double retval = 0;
 
-#ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
 
   u.x = x;
   m = u.i[HIGH_HALF];
@@ -560,8 +283,8 @@ __cos (double x)
 
   else if (k < 0x3feb6000)
     {				/* 2^-27 < |x| < 0.855469 */
-      res = do_cos (x, 0, &cor);
-      retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
+      /* Max ULP is 0.51.  */
+      retval = do_cos (x, 0);
     }				/*   else  if (k < 0x3feb6000)    */
 
   else if (k < 0x400368fd)
@@ -569,43 +292,23 @@ __cos (double x)
       y = hp0 - fabs (x);
       a = y + hp1;
       da = (y - a) + hp1;
-      xx = a * a;
-      if (xx < 0.01588)
-	{
-	  res = TAYLOR_SIN (xx, a, da, cor);
-	  cor = 1.02 * cor + __copysign (1.0e-31, cor);
-	  retval = (res == res + cor) ? res : sloww (a, da, x, true);
-	}
-      else
-	{
-	  res = do_sin (a, da, &cor);
-	  cor = 1.035 * cor + __copysign (1.0e-31, cor);
-	  retval = ((res == res + cor) ? __copysign (res, a)
-		    : sloww1 (a, da, x, true));
-	}
-
+      /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise.
+	 Range reduction uses 106 bits here which is sufficient.  */
+      retval = do_sin (a, da);
     }				/*   else  if (k < 0x400368fd)    */
 
-
-#ifndef IN_SINCOS
   else if (k < 0x419921FB)
     {				/* 2.426265<|x|< 105414350 */
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, true);
+      n = reduce_sincos (x, &a, &da);
+      retval = do_sincos (a, da, n + 1);
     }				/*   else  if (k <  0x419921FB )    */
 
-  else if (k < 0x42F00000)
-    {
-      double a, da;
-
-      int4 n = reduce_sincos_2 (x, &a, &da);
-      retval = do_sincos_2 (a, da, x, n, true);
-    }				/*   else  if (k <  0x42F00000 )    */
-
-  /* 281474976710656 <|x| <2^1024 */
+  /* 105414350 <|x| <2^1024 */
   else if (k < 0x7ff00000)
-    retval = reduce_and_compute (x, true);
+    {
+      n = __branred (x, &a, &da);
+      retval = do_sincos (a, da, n + 1);
+    }
 
   else
     {
@@ -613,304 +316,10 @@ __cos (double x)
 	__set_errno (EDOM);
       retval = x / x;		/* |x| > 2^1024 */
     }
-#endif
 
   return retval;
 }
 
-/************************************************************************/
-/*  Routine compute sin(x) for  2^-26 < |x|< 0.25 by  Taylor with more   */
-/* precision  and if still doesn't accurate enough by mpsin   or dubsin */
-/************************************************************************/
-
-static inline double
-__always_inline
-slow (double x)
-{
-  double res, cor, w[2];
-  res = TAYLOR_SLOW (x, 0, cor);
-  if (res == res + 1.0007 * cor)
-    return res;
-
-  __dubsin (fabs (x), 0, w);
-  if (w[0] == w[0] + 1.000000001 * w[1])
-    return __copysign (w[0], x);
-
-  return __copysign (__mpsin (fabs (x), 0, false), x);
-}
-
-/*******************************************************************************/
-/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
-/* and if result still doesn't accurate enough by mpsin   or dubsin            */
-/*******************************************************************************/
-
-static inline double
-__always_inline
-slow1 (double x)
-{
-  double w[2], cor, res;
-
-  res = do_sin_slow (x, 0, 0, &cor);
-  if (res == res + cor)
-    return res;
-
-  __dubsin (fabs (x), 0, w);
-  if (w[0] == w[0] + 1.000000005 * w[1])
-    return w[0];
-
-  return __mpsin (fabs (x), 0, false);
-}
-
-/**************************************************************************/
-/*  Routine compute sin(x) for   0.855469  <|x|<2.426265  by  __sincostab.tbl  */
-/* and if result still doesn't accurate enough by mpsin   or dubsin       */
-/**************************************************************************/
-static inline double
-__always_inline
-slow2 (double x)
-{
-  double w[2], y, y1, y2, cor, res;
-
-  double t = hp0 - fabs (x);
-  res = do_cos_slow (t, hp1, 0, &cor);
-  if (res == res + cor)
-    return res;
-
-  y = fabs (x) - hp0;
-  y1 = y - hp1;
-  y2 = (y - y1) - hp1;
-  __docos (y1, y2, w);
-  if (w[0] == w[0] + 1.000000005 * w[1])
-    return w[0];
-
-  return __mpsin (fabs (x), 0, false);
-}
-
-/* Compute sin(x + dx) where X is small enough to use Taylor series around zero
-   and (x + dx) in the first or third quarter of the unit circle.  ORIG is the
-   original value of X for computing error of the result.  If the result is not
-   accurate enough, the routine calls mpsin or dubsin.  SHIFT_QUADRANT rotates
-   the unit circle by 1 to compute the cosine instead of sine.  */
-static inline double
-__always_inline
-sloww (double x, double dx, double orig, bool shift_quadrant)
-{
-  double y, t, res, cor, w[2], a, da, xn;
-  mynumber v;
-  int4 n;
-  res = TAYLOR_SLOW (x, dx, cor);
-
-  double eps = fabs (orig) * 3.1e-30;
-
-  cor = 1.0005 * cor + __copysign (eps, cor);
-
-  if (res == res + cor)
-    return res;
-
-  a = fabs (x);
-  da = (x > 0) ? dx : -dx;
-  __dubsin (a, da, w);
-  eps = fabs (orig) * 1.1e-30;
-  cor = 1.000000001 * w[1] + __copysign (eps, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return __copysign (w[0], x);
-
-  t = (orig * hpinv + toint);
-  xn = t - toint;
-  v.x = t;
-  y = (orig - xn * mp1) - xn * mp2;
-  n = (v.i[LOW_HALF] + shift_quadrant) & 3;
-  da = xn * pp3;
-  t = y - da;
-  da = (y - t) - da;
-  y = xn * pp4;
-  a = t - y;
-  da = ((t - a) - y) + da;
-
-  if (n & 2)
-    {
-      a = -a;
-      da = -da;
-    }
-  x = fabs (a);
-  dx = (a > 0) ? da : -da;
-  __dubsin (x, dx, w);
-  eps = fabs (orig) * 1.1e-40;
-  cor = 1.000000001 * w[1] + __copysign (eps, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return __copysign (w[0], a);
-
-  return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
-}
-
-/* Compute sin(x + dx) where X is in the first or third quarter of the unit
-   circle.  ORIG is the original value of X for computing error of the result.
-   If the result is not accurate enough, the routine calls mpsin or dubsin.
-   SHIFT_QUADRANT rotates the unit circle by 1 to compute the cosine instead of
-   sine.  */
-static inline double
-__always_inline
-sloww1 (double x, double dx, double orig, bool shift_quadrant)
-{
-  double w[2], cor, res;
-
-  res = do_sin_slow (x, dx, 3.1e-30 * fabs (orig), &cor);
-
-  if (res == res + cor)
-    return __copysign (res, x);
-
-  dx = (x > 0 ? dx : -dx);
-  __dubsin (fabs (x), dx, w);
-
-  double eps = 1.1e-30 * fabs (orig);
-  cor = 1.000000005 * w[1] + __copysign (eps, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return __copysign (w[0], x);
-
-  return shift_quadrant ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
-}
-
-/***************************************************************************/
-/*  Routine compute sin(x+dx)   (Double-Length number) where x in second or */
-/*  fourth quarter of unit circle.Routine receive also  the  original value */
-/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
-/* accurate enough routine calls  mpsin1   or dubsin                       */
-/***************************************************************************/
-
-static inline double
-__always_inline
-sloww2 (double x, double dx, double orig, int n)
-{
-  double w[2], cor, res;
-
-  res = do_cos_slow (x, dx, 3.1e-30 * fabs (orig), &cor);
-
-  if (res == res + cor)
-    return (n & 2) ? -res : res;
-
-  dx = x > 0 ? dx : -dx;
-  __docos (fabs (x), dx, w);
-
-  double eps = 1.1e-30 * fabs (orig);
-  cor = 1.000000005 * w[1] + __copysign (eps, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return (n & 2) ? -w[0] : w[0];
-
-  return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
-}
-
-/***************************************************************************/
-/*  Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x   */
-/* is small enough to use Taylor series around zero and   (x+dx)            */
-/* in first or third quarter of unit circle.Routine receive also            */
-/* (right argument) the  original   value of x for computing error of      */
-/* result.And if result not accurate enough routine calls other routines    */
-/***************************************************************************/
-
-static inline double
-__always_inline
-bsloww (double x, double dx, double orig, int n)
-{
-  double res, cor, w[2], a, da;
-
-  res = TAYLOR_SLOW (x, dx, cor);
-  cor = 1.0005 * cor + __copysign (1.1e-24, cor);
-  if (res == res + cor)
-    return res;
-
-  a = fabs (x);
-  da = (x > 0) ? dx : -dx;
-  __dubsin (a, da, w);
-  cor = 1.000000001 * w[1] + __copysign (1.1e-24, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return __copysign (w[0], x);
-
-  return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
-}
-
-/***************************************************************************/
-/*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
-/* in first or third quarter of unit circle.Routine receive also            */
-/* (right argument) the original  value of x for computing error of result.*/
-/* And if result not  accurate enough routine calls  other routines         */
-/***************************************************************************/
-
-static inline double
-__always_inline
-bsloww1 (double x, double dx, double orig, int n)
-{
-  double w[2], cor, res;
-
-  res = do_sin_slow (x, dx, 1.1e-24, &cor);
-  if (res == res + cor)
-    return (x > 0) ? res : -res;
-
-  dx = (x > 0) ? dx : -dx;
-  __dubsin (fabs (x), dx, w);
-
-  cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return __copysign (w[0], x);
-
-  return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
-}
-
-/***************************************************************************/
-/*  Routine compute sin(x+dx)  or cos(x+dx) (Double-Length number) where x  */
-/* in second or fourth quarter of unit circle.Routine receive also  the     */
-/* original value and quarter(n= 1or 3)of x for computing error of result.  */
-/* And if result not accurate enough routine calls  other routines          */
-/***************************************************************************/
-
-static inline double
-__always_inline
-bsloww2 (double x, double dx, double orig, int n)
-{
-  double w[2], cor, res;
-
-  res = do_cos_slow (x, dx, 1.1e-24, &cor);
-  if (res == res + cor)
-    return (n & 2) ? -res : res;
-
-  dx = (x > 0) ? dx : -dx;
-  __docos (fabs (x), dx, w);
-
-  cor = 1.000000005 * w[1] + __copysign (1.1e-24, w[1]);
-
-  if (w[0] == w[0] + cor)
-    return (n & 2) ? -w[0] : w[0];
-
-  return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
-}
-
-/************************************************************************/
-/*  Routine compute cos(x) for  2^-27 < |x|< 0.25 by  Taylor with more   */
-/* precision  and if still doesn't accurate enough by mpcos   or docos  */
-/************************************************************************/
-
-static inline double
-__always_inline
-cslow2 (double x)
-{
-  double w[2], cor, res;
-
-  res = do_cos_slow (x, 0, 0, &cor);
-  if (res == res + cor)
-    return res;
-
-  __docos (fabs (x), 0, w);
-  if (w[0] == w[0] + 1.000000005 * w[1])
-    return w[0];
-
-  return __mpcos (x, 0, false);
-}
-
 #ifndef __cos
 weak_alias (__cos, cos)
 # ifdef NO_LONG_DOUBLE
@@ -925,3 +334,5 @@ strong_alias (__sin, __sinl)
 weak_alias (__sin, sinl)
 # endif
 #endif
+
+#endif
Index: glibc-2.26/sysdeps/ieee754/dbl-64/s_sincos.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/s_sincos.c
+++ glibc-2.26/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -22,42 +22,9 @@
 
 #include <math_private.h>
 
-#define __sin __sin_local
-#define __cos __cos_local
-#define IN_SINCOS 1
+#define IN_SINCOS
 #include "s_sin.c"
 
-/* Consolidated version of reduce_and_compute in s_sin.c that does range
-   reduction only once and computes sin and cos together.  */
-static inline void
-__always_inline
-reduce_and_compute_sincos (double x, double *sinx, double *cosx)
-{
-  double a, da;
-  unsigned int n = __branred (x, &a, &da);
-
-  n = n & 3;
-
-  if (n == 1 || n == 2)
-    {
-      a = -a;
-      da = -da;
-    }
-
-  if (n & 1)
-    {
-      double *temp = cosx;
-      cosx = sinx;
-      sinx = temp;
-    }
-
-  if (a * a < 0.01588)
-    *sinx = bsloww (a, da, x, n);
-  else
-    *sinx = bsloww1 (a, da, x, n);
-  *cosx = bsloww2 (a, da, x, n);
-}
-
 void
 __sincos (double x, double *sinx, double *cosx)
 {
@@ -67,37 +34,62 @@ __sincos (double x, double *sinx, double
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 
   u.x = x;
-  k = 0x7fffffff & u.i[HIGH_HALF];
+  k = u.i[HIGH_HALF] & 0x7fffffff;
 
   if (k < 0x400368fd)
     {
-      *sinx = __sin_local (x);
-      *cosx = __cos_local (x);
-      return;
-    }
-  if (k < 0x419921FB)
-    {
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-
-      *sinx = do_sincos_1 (a, da, x, n, false);
-      *cosx = do_sincos_1 (a, da, x, n, true);
-
-      return;
-    }
-  if (k < 0x42F00000)
-    {
-      double a, da;
-      int4 n = reduce_sincos_2 (x, &a, &da);
-
-      *sinx = do_sincos_2 (a, da, x, n, false);
-      *cosx = do_sincos_2 (a, da, x, n, true);
-
+      double a, da, y;
+      /* |x| < 2^-27 => cos (x) = 1, sin (x) = x.  */
+      if (k < 0x3e400000)
+	{
+	  if (k < 0x3e500000)
+	    math_check_force_underflow (x);
+	  *sinx = x;
+	  *cosx = 1.0;
+	  return;
+	}
+      /* |x| < 0.855469.  */
+      else if (k < 0x3feb6000)
+	{
+	  *sinx = do_sin (x, 0);
+	  *cosx = do_cos (x, 0);
+	  return;
+	}
+
+      /* |x| < 2.426265.  */
+      y = hp0 - fabs (x);
+      a = y + hp1;
+      da = (y - a) + hp1;
+      *sinx = __copysign (do_cos (a, da), x);
+      *cosx = do_sin (a, da);
       return;
     }
+  /* |x| < 2^1024.  */
   if (k < 0x7ff00000)
     {
-      reduce_and_compute_sincos (x, sinx, cosx);
+      double a, da, xx;
+      unsigned int n;
+
+      /* If |x| < 105414350 use simple range reduction.  */
+      n = k < 0x419921FB ? reduce_sincos (x, &a, &da) : __branred (x, &a, &da);
+      n = n & 3;
+
+      if (n == 1 || n == 2)
+	{
+	  a = -a;
+	  da = -da;
+	}
+
+      if (n & 1)
+	{
+	  double *temp = cosx;
+	  cosx = sinx;
+	  sinx = temp;
+	}
+
+      *sinx = do_sin (a, da);
+      xx = do_cos (a, da);
+      *cosx = (n & 2) ? -xx : xx;
       return;
     }
 
Index: glibc-2.26/sysdeps/ieee754/dbl-64/slowexp.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/slowexp.c
+++ /dev/null
@@ -1,86 +0,0 @@
-/*
- * 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:slowexp.c                                                 */
-/*                                                                        */
-/*  FUNCTION:slowexp                                                      */
-/*                                                                        */
-/*  FILES NEEDED:mpa.h                                                    */
-/*               mpa.c mpexp.c                                            */
-/*                                                                        */
-/*Converting from double precision to Multi-precision and calculating     */
-/* e^x                                                                    */
-/**************************************************************************/
-#include <math_private.h>
-
-#include <stap-probe.h>
-
-#ifndef USE_LONG_DOUBLE_FOR_MP
-# include "mpa.h"
-void __mpexp (mp_no *x, mp_no *y, int p);
-#endif
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-/*Converting from double precision to Multi-precision and calculating  e^x */
-double
-SECTION
-__slowexp (double x)
-{
-#ifndef USE_LONG_DOUBLE_FOR_MP
-  double w, z, res, eps = 3.0e-26;
-  int p;
-  mp_no mpx, mpy, mpz, mpw, mpeps, mpcor;
-
-  /* Use the multiple precision __MPEXP function to compute the exponential
-     First at 144 bits and if it is not accurate enough, at 768 bits.  */
-  p = 6;
-  __dbl_mp (x, &mpx, p);
-  __mpexp (&mpx, &mpy, p);
-  __dbl_mp (eps, &mpeps, p);
-  __mul (&mpeps, &mpy, &mpcor, p);
-  __add (&mpy, &mpcor, &mpw, p);
-  __sub (&mpy, &mpcor, &mpz, p);
-  __mp_dbl (&mpw, &w, p);
-  __mp_dbl (&mpz, &z, p);
-  if (w == z)
-    {
-      /* Track how often we get to the slow exp code plus
-	 its input/output values.  */
-      LIBC_PROBE (slowexp_p6, 2, &x, &w);
-      return w;
-    }
-  else
-    {
-      p = 32;
-      __dbl_mp (x, &mpx, p);
-      __mpexp (&mpx, &mpy, p);
-      __mp_dbl (&mpy, &res, p);
-
-      /* Track how often we get to the uber-slow exp code plus
-	 its input/output values.  */
-      LIBC_PROBE (slowexp_p32, 2, &x, &res);
-      return res;
-    }
-#else
-  return (double) __ieee754_expl((long double)x);
-#endif
-}
Index: glibc-2.26/sysdeps/ieee754/dbl-64/slowpow.c
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/slowpow.c
+++ /dev/null
@@ -1,125 +0,0 @@
-/*
- * 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:slowpow.c                                                 */
-/*                                                                       */
-/* FUNCTION:slowpow                                                      */
-/*                                                                       */
-/*FILES NEEDED:mpa.h                                                     */
-/*             mpa.c mpexp.c mplog.c halfulp.c                           */
-/*                                                                       */
-/* Given two IEEE double machine numbers y,x , routine  computes the     */
-/* correctly  rounded (to nearest) value of x^y. Result calculated  by   */
-/* multiplication (in halfulp.c) or if result isn't accurate enough      */
-/* then routine converts x and y into multi-precision doubles     and    */
-/* calls to mpexp routine                                                */
-/*************************************************************************/
-
-#include "mpa.h"
-#include <math_private.h>
-
-#include <stap-probe.h>
-
-#ifndef SECTION
-# define SECTION
-#endif
-
-void __mpexp (mp_no *x, mp_no *y, int p);
-void __mplog (mp_no *x, mp_no *y, int p);
-double ulog (double);
-double __halfulp (double x, double y);
-
-double
-SECTION
-__slowpow (double x, double y, double z)
-{
-  double res, res1;
-  mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1;
-  static const mp_no eps = {-3, {1.0, 4.0}};
-  int p;
-
-  /* __HALFULP returns -10 or X^Y.  */
-  res = __halfulp (x, y);
-
-  /* Return if the result was computed by __HALFULP.  */
-  if (res >= 0)
-    return res;
-
-  /* Compute pow as long double.  This is currently only used by powerpc, where
-     one may get 106 bits of accuracy.  */
-#ifdef USE_LONG_DOUBLE_FOR_MP
-  long double ldw, ldz, ldpp;
-  static const long double ldeps = 0x4.0p-96;
-
-  ldz = __ieee754_logl ((long double) x);
-  ldw = (long double) y *ldz;
-  ldpp = __ieee754_expl (ldw);
-  res = (double) (ldpp + ldeps);
-  res1 = (double) (ldpp - ldeps);
-
-  /* Return the result if it is accurate enough.  */
-  if (res == res1)
-    return res;
-#endif
-
-  /* Or else, calculate using multiple precision.  P = 10 implies accuracy of
-     240 bits accuracy, since MP_NO has a radix of 2^24.  */
-  p = 10;
-  __dbl_mp (x, &mpx, p);
-  __dbl_mp (y, &mpy, p);
-  __dbl_mp (z, &mpz, p);
-
-  /* z = x ^ y
-     log (z) = y * log (x)
-     z = exp (y * log (x))  */
-  __mplog (&mpx, &mpz, p);
-  __mul (&mpy, &mpz, &mpw, p);
-  __mpexp (&mpw, &mpp, p);
-
-  /* Add and subtract EPS to ensure that the result remains unchanged, i.e. we
-     have last bit accuracy.  */
-  __add (&mpp, &eps, &mpr, p);
-  __mp_dbl (&mpr, &res, p);
-  __sub (&mpp, &eps, &mpr1, p);
-  __mp_dbl (&mpr1, &res1, p);
-  if (res == res1)
-    {
-      /* Track how often we get to the slow pow code plus
-	 its input/output values.  */
-      LIBC_PROBE (slowpow_p10, 4, &x, &y, &z, &res);
-      return res;
-    }
-
-  /* If we don't, then we repeat using a higher precision.  768 bits of
-     precision ought to be enough for anybody.  */
-  p = 32;
-  __dbl_mp (x, &mpx, p);
-  __dbl_mp (y, &mpy, p);
-  __dbl_mp (z, &mpz, p);
-  __mplog (&mpx, &mpz, p);
-  __mul (&mpy, &mpz, &mpw, p);
-  __mpexp (&mpw, &mpp, p);
-  __mp_dbl (&mpp, &res, p);
-
-  /* Track how often we get to the uber-slow pow code plus
-     its input/output values.  */
-  LIBC_PROBE (slowpow_p32, 4, &x, &y, &z, &res);
-
-  return res;
-}
Index: glibc-2.26/sysdeps/ieee754/dbl-64/uexp.h
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/uexp.h
+++ glibc-2.26/sysdeps/ieee754/dbl-64/uexp.h
@@ -29,8 +29,7 @@
 
 #include "mydefs.h"
 
-const static double zero = 0.0, hhuge = 1.0e300, tiny = 1.0e-300,
-err_0 = 1.000014, err_1 = 0.000016;
+const static double zero = 0.0, hhuge = 1.0e300, tiny = 1.0e-300;
 const static int4 bigint = 0x40862002,
              badint = 0x40876000,smallint = 0x3C8fffff;
 const static int4 hugeint = 0x7FFFFFFF, infint = 0x7ff00000;
Index: glibc-2.26/sysdeps/ieee754/dbl-64/ulog.h
===================================================================
--- glibc-2.26.orig/sysdeps/ieee754/dbl-64/ulog.h
+++ glibc-2.26/sysdeps/ieee754/dbl-64/ulog.h
@@ -42,43 +42,6 @@
 /**/ b6             = {{0x3fbc71c5, 0x25db58ac} }, /*  0.111... */
 /**/ b7             = {{0xbfb9a4ac, 0x11a2a61c} }, /* -0.100... */
 /**/ b8             = {{0x3fb75077, 0x0df2b591} }, /*  0.091... */
-  /* polynomial III */
-#if 0
-/**/ c1             = {{0x3ff00000, 0x00000000} }, /*  1        */
-#endif
-/**/ c2             = {{0xbfe00000, 0x00000000} }, /* -1/2      */
-/**/ c3             = {{0x3fd55555, 0x55555555} }, /*  1/3      */
-/**/ c4             = {{0xbfd00000, 0x00000000} }, /* -1/4      */
-/**/ c5             = {{0x3fc99999, 0x9999999a} }, /*  1/5      */
-  /* polynomial IV */
-/**/ d2             = {{0xbfe00000, 0x00000000} }, /* -1/2      */
-/**/ dd2            = {{0x00000000, 0x00000000} }, /* -1/2-d2   */
-/**/ d3             = {{0x3fd55555, 0x55555555} }, /*  1/3      */
-/**/ dd3            = {{0x3c755555, 0x55555555} }, /*  1/3-d3   */
-/**/ d4             = {{0xbfd00000, 0x00000000} }, /* -1/4      */
-/**/ dd4            = {{0x00000000, 0x00000000} }, /* -1/4-d4   */
-/**/ d5             = {{0x3fc99999, 0x9999999a} }, /*  1/5      */
-/**/ dd5            = {{0xbc699999, 0x9999999a} }, /*  1/5-d5   */
-/**/ d6             = {{0xbfc55555, 0x55555555} }, /* -1/6      */
-/**/ dd6            = {{0xbc655555, 0x55555555} }, /* -1/6-d6   */
-/**/ d7             = {{0x3fc24924, 0x92492492} }, /*  1/7      */
-/**/ dd7            = {{0x3c624924, 0x92492492} }, /*  1/7-d7   */
-/**/ d8             = {{0xbfc00000, 0x00000000} }, /* -1/8      */
-/**/ dd8            = {{0x00000000, 0x00000000} }, /* -1/8-d8   */
-/**/ d9             = {{0x3fbc71c7, 0x1c71c71c} }, /*  1/9      */
-/**/ dd9            = {{0x3c5c71c7, 0x1c71c71c} }, /*  1/9-d9   */
-/**/ d10            = {{0xbfb99999, 0x9999999a} }, /* -1/10     */
-/**/ dd10           = {{0x3c599999, 0x9999999a} }, /* -1/10-d10 */
-/**/ d11            = {{0x3fb745d1, 0x745d1746} }, /*  1/11     */
-/**/ d12            = {{0xbfb55555, 0x55555555} }, /* -1/12     */
-/**/ d13            = {{0x3fb3b13b, 0x13b13b14} }, /*  1/13     */
-/**/ d14            = {{0xbfb24924, 0x92492492} }, /* -1/14     */
-/**/ d15            = {{0x3fb11111, 0x11111111} }, /*  1/15     */
-/**/ d16            = {{0xbfb00000, 0x00000000} }, /* -1/16     */
-/**/ d17            = {{0x3fae1e1e, 0x1e1e1e1e} }, /*  1/17     */
-/**/ d18            = {{0xbfac71c7, 0x1c71c71c} }, /* -1/18     */
-/**/ d19            = {{0x3faaf286, 0xbca1af28} }, /*  1/19     */
-/**/ d20            = {{0xbfa99999, 0x9999999a} }, /* -1/20     */
   /* constants    */
 /**/ sqrt_2         = {{0x3ff6a09e, 0x667f3bcc} }, /* sqrt(2)   */
 /**/ h1             = {{0x3fd2e000, 0x00000000} }, /* 151/2**9  */
@@ -87,14 +50,6 @@
 /**/ delv           = {{0x3ef00000, 0x00000000} }, /* 1/2**16   */
 /**/ ln2a           = {{0x3fe62e42, 0xfefa3800} }, /* ln(2) 43 bits */
 /**/ ln2b           = {{0x3d2ef357, 0x93c76730} }, /* ln(2)-ln2a    */
-/**/ e1             = {{0x3bbcc868, 0x00000000} }, /* 6.095e-21     */
-/**/ e2             = {{0x3c1138ce, 0x00000000} }, /* 2.334e-19     */
-/**/ e3             = {{0x3aa1565d, 0x00000000} }, /* 2.801e-26     */
-/**/ e4             = {{0x39809d88, 0x00000000} }, /* 1.024e-31     */
-/**/ e[M]           ={{{0x37da223a, 0x00000000} }, /* 1.2e-39       */
-/**/                  {{0x35c851c4, 0x00000000} }, /* 1.3e-49       */
-/**/                  {{0x2ab85e51, 0x00000000} }, /* 6.8e-103      */
-/**/                  {{0x17383827, 0x00000000} }},/* 8.1e-197      */
 /**/ two54          = {{0x43500000, 0x00000000} }, /* 2**54         */
 /**/ u03            = {{0x3f9eb851, 0xeb851eb8} }; /* 0.03          */
 
@@ -114,43 +69,6 @@
 /**/ b6             = {{0x25db58ac, 0x3fbc71c5} }, /*  0.111... */
 /**/ b7             = {{0x11a2a61c, 0xbfb9a4ac} }, /* -0.100... */
 /**/ b8             = {{0x0df2b591, 0x3fb75077} }, /*  0.091... */
-  /* polynomial III */
-#if 0
-/**/ c1             = {{0x00000000, 0x3ff00000} }, /*  1        */
-#endif
-/**/ c2             = {{0x00000000, 0xbfe00000} }, /* -1/2      */
-/**/ c3             = {{0x55555555, 0x3fd55555} }, /*  1/3      */
-/**/ c4             = {{0x00000000, 0xbfd00000} }, /* -1/4      */
-/**/ c5             = {{0x9999999a, 0x3fc99999} }, /*  1/5      */
-  /* polynomial IV */
-/**/ d2             = {{0x00000000, 0xbfe00000} }, /* -1/2      */
-/**/ dd2            = {{0x00000000, 0x00000000} }, /* -1/2-d2   */
-/**/ d3             = {{0x55555555, 0x3fd55555} }, /*  1/3      */
-/**/ dd3            = {{0x55555555, 0x3c755555} }, /*  1/3-d3   */
-/**/ d4             = {{0x00000000, 0xbfd00000} }, /* -1/4      */
-/**/ dd4            = {{0x00000000, 0x00000000} }, /* -1/4-d4   */
-/**/ d5             = {{0x9999999a, 0x3fc99999} }, /*  1/5      */
-/**/ dd5            = {{0x9999999a, 0xbc699999} }, /*  1/5-d5   */
-/**/ d6             = {{0x55555555, 0xbfc55555} }, /* -1/6      */
-/**/ dd6            = {{0x55555555, 0xbc655555} }, /* -1/6-d6   */
-/**/ d7             = {{0x92492492, 0x3fc24924} }, /*  1/7      */
-/**/ dd7            = {{0x92492492, 0x3c624924} }, /*  1/7-d7   */
-/**/ d8             = {{0x00000000, 0xbfc00000} }, /* -1/8      */
-/**/ dd8            = {{0x00000000, 0x00000000} }, /* -1/8-d8   */
-/**/ d9             = {{0x1c71c71c, 0x3fbc71c7} }, /*  1/9      */
-/**/ dd9            = {{0x1c71c71c, 0x3c5c71c7} }, /*  1/9-d9   */
-/**/ d10            = {{0x9999999a, 0xbfb99999} }, /* -1/10     */
-/**/ dd10           = {{0x9999999a, 0x3c599999} }, /* -1/10-d10 */
-/**/ d11            = {{0x745d1746, 0x3fb745d1} }, /*  1/11     */
-/**/ d12            = {{0x55555555, 0xbfb55555} }, /* -1/12     */
-/**/ d13            = {{0x13b13b14, 0x3fb3b13b} }, /*  1/13     */
-/**/ d14            = {{0x92492492, 0xbfb24924} }, /* -1/14     */
-/**/ d15            = {{0x11111111, 0x3fb11111} }, /*  1/15     */
-/**/ d16            = {{0x00000000, 0xbfb00000} }, /* -1/16     */
-/**/ d17            = {{0x1e1e1e1e, 0x3fae1e1e} }, /*  1/17     */
-/**/ d18            = {{0x1c71c71c, 0xbfac71c7} }, /* -1/18     */
-/**/ d19            = {{0xbca1af28, 0x3faaf286} }, /*  1/19     */
-/**/ d20            = {{0x9999999a, 0xbfa99999} }, /* -1/20     */
   /* constants    */
 /**/ sqrt_2         = {{0x667f3bcc, 0x3ff6a09e} }, /* sqrt(2)   */
 /**/ h1             = {{0x00000000, 0x3fd2e000} }, /* 151/2**9  */
@@ -159,14 +77,6 @@
 /**/ delv           = {{0x00000000, 0x3ef00000} }, /* 1/2**16   */
 /**/ ln2a           = {{0xfefa3800, 0x3fe62e42} }, /* ln(2) 43 bits */
 /**/ ln2b           = {{0x93c76730, 0x3d2ef357} }, /* ln(2)-ln2a    */
-/**/ e1             = {{0x00000000, 0x3bbcc868} }, /* 6.095e-21     */
-/**/ e2             = {{0x00000000, 0x3c1138ce} }, /* 2.334e-19     */
-/**/ e3             = {{0x00000000, 0x3aa1565d} }, /* 2.801e-26     */
-/**/ e4             = {{0x00000000, 0x39809d88} }, /* 1.024e-31     */
-/**/ e[M]           ={{{0x00000000, 0x37da223a} }, /* 1.2e-39       */
-/**/                  {{0x00000000, 0x35c851c4} }, /* 1.3e-49       */
-/**/                  {{0x00000000, 0x2ab85e51} }, /* 6.8e-103      */
-/**/                  {{0x00000000, 0x17383827} }},/* 8.1e-197      */
 /**/ two54          = {{0x00000000, 0x43500000} }, /* 2**54         */
 /**/ u03            = {{0xeb851eb8, 0x3f9eb851} }; /* 0.03          */
 
@@ -178,10 +88,6 @@
 #define  DEL_V     delv.d
 #define  LN2A      ln2a.d
 #define  LN2B      ln2b.d
-#define  E1        e1.d
-#define  E2        e2.d
-#define  E3        e3.d
-#define  E4        e4.d
 #define  U03       u03.d
 
 #endif
Index: glibc-2.26/sysdeps/m68k/m680x0/fpu/halfulp.c
===================================================================
--- glibc-2.26.orig/sysdeps/m68k/m680x0/fpu/halfulp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/m68k/m680x0/fpu/slowexp.c
===================================================================
--- glibc-2.26.orig/sysdeps/m68k/m680x0/fpu/slowexp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/m68k/m680x0/fpu/slowpow.c
===================================================================
--- glibc-2.26.orig/sysdeps/m68k/m680x0/fpu/slowpow.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
Index: glibc-2.26/sysdeps/powerpc/power4/fpu/Makefile
===================================================================
--- glibc-2.26.orig/sysdeps/powerpc/power4/fpu/Makefile
+++ glibc-2.26/sysdeps/powerpc/power4/fpu/Makefile
@@ -2,6 +2,4 @@
 
 ifeq ($(subdir),math)
 CFLAGS-mpa.c += --param max-unroll-times=4 -funroll-loops -fpeel-loops
-CPPFLAGS-slowpow.c += -DUSE_LONG_DOUBLE_FOR_MP=1
-CPPFLAGS-slowexp.c += -DUSE_LONG_DOUBLE_FOR_MP=1
 endif
Index: glibc-2.26/sysdeps/x86_64/fpu/libm-test-ulps
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/libm-test-ulps
+++ glibc-2.26/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1262,7 +1262,9 @@ ildouble: 1
 ldouble: 1
 
 Function: "cos":
+double: 1
 float128: 1
+idouble: 1
 ifloat128: 1
 ildouble: 1
 ldouble: 1
@@ -2464,8 +2466,10 @@ Function: "log_vlen8_avx2":
 float: 2
 
 Function: "pow":
+double: 1
 float: 1
 float128: 2
+idouble: 1
 ifloat: 1
 ifloat128: 2
 ildouble: 1
@@ -2552,7 +2556,9 @@ Function: "pow_vlen8_avx2":
 float: 3
 
 Function: "sin":
+double: 1
 float128: 1
+idouble: 1
 ifloat128: 1
 ildouble: 1
 ldouble: 1
@@ -2602,7 +2608,9 @@ Function: "sin_vlen8_avx2":
 float: 1
 
 Function: "sincos":
+double: 1
 float128: 1
+idouble: 1
 ifloat128: 1
 ildouble: 1
 ldouble: 1
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/Makefile
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/Makefile
+++ glibc-2.26/sysdeps/x86_64/fpu/multiarch/Makefile
@@ -4,9 +4,9 @@ libm-sysdep_routines += s_floor-c s_ceil
 
 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 slowexp-fma4 slowpow-fma4 \
+			mplog-fma4 mpa-fma4 \
 			sincos32-fma4 doasin-fma4 dosincos-fma4 \
-			halfulp-fma4 mpexp-fma4 \
+			mpexp-fma4 \
 			mpatan2-fma4 mpatan-fma4 mpsqrt-fma4 mptan-fma4
 
 CFLAGS-doasin-fma4.c = -mfma4
@@ -16,7 +16,6 @@ CFLAGS-e_atan2-fma4.c = -mfma4
 CFLAGS-e_exp-fma4.c = -mfma4
 CFLAGS-e_log-fma4.c = -mfma4
 CFLAGS-e_pow-fma4.c = -mfma4 $(config-cflags-nofma)
-CFLAGS-halfulp-fma4.c = -mfma4
 CFLAGS-mpa-fma4.c = -mfma4
 CFLAGS-mpatan-fma4.c = -mfma4
 CFLAGS-mpatan2-fma4.c = -mfma4
@@ -26,14 +25,12 @@ CFLAGS-mpsqrt-fma4.c = -mfma4
 CFLAGS-mptan-fma4.c = -mfma4
 CFLAGS-s_atan-fma4.c = -mfma4
 CFLAGS-sincos32-fma4.c = -mfma4
-CFLAGS-slowexp-fma4.c = -mfma4
-CFLAGS-slowpow-fma4.c = -mfma4
 CFLAGS-s_sin-fma4.c = -mfma4
 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 slowexp-avx \
+			mplog-avx mpa-avx \
 			mpexp-avx
 
 CFLAGS-e_atan2-avx.c = -msse2avx -DSSE2AVX
@@ -44,7 +41,6 @@ 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-slowexp-avx.c = -msse2avx -DSSE2AVX
 CFLAGS-s_tan-avx.c = -msse2avx -DSSE2AVX
 endif
 
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_exp-avx.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/e_exp-avx.c
+++ glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_exp-avx.c
@@ -1,6 +1,5 @@
 #define __ieee754_exp __ieee754_exp_avx
 #define __exp1 __exp1_avx
-#define __slowexp __slowexp_avx
 #define SECTION __attribute__ ((section (".text.avx")))
 
 #include <sysdeps/ieee754/dbl-64/e_exp.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c
+++ glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c
@@ -1,6 +1,5 @@
 #define __ieee754_exp __ieee754_exp_fma4
 #define __exp1 __exp1_fma4
-#define __slowexp __slowexp_fma4
 #define SECTION __attribute__ ((section (".text.fma4")))
 
 #include <sysdeps/ieee754/dbl-64/e_exp.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c
+++ glibc-2.26/sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c
@@ -1,6 +1,5 @@
 #define __ieee754_pow __ieee754_pow_fma4
 #define __exp1 __exp1_fma4
-#define __slowpow __slowpow_fma4
 #define SECTION __attribute__ ((section (".text.fma4")))
 
 #include <sysdeps/ieee754/dbl-64/e_pow.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c
+++ /dev/null
@@ -1,4 +0,0 @@
-#define __halfulp __halfulp_fma4
-#define SECTION __attribute__ ((section (".text.fma4")))
-
-#include <sysdeps/ieee754/dbl-64/halfulp.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/slowexp-avx.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/slowexp-avx.c
+++ /dev/null
@@ -1,9 +0,0 @@
-#define __slowexp __slowexp_avx
-#define __add __add_avx
-#define __dbl_mp __dbl_mp_avx
-#define __mpexp __mpexp_avx
-#define __mul __mul_avx
-#define __sub __sub_avx
-#define SECTION __attribute__ ((section (".text.avx")))
-
-#include <sysdeps/ieee754/dbl-64/slowexp.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c
+++ /dev/null
@@ -1,9 +0,0 @@
-#define __slowexp __slowexp_fma4
-#define __add __add_fma4
-#define __dbl_mp __dbl_mp_fma4
-#define __mpexp __mpexp_fma4
-#define __mul __mul_fma4
-#define __sub __sub_fma4
-#define SECTION __attribute__ ((section (".text.fma4")))
-
-#include <sysdeps/ieee754/dbl-64/slowexp.c>
Index: glibc-2.26/sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c
===================================================================
--- glibc-2.26.orig/sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c
+++ /dev/null
@@ -1,11 +0,0 @@
-#define __slowpow __slowpow_fma4
-#define __add __add_fma4
-#define __dbl_mp __dbl_mp_fma4
-#define __mpexp __mpexp_fma4
-#define __mplog __mplog_fma4
-#define __mul __mul_fma4
-#define __sub __sub_fma4
-#define __halfulp __halfulp_fma4
-#define SECTION __attribute__ ((section (".text.fma4")))
-
-#include <sysdeps/ieee754/dbl-64/slowpow.c>

Places

File math-remove-slow-path.patch of Package glibc (Revision e53058db315273b86de539470702f9df)

Places