Sign Up
Log In
Log In
or
Sign Up
Places
All Projects
Status Monitor
Collapse sidebar
SUSE:SLE-15-SP4:GA
glibc.9959
math-remove-slow-path.patch
Overview
Repositories
Revisions
Requests
Users
Attributes
Meta
File math-remove-slow-path.patch of Package glibc.9959
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>
Locations
Projects
Search
Status Monitor
Help
OpenBuildService.org
Documentation
API Documentation
Code of Conduct
Contact
Support
@OBShq
Terms
openSUSE Build Service is sponsored by
The Open Build Service is an
openSUSE project
.
Sign Up
Log In
Places
Places
All Projects
Status Monitor
