diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64')
157 files changed, 1790 insertions, 2460 deletions
diff --git a/sysdeps/ieee754/dbl-64/Makefile b/sysdeps/ieee754/dbl-64/Makefile index 5557c75b45..c965982fa5 100644 --- a/sysdeps/ieee754/dbl-64/Makefile +++ b/sysdeps/ieee754/dbl-64/Makefile @@ -1,6 +1,6 @@ ifeq ($(subdir),math) # branred depends on precise IEEE double rounding -CFLAGS-branred.c = $(config-cflags-nofma) -CFLAGS-e_sqrt.c = $(config-cflags-nofma) -CFLAGS-e_pow.c = $(config-cflags-nofma) +CFLAGS-branred.c += $(config-cflags-nofma) +CFLAGS-e_sqrt.c += $(config-cflags-nofma) +CFLAGS-e_pow.c += $(config-cflags-nofma) endif diff --git a/sysdeps/ieee754/dbl-64/MathLib.h b/sysdeps/ieee754/dbl-64/MathLib.h index 94182e45c1..671d78f72b 100644 --- a/sysdeps/ieee754/dbl-64/MathLib.h +++ b/sysdeps/ieee754/dbl-64/MathLib.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/asincos.tbl b/sysdeps/ieee754/dbl-64/asincos.tbl index 8da9ff3c07..bb9b7c7c07 100644 --- a/sysdeps/ieee754/dbl-64/asincos.tbl +++ b/sysdeps/ieee754/dbl-64/asincos.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/atnat.h b/sysdeps/ieee754/dbl-64/atnat.h index a9649c57ea..ff8671ac48 100644 --- a/sysdeps/ieee754/dbl-64/atnat.h +++ b/sysdeps/ieee754/dbl-64/atnat.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/atnat2.h b/sysdeps/ieee754/dbl-64/atnat2.h index d5f4ab946d..1c7d3cb760 100644 --- a/sysdeps/ieee754/dbl-64/atnat2.h +++ b/sysdeps/ieee754/dbl-64/atnat2.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/branred.c b/sysdeps/ieee754/dbl-64/branred.c index ca9dd90aff..b1490cce2d 100644 --- a/sysdeps/ieee754/dbl-64/branred.c +++ b/sysdeps/ieee754/dbl-64/branred.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/branred.h b/sysdeps/ieee754/dbl-64/branred.h index 7b6b5473b6..0c0f054e62 100644 --- a/sysdeps/ieee754/dbl-64/branred.h +++ b/sysdeps/ieee754/dbl-64/branred.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dbl2mpn.c b/sysdeps/ieee754/dbl-64/dbl2mpn.c index cb5a070914..8d6e8a1836 100644 --- a/sysdeps/ieee754/dbl-64/dbl2mpn.c +++ b/sysdeps/ieee754/dbl-64/dbl2mpn.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1993-2016 Free Software Foundation, Inc. +/* Copyright (C) 1993-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/dbl-64/dla.h b/sysdeps/ieee754/dbl-64/dla.h index d21c47a68f..5196759ca7 100644 --- a/sysdeps/ieee754/dbl-64/dla.h +++ b/sysdeps/ieee754/dbl-64/dla.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -57,6 +57,10 @@ z=(x)-(y); zz=(fabs(x)>fabs(y)) ? (((x)-(z))-(y)) : ((x)-((y)+(z))); +#ifdef __FP_FAST_FMA +# define DLA_FMS(x, y, z) __builtin_fma (x, y, -(z)) +#endif + /* Exact multiplication of two single-length floating point numbers, */ /* Veltkamp. The macro produces a double-length number (z,zz) that */ /* satisfies z+zz = x*y exactly. p,hx,tx,hy,ty are temporary */ diff --git a/sysdeps/ieee754/dbl-64/doasin.c b/sysdeps/ieee754/dbl-64/doasin.c index 903b5484ba..9355333780 100644 --- a/sysdeps/ieee754/dbl-64/doasin.c +++ b/sysdeps/ieee754/dbl-64/doasin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/doasin.h b/sysdeps/ieee754/dbl-64/doasin.h index 044d597f9c..240c7ec67c 100644 --- a/sysdeps/ieee754/dbl-64/doasin.h +++ b/sysdeps/ieee754/dbl-64/doasin.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dosincos.c b/sysdeps/ieee754/dbl-64/dosincos.c index 931975e6c3..689b453f8d 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.c +++ b/sysdeps/ieee754/dbl-64/dosincos.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/dosincos.h b/sysdeps/ieee754/dbl-64/dosincos.h index 02772592e8..23858029e9 100644 --- a/sysdeps/ieee754/dbl-64/dosincos.h +++ b/sysdeps/ieee754/dbl-64/dosincos.h @@ -2,7 +2,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/e_acosh.c b/sysdeps/ieee754/dbl-64/e_acosh.c index c1f3590f75..fe0c375f00 100644 --- a/sysdeps/ieee754/dbl-64/e_acosh.c +++ b/sysdeps/ieee754/dbl-64/e_acosh.c @@ -36,7 +36,7 @@ __ieee754_acosh (double x) { double t; int32_t hx; - u_int32_t lx; + uint32_t lx; EXTRACT_WORDS (hx, lx, x); if (hx < 0x3ff00000) /* x < 1 */ { @@ -58,12 +58,12 @@ __ieee754_acosh (double x) else if (hx > 0x40000000) /* 2**28 > x > 2 */ { t = x * x; - return __ieee754_log (2.0 * x - one / (x + __ieee754_sqrt (t - one))); + return __ieee754_log (2.0 * x - one / (x + sqrt (t - one))); } else /* 1<x<2 */ { t = x - one; - return __log1p (t + __ieee754_sqrt (2.0 * t + t * t)); + return __log1p (t + sqrt (2.0 * t + t * t)); } } strong_alias (__ieee754_acosh, __acosh_finite) diff --git a/sysdeps/ieee754/dbl-64/e_asin.c b/sysdeps/ieee754/dbl-64/e_asin.c index d8c012d1d3..6bf56945a6 100644 --- a/sysdeps/ieee754/dbl-64/e_asin.c +++ b/sysdeps/ieee754/dbl-64/e_asin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -42,6 +42,7 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> #ifndef SECTION # define SECTION @@ -323,7 +324,7 @@ __ieee754_asin(double x){ /*---------------------------- |x|>=1 -------------------------------*/ else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?hp0.x:-hp0.x; else - if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x; + if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x; else { u.i[HIGH_HALF]=0x7ff00000; v.i[HIGH_HALF]=0x7ff00000; @@ -633,7 +634,7 @@ __ieee754_acos(double x) else if (k==0x3ff00000 && u.i[LOW_HALF]==0) return (m>0)?0:2.0*hp0.x; else - if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x; + if (k>0x7ff00000 || (k == 0x7ff00000 && u.i[LOW_HALF] != 0)) return x + x; else { u.i[HIGH_HALF]=0x7ff00000; v.i[HIGH_HALF]=0x7ff00000; diff --git a/sysdeps/ieee754/dbl-64/e_atan2.c b/sysdeps/ieee754/dbl-64/e_atan2.c index 22e8fb8fef..7295067507 100644 --- a/sysdeps/ieee754/dbl-64/e_atan2.c +++ b/sysdeps/ieee754/dbl-64/e_atan2.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -44,6 +44,7 @@ #include <fenv.h> #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> #include <stap-probe.h> @@ -91,7 +92,7 @@ __ieee754_atan2 (double y, double x) if ((ux & 0x7ff00000) == 0x7ff00000) { if (((ux & 0x000fffff) | dx) != 0x00000000) - return x + x; + return x + y; } num.d = y; uy = num.i[HIGH_HALF]; diff --git a/sysdeps/ieee754/dbl-64/e_atanh.c b/sysdeps/ieee754/dbl-64/e_atanh.c index 395118e92a..da4da8270c 100644 --- a/sysdeps/ieee754/dbl-64/e_atanh.c +++ b/sysdeps/ieee754/dbl-64/e_atanh.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2016 Free Software Foundation, Inc. +/* Copyright (C) 2011-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -38,7 +38,9 @@ #include <float.h> #include <inttypes.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> static const double huge = 1e300; diff --git a/sysdeps/ieee754/dbl-64/e_cosh.c b/sysdeps/ieee754/dbl-64/e_cosh.c index 52a5d5007d..ae2180aa89 100644 --- a/sysdeps/ieee754/dbl-64/e_cosh.c +++ b/sysdeps/ieee754/dbl-64/e_cosh.c @@ -32,6 +32,7 @@ */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> static const double one = 1.0, half = 0.5, huge = 1.0e300; @@ -41,7 +42,7 @@ __ieee754_cosh (double x) { double t, w; int32_t ix; - u_int32_t lx; + uint32_t lx; /* High word of |x|. */ GET_HIGH_WORD (ix, x); @@ -71,7 +72,7 @@ __ieee754_cosh (double x) /* |x| in [log(maxdouble), overflowthresold] */ GET_LOW_WORD (lx, x); - if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (u_int32_t) 0x8fb9f87d))) + if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (uint32_t) 0x8fb9f87d))) { w = __ieee754_exp (half * fabs (x)); t = half * w; diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c index ad1bc84625..ddd2bcb1c2 100644 --- a/sysdeps/ieee754/dbl-64/e_exp.c +++ b/sysdeps/ieee754/dbl-64/e_exp.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -23,10 +23,9 @@ /* exp1 */ /* */ /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */ -/* mpa.c mpexp.x slowexp.c */ /* */ /* 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. */ /* */ @@ -38,25 +37,25 @@ #include "mydefs.h" #include "MathLib.h" #include "uexp.tbl" +#include <math-barriers.h> #include <math_private.h> #include <fenv.h> #include <float.h> +#include "eexp.tbl" #ifndef SECTION # 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) { double bexp, t, eps, del, base, y, al, bet, res, rem, cor; + double z; mynumber junk1, junk2, binexp = {{0, 0}}; int4 i, j, m, n, ex; + int4 k; double retval; { @@ -66,7 +65,42 @@ __ieee754_exp (double x) m = junk1.i[HIGH_HALF]; n = m & hugeint; - if (n > smallint && n < bigint) + if (n < 0x3ff0a2b2) /* |x| < 1.03972053527832 */ + { + if (n < 0x3f862e42) /* |x| < 3/2 ln 2 */ + { + if (n < 0x3ed00000) /* |x| < 1/64 ln 2 */ + { + if (n < 0x3e300000) /* |x| < 2^18 */ + { + retval = one + junk1.x; + goto ret; + } + retval = one + junk1.x * (one + half * junk1.x); + goto ret; + } + t = junk1.x * junk1.x; + retval = junk1.x + (t * (half + junk1.x * t2) + + (t * t) * (t3 + junk1.x * t4 + t * t5)); + retval = one + retval; + goto ret; + } + + /* Find the multiple of 2^-6 nearest x. */ + k = n >> 20; + j = (0x00100000 | (n & 0x000fffff)) >> (0x40c - k); + j = (j - 1) & ~1; + if (m < 0) + j += 134; + z = junk1.x - TBL2[j]; + t = z * z; + retval = z + (t * (half + (z * t2)) + + (t * t) * (t3 + z * t4 + t * t5)); + retval = TBL2[j + 1] + TBL2[j + 1] * retval; + goto ret; + } + + if (n < bigint) /* && |x| >= 1.03972053527832 */ { y = x * log2e.x + three51.x; bexp = y - three51.x; /* multiply the result by 2**bexp */ @@ -93,22 +127,9 @@ __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 */ - } - - if (n <= smallint) - { - retval = 1.0; + /* Maximum relative error is 7.8e-22 (70.1 bits). + Maximum ULP error is 0.500007. */ + retval = res * binexp.x; goto ret; } @@ -166,38 +187,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 +215,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 +237,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 +250,7 @@ __exp1 (double x, double xx, double error) 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 +278,9 @@ __exp1 (double x, double xx, double error) 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 +318,7 @@ __exp1 (double x, double xx, double error) cor = (al - res) + rem; if (m >> 31) { + /* x < 0. */ ex = junk1.i[LOW_HALF]; if (res < 1.0) { @@ -328,34 +329,25 @@ __exp1 (double x, double xx, double error) 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; } } diff --git a/sysdeps/ieee754/dbl-64/e_exp10.c b/sysdeps/ieee754/dbl-64/e_exp10.c index e3b39f7e31..0721143a69 100644 --- a/sysdeps/ieee754/dbl-64/e_exp10.c +++ b/sysdeps/ieee754/dbl-64/e_exp10.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2012-2016 Free Software Foundation, Inc. +/* Copyright (C) 2012-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c index 7402bd7d89..c45bb44744 100644 --- a/sysdeps/ieee754/dbl-64/e_exp2.c +++ b/sysdeps/ieee754/dbl-64/e_exp2.c @@ -1,5 +1,5 @@ /* Double-precision floating point 2^x. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> @@ -29,7 +29,9 @@ #include <math.h> #include <fenv.h> #include <inttypes.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> #include "t_exp2.h" diff --git a/sysdeps/ieee754/dbl-64/e_fmod.c b/sysdeps/ieee754/dbl-64/e_fmod.c index e82b302200..1a8c14dc2a 100644 --- a/sysdeps/ieee754/dbl-64/e_fmod.c +++ b/sysdeps/ieee754/dbl-64/e_fmod.c @@ -24,7 +24,7 @@ double __ieee754_fmod (double x, double y) { int32_t n, hx, hy, hz, ix, iy, sx, i; - u_int32_t lx, ly, lz; + uint32_t lx, ly, lz; EXTRACT_WORDS (hx, lx, x); EXTRACT_WORDS (hy, ly, y); @@ -41,7 +41,7 @@ __ieee754_fmod (double x, double y) if ((hx < hy) || (lx < ly)) return x; /* |x|<|y| return x */ if (lx == ly) - return Zero[(u_int32_t) sx >> 31]; /* |x|=|y| return x*0*/ + return Zero[(uint32_t) sx >> 31]; /* |x|=|y| return x*0*/ } /* determine ix = ilogb(x) */ @@ -125,7 +125,7 @@ __ieee754_fmod (double x, double y) else { if ((hz | lz) == 0) /* return sign(x)*0 */ - return Zero[(u_int32_t) sx >> 31]; + return Zero[(uint32_t) sx >> 31]; hx = hz + hz + (lz >> 31); lx = lz + lz; } } @@ -138,7 +138,7 @@ __ieee754_fmod (double x, double y) /* convert back to floating value and restore the sign */ if ((hx | lx) == 0) /* return sign(x)*0 */ - return Zero[(u_int32_t) sx >> 31]; + return Zero[(uint32_t) sx >> 31]; while (hx < 0x00100000) /* normalize x */ { hx = hx + hx + (lx >> 31); lx = lx + lx; @@ -154,7 +154,7 @@ __ieee754_fmod (double x, double y) n = -1022 - iy; if (n <= 20) { - lx = (lx >> n) | ((u_int32_t) hx << (32 - n)); + lx = (lx >> n) | ((uint32_t) hx << (32 - n)); hx >>= n; } else if (n <= 31) diff --git a/sysdeps/ieee754/dbl-64/e_gamma_r.c b/sysdeps/ieee754/dbl-64/e_gamma_r.c index 2d156850d5..2744549cbd 100644 --- a/sysdeps/ieee754/dbl-64/e_gamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_gamma_r.c @@ -1,5 +1,5 @@ /* Implementation of gamma function according to ISO C. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -18,7 +18,9 @@ <http://www.gnu.org/licenses/>. */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> #include <float.h> /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's @@ -98,7 +100,7 @@ gamma_positive (double x, int *exp2_adj) double ret = (__ieee754_pow (x_adj_mant, x_adj) * __ieee754_exp2 (x_adj_log2 * x_adj_frac) * __ieee754_exp (-x_adj) - * __ieee754_sqrt (2 * M_PI / x_adj) + * sqrt (2 * M_PI / x_adj) / prod); exp_adj += x_eps * __ieee754_log (x_adj); double bsum = gamma_coeff[NCOEFF - 1]; @@ -114,7 +116,7 @@ double __ieee754_gamma_r (double x, int *signgamp) { int32_t hx; - u_int32_t lx; + uint32_t lx; double ret; EXTRACT_WORDS (hx, lx, x); @@ -126,7 +128,7 @@ __ieee754_gamma_r (double x, int *signgamp) return 1.0 / x; } if (__builtin_expect (hx < 0, 0) - && (u_int32_t) hx < 0xfff00000 && __rint (x) == x) + && (uint32_t) hx < 0xfff00000 && __rint (x) == x) { /* Return value for integer x < 0 is NaN with invalid exception. */ *signgamp = 0; diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c index f142c450a2..a2c33cc4ed 100644 --- a/sysdeps/ieee754/dbl-64/e_hypot.c +++ b/sysdeps/ieee754/dbl-64/e_hypot.c @@ -44,6 +44,7 @@ #include <math.h> #include <math_private.h> +#include <math-underflow.h> double __ieee754_hypot (double x, double y) @@ -74,8 +75,10 @@ __ieee754_hypot (double x, double y) { if (ha >= 0x7ff00000) /* Inf or NaN */ { - u_int32_t low; + uint32_t low; w = a + b; /* for sNaN */ + if (issignaling (a) || issignaling (b)) + return w; GET_LOW_WORD (low, a); if (((ha & 0xfffff) | low) == 0) w = a; @@ -93,7 +96,7 @@ __ieee754_hypot (double x, double y) { if (hb <= 0x000fffff) /* subnormal b or 0 */ { - u_int32_t low; + uint32_t low; GET_LOW_WORD (low, b); if ((hb | low) == 0) return a; @@ -130,7 +133,7 @@ __ieee754_hypot (double x, double y) t1 = 0; SET_HIGH_WORD (t1, ha); t2 = a - t1; - w = __ieee754_sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1))); + w = sqrt (t1 * t1 - (b * (-b) - t2 * (a + t1))); } else { @@ -141,11 +144,11 @@ __ieee754_hypot (double x, double y) t1 = 0; SET_HIGH_WORD (t1, ha + 0x00100000); t2 = a - t1; - w = __ieee754_sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b))); + w = sqrt (t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b))); } if (k != 0) { - u_int32_t high; + uint32_t high; t1 = 1.0; GET_HIGH_WORD (high, t1); SET_HIGH_WORD (t1, high + (k << 20)); diff --git a/sysdeps/ieee754/dbl-64/e_j0.c b/sysdeps/ieee754/dbl-64/e_j0.c index 9f25aa855e..7f5919910d 100644 --- a/sysdeps/ieee754/dbl-64/e_j0.c +++ b/sysdeps/ieee754/dbl-64/e_j0.c @@ -59,6 +59,7 @@ */ #include <math.h> +#include <math-barriers.h> #include <math_private.h> static double pzero (double), qzero (double); @@ -109,11 +110,11 @@ __ieee754_j0 (double x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ if (ix > 0x48000000) - z = (invsqrtpi * cc) / __ieee754_sqrt (x); + z = (invsqrtpi * cc) / sqrt (x); else { u = pzero (x); v = qzero (x); - z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (x); + z = invsqrtpi * (u * cc - v * ss) / sqrt (x); } return z; } @@ -169,7 +170,7 @@ __ieee754_y0 (double x) if (ix >= 0x7ff00000) return one / (x + x * x); if ((ix | lx) == 0) - return -HUGE_VAL + x; /* -inf and overflow exception. */ + return -1 / zero; /* -inf and divide by zero exception. */ if (hx < 0) return zero / (zero * x); if (ix >= 0x40000000) /* |x| >= 2.0 */ @@ -200,11 +201,11 @@ __ieee754_y0 (double x) ss = z / cc; } if (ix > 0x48000000) - z = (invsqrtpi * ss) / __ieee754_sqrt (x); + z = (invsqrtpi * ss) / sqrt (x); else { u = pzero (x); v = qzero (x); - z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x); + z = invsqrtpi * (u * ss + v * cc) / sqrt (x); } return z; } diff --git a/sysdeps/ieee754/dbl-64/e_j1.c b/sysdeps/ieee754/dbl-64/e_j1.c index 4827fbf3d3..734f3ca64a 100644 --- a/sysdeps/ieee754/dbl-64/e_j1.c +++ b/sysdeps/ieee754/dbl-64/e_j1.c @@ -61,7 +61,9 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static double pone (double), qone (double); @@ -112,11 +114,11 @@ __ieee754_j1 (double x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ if (ix > 0x48000000) - z = (invsqrtpi * cc) / __ieee754_sqrt (y); + z = (invsqrtpi * cc) / sqrt (y); else { u = pone (y); v = qone (y); - z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrt (y); + z = invsqrtpi * (u * cc - v * ss) / sqrt (y); } if (hx < 0) return -z; @@ -174,7 +176,7 @@ __ieee754_y1 (double x) if (__glibc_unlikely (ix >= 0x7ff00000)) return one / (x + x * x); if (__glibc_unlikely ((ix | lx) == 0)) - return -HUGE_VAL + x; + return -1 / zero; /* -inf and divide by zero exception. */ /* -inf and overflow exception. */; if (__glibc_unlikely (hx < 0)) return zero / (zero * x); @@ -203,11 +205,11 @@ __ieee754_y1 (double x) * to compute the worse one. */ if (ix > 0x48000000) - z = (invsqrtpi * ss) / __ieee754_sqrt (x); + z = (invsqrtpi * ss) / sqrt (x); else { u = pone (x); v = qone (x); - z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrt (x); + z = invsqrtpi * (u * ss + v * cc) / sqrt (x); } return z; } diff --git a/sysdeps/ieee754/dbl-64/e_jn.c b/sysdeps/ieee754/dbl-64/e_jn.c index 3fecf82f10..9181b22bb8 100644 --- a/sysdeps/ieee754/dbl-64/e_jn.c +++ b/sysdeps/ieee754/dbl-64/e_jn.c @@ -39,7 +39,9 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static const double invsqrtpi = 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ @@ -61,7 +63,7 @@ __ieee754_jn (int n, double x) EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; /* if J(n,NaN) is NaN */ - if (__glibc_unlikely ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000)) + if (__glibc_unlikely ((ix | ((uint32_t) (lx | -lx)) >> 31) > 0x7ff00000)) return x + x; if (n < 0) { @@ -107,7 +109,7 @@ __ieee754_jn (int n, double x) case 2: temp = -c - s; break; case 3: temp = c - s; break; } - b = invsqrtpi * temp / __ieee754_sqrt (x); + b = invsqrtpi * temp / sqrt (x); } else { @@ -266,13 +268,8 @@ __ieee754_yn (int n, double x) EXTRACT_WORDS (hx, lx, x); ix = 0x7fffffff & hx; /* if Y(n,NaN) is NaN */ - if (__glibc_unlikely ((ix | ((u_int32_t) (lx | -lx)) >> 31) > 0x7ff00000)) + if (__glibc_unlikely ((ix | ((uint32_t) (lx | -lx)) >> 31) > 0x7ff00000)) return x + x; - if (__glibc_unlikely ((ix | lx) == 0)) - return -HUGE_VAL + x; - /* -inf and overflow exception. */; - if (__glibc_unlikely (hx < 0)) - return zero / (zero * x); sign = 1; if (n < 0) { @@ -281,6 +278,11 @@ __ieee754_yn (int n, double x) } if (n == 0) return (__ieee754_y0 (x)); + if (__glibc_unlikely ((ix | lx) == 0)) + return -sign / zero; + /* -inf and overflow exception. */; + if (__glibc_unlikely (hx < 0)) + return zero / (zero * x); { SET_RESTORE_ROUND (FE_TONEAREST); if (n == 1) @@ -314,11 +316,11 @@ __ieee754_yn (int n, double x) case 2: temp = -s + c; break; case 3: temp = s + c; break; } - b = invsqrtpi * temp / __ieee754_sqrt (x); + b = invsqrtpi * temp / sqrt (x); } else { - u_int32_t high; + uint32_t high; a = __ieee754_y0 (x); b = __ieee754_y1 (x); /* quit if b is -inf */ diff --git a/sysdeps/ieee754/dbl-64/e_lgamma_r.c b/sysdeps/ieee754/dbl-64/e_lgamma_r.c index 15154c0f43..17717d915f 100644 --- a/sysdeps/ieee754/dbl-64/e_lgamma_r.c +++ b/sysdeps/ieee754/dbl-64/e_lgamma_r.c @@ -77,9 +77,10 @@ * */ -#include <libc-internal.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libc-diag.h> static const double two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ @@ -225,7 +226,7 @@ __ieee754_lgamma_r(double x, int *signgamp) if(hx<0) { if(__builtin_expect(ix>=0x43300000, 0)) /* |x|>=2**52, must be -integer */ - return x/zero; + return fabs (x)/zero; if (x < -2.0 && x > -28.0) return __lgamma_neg (x, signgamp); t = sin_pi(x); diff --git a/sysdeps/ieee754/dbl-64/e_log.c b/sysdeps/ieee754/dbl-64/e_log.c index 9917dc236f..2483dd8551 100644 --- a/sysdeps/ieee754/dbl-64/e_log.c +++ b/sysdeps/ieee754/dbl-64/e_log.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -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); + /* 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; - /* End stage II, case abs(x-1) < 0.03 */ - if ((y = b + (bb + b * E4)) == b + (bb - b * E4)) - return y; - goto stage_n; - - /*--- 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 diff --git a/sysdeps/ieee754/dbl-64/e_log10.c b/sysdeps/ieee754/dbl-64/e_log10.c index df59d9dce4..677cbc4df8 100644 --- a/sysdeps/ieee754/dbl-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/e_log10.c @@ -57,7 +57,7 @@ __ieee754_log10 (double x) { double y, z; int32_t i, k, hx; - u_int32_t lx; + uint32_t lx; EXTRACT_WORDS (hx, lx, x); @@ -65,7 +65,7 @@ __ieee754_log10 (double x) if (hx < 0x00100000) { /* x < 2**-1022 */ if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0)) - return -two54 / (x - x); /* log(+-0)=-inf */ + return -two54 / fabs (x); /* log(+-0)=-inf */ if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; @@ -75,7 +75,7 @@ __ieee754_log10 (double x) if (__glibc_unlikely (hx >= 0x7ff00000)) return x + x; k += (hx >> 20) - 1023; - i = ((u_int32_t) k & 0x80000000) >> 31; + i = ((uint32_t) k & 0x80000000) >> 31; hx = (hx & 0x000fffff) | ((0x3ff - i) << 20); y = (double) (k + i); if (FIX_INT_FP_CONVERT_ZERO && y == 0.0) diff --git a/sysdeps/ieee754/dbl-64/e_log2.c b/sysdeps/ieee754/dbl-64/e_log2.c index bc6a34192a..e4a6aff9a3 100644 --- a/sysdeps/ieee754/dbl-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/e_log2.c @@ -75,7 +75,7 @@ __ieee754_log2 (double x) { double hfsq, f, s, z, R, w, t1, t2, dk; int32_t k, hx, i, j; - u_int32_t lx; + uint32_t lx; EXTRACT_WORDS (hx, lx, x); @@ -83,7 +83,7 @@ __ieee754_log2 (double x) if (hx < 0x00100000) { /* x < 2**-1022 */ if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0)) - return -two54 / (x - x); /* log(+-0)=-inf */ + return -two54 / fabs (x); /* log(+-0)=-inf */ if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 663fa392c2..96d5b23ccc 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -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. */ @@ -42,6 +38,7 @@ #include "MathLib.h" #include "upow.tbl" #include <math_private.h> +#include <math-underflow.h> #include <fenv.h> #ifndef SECTION @@ -50,11 +47,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 +57,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; @@ -73,8 +67,9 @@ __ieee754_pow (double x, double y) { /* of y */ qx = u.i[HIGH_HALF] & 0x7fffffff; /* Is x a NaN? */ - if (((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000)) - return x; + if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000)) + && (y != 0 || issignaling (x))) + return x + x; if (y == 1.0) return x; if (y == 2.0) @@ -99,7 +94,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; @@ -110,9 +105,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)) @@ -128,7 +130,7 @@ __ieee754_pow (double x, double y) { if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0) || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000) /* NaN */ - return y; + return y + y; if (fabs (y) > 1.0e20) return (y > 0) ? 0 : 1.0 / 0.0; k = checkint (y); @@ -142,9 +144,9 @@ __ieee754_pow (double x, double y) qy = v.i[HIGH_HALF] & 0x7fffffff; /* no sign */ if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0)) /* NaN */ - return x; + return x + y; if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0)) /* NaN */ - return x == 1.0 ? 1.0 : y; + return x == 1.0 && !issignaling (y) ? 1.0 : y + y; /* if x<0 */ if (u.i[HIGH_HALF] < 0) @@ -217,33 +219,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; @@ -259,9 +239,7 @@ log1 (double x, double *delta, double *error) 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; @@ -283,7 +261,7 @@ log1 (double x, double *delta, double *error) 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 */ { @@ -295,8 +273,8 @@ log1 (double x, double *delta, double *error) * (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 @@ -315,12 +293,12 @@ log1 (double x, double *delta, double *error) 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; @@ -333,113 +311,13 @@ log1 (double x, double *delta, double *error) 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 @@ -451,7 +329,8 @@ checkint (double x) int4 i[2]; double x; } u; - int k, m, n; + int k; + unsigned int m, n; u.x = x; m = u.i[HIGH_HALF] & 0x7fffffff; /* no sign */ if (m >= 0x7ff00000) @@ -466,15 +345,15 @@ checkint (double x) return (n & 1) ? -1 : 1; /* odd or even */ if (k > 20) { - if (n << (k - 20)) + if (n << (k - 20) != 0) return 0; /* if not integer */ - return (n << (k - 21)) ? -1 : 1; + return (n << (k - 21) != 0) ? -1 : 1; } if (n) return 0; /*if not integer */ if (k == 20) return (m & 1) ? -1 : 1; - if (m << (k + 12)) + if (m << (k + 12) != 0) return 0; - return (m << (k + 11)) ? -1 : 1; + return (m << (k + 11) != 0) ? -1 : 1; } diff --git a/sysdeps/ieee754/dbl-64/e_rem_pio2.c b/sysdeps/ieee754/dbl-64/e_rem_pio2.c index 2f55ca294b..81a3d073d4 100644 --- a/sysdeps/ieee754/dbl-64/e_rem_pio2.c +++ b/sysdeps/ieee754/dbl-64/e_rem_pio2.c @@ -74,7 +74,7 @@ __ieee754_rem_pio2 (double x, double *y) double z, w, t, r, fn; double tx[3]; int32_t e0, i, j, nx, n, ix, hx; - u_int32_t low; + uint32_t low; GET_HIGH_WORD (hx, x); /* high word of x */ ix = hx & 0x7fffffff; @@ -130,7 +130,7 @@ __ieee754_rem_pio2 (double x, double *y) } else { - u_int32_t high; + uint32_t high; j = ix >> 20; y[0] = r - w; GET_HIGH_WORD (high, y[0]); diff --git a/sysdeps/ieee754/dbl-64/e_remainder.c b/sysdeps/ieee754/dbl-64/e_remainder.c index a445e74b5c..2e7f0ac1f4 100644 --- a/sysdeps/ieee754/dbl-64/e_remainder.c +++ b/sysdeps/ieee754/dbl-64/e_remainder.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/e_sinh.c b/sysdeps/ieee754/dbl-64/e_sinh.c index 8479bdd9b8..c4e34211ac 100644 --- a/sysdeps/ieee754/dbl-64/e_sinh.c +++ b/sysdeps/ieee754/dbl-64/e_sinh.c @@ -34,7 +34,9 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> static const double one = 1.0, shuge = 1.0e307; @@ -43,7 +45,7 @@ __ieee754_sinh (double x) { double t, w, h; int32_t ix, jx; - u_int32_t lx; + uint32_t lx; /* High word of |x|. */ GET_HIGH_WORD (jx, x); @@ -77,7 +79,7 @@ __ieee754_sinh (double x) /* |x| in [log(maxdouble), overflowthresold] */ GET_LOW_WORD (lx, x); - if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (u_int32_t) 0x8fb9f87d))) + if (ix < 0x408633ce || ((ix == 0x408633ce) && (lx <= (uint32_t) 0x8fb9f87d))) { w = __ieee754_exp (0.5 * fabs (x)); t = h * w; diff --git a/sysdeps/ieee754/dbl-64/e_sqrt.c b/sysdeps/ieee754/dbl-64/e_sqrt.c index 8304a2bb63..f70623e17b 100644 --- a/sysdeps/ieee754/dbl-64/e_sqrt.c +++ b/sysdeps/ieee754/dbl-64/e_sqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -37,6 +37,7 @@ #include <dla.h> #include "MathLib.h" #include "root.tbl" +#include <math-barriers.h> #include <math_private.h> /*********************************************************************/ diff --git a/sysdeps/ieee754/dbl-64/eexp.tbl b/sysdeps/ieee754/dbl-64/eexp.tbl new file mode 100644 index 0000000000..4ee6040638 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/eexp.tbl @@ -0,0 +1,172 @@ +/* EXP function tables - for use in computing double precision exponential + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +/* For i = 0, ..., 66, + TBL2[2*i] is a double precision number near (i+1)*2^-6, and + TBL2[2*i+1] = exp(TBL2[2*i]) to within a relative error less + than 2^-60. + + For i = 67, ..., 133, + TBL2[2*i] is a double precision number near -(i+1)*2^-6, and + TBL2[2*i+1] = exp(TBL2[2*i]) to within a relative error less + than 2^-60. */ + +static const double TBL2[268] = { + 0x1.ffffffffffc82p-7, 0x1.04080ab55de32p+0, + 0x1.fffffffffffdbp-6, 0x1.08205601127ecp+0, + 0x1.80000000000a0p-5, 0x1.0c49236829e91p+0, + 0x1.fffffffffff79p-5, 0x1.1082b577d34e9p+0, + 0x1.3fffffffffffcp-4, 0x1.14cd4fc989cd6p+0, + 0x1.8000000000060p-4, 0x1.192937074e0d4p+0, + 0x1.c000000000061p-4, 0x1.1d96b0eff0e80p+0, + 0x1.fffffffffffd6p-4, 0x1.2216045b6f5cap+0, + 0x1.1ffffffffff58p-3, 0x1.26a7793f6014cp+0, + 0x1.3ffffffffff75p-3, 0x1.2b4b58b372c65p+0, + 0x1.5ffffffffff00p-3, 0x1.3001ecf601ad1p+0, + 0x1.8000000000020p-3, 0x1.34cb8170b583ap+0, + 0x1.9ffffffffa629p-3, 0x1.39a862bd3b344p+0, + 0x1.c00000000000fp-3, 0x1.3e98deaa11dcep+0, + 0x1.e00000000007fp-3, 0x1.439d443f5f16dp+0, + 0x1.0000000000072p-2, 0x1.48b5e3c3e81abp+0, + 0x1.0fffffffffecap-2, 0x1.4de30ec211dfbp+0, + 0x1.1ffffffffff8fp-2, 0x1.5325180cfacd2p+0, + 0x1.300000000003bp-2, 0x1.587c53c5a7b04p+0, + 0x1.4000000000034p-2, 0x1.5de9176046007p+0, + 0x1.4ffffffffff89p-2, 0x1.636bb9a98322fp+0, + 0x1.5ffffffffffe7p-2, 0x1.690492cbf942ap+0, + 0x1.6ffffffffff78p-2, 0x1.6eb3fc55b1e45p+0, + 0x1.7ffffffffff65p-2, 0x1.747a513dbef32p+0, + 0x1.8ffffffffffd5p-2, 0x1.7a57ede9ea22ep+0, + 0x1.9ffffffffff6ep-2, 0x1.804d30347b50fp+0, + 0x1.affffffffffc3p-2, 0x1.865a7772164aep+0, + 0x1.c000000000053p-2, 0x1.8c802477b0030p+0, + 0x1.d00000000004dp-2, 0x1.92be99a09bf1ep+0, + 0x1.e000000000096p-2, 0x1.99163ad4b1e08p+0, + 0x1.efffffffffefap-2, 0x1.9f876d8e8c4fcp+0, + 0x1.fffffffffffd0p-2, 0x1.a61298e1e0688p+0, + 0x1.0800000000002p-1, 0x1.acb82581eee56p+0, + 0x1.100000000001fp-1, 0x1.b3787dc80f979p+0, + 0x1.17ffffffffff8p-1, 0x1.ba540dba56e4fp+0, + 0x1.1fffffffffffap-1, 0x1.c14b431256441p+0, + 0x1.27fffffffffc4p-1, 0x1.c85e8d43f7c9bp+0, + 0x1.2fffffffffffdp-1, 0x1.cf8e5d84758a6p+0, + 0x1.380000000001fp-1, 0x1.d6db26d16cd84p+0, + 0x1.3ffffffffffd8p-1, 0x1.de455df80e39bp+0, + 0x1.4800000000052p-1, 0x1.e5cd799c6a59cp+0, + 0x1.4ffffffffffc8p-1, 0x1.ed73f240dc10cp+0, + 0x1.5800000000013p-1, 0x1.f539424d90f71p+0, + 0x1.5ffffffffffbcp-1, 0x1.fd1de6182f885p+0, + 0x1.680000000002dp-1, 0x1.02912df5ce741p+1, + 0x1.7000000000040p-1, 0x1.06a39207f0a2ap+1, + 0x1.780000000004fp-1, 0x1.0ac660691652ap+1, + 0x1.7ffffffffff6fp-1, 0x1.0ef9db467dcabp+1, + 0x1.87fffffffffe5p-1, 0x1.133e45d82e943p+1, + 0x1.9000000000035p-1, 0x1.1793e4652cc6dp+1, + 0x1.97fffffffffb3p-1, 0x1.1bfafc47bda48p+1, + 0x1.a000000000000p-1, 0x1.2073d3f1bd518p+1, + 0x1.a80000000004ap-1, 0x1.24feb2f105ce2p+1, + 0x1.affffffffffedp-1, 0x1.299be1f3e7f11p+1, + 0x1.b7ffffffffffbp-1, 0x1.2e4baacdb6611p+1, + 0x1.c00000000001dp-1, 0x1.330e587b62b39p+1, + 0x1.c800000000079p-1, 0x1.37e437282d538p+1, + 0x1.cffffffffff51p-1, 0x1.3ccd943268248p+1, + 0x1.d7fffffffff74p-1, 0x1.41cabe304cadcp+1, + 0x1.e000000000011p-1, 0x1.46dc04f4e5343p+1, + 0x1.e80000000001ep-1, 0x1.4c01b9950a124p+1, + 0x1.effffffffff9ep-1, 0x1.513c2e6c73196p+1, + 0x1.f7fffffffffedp-1, 0x1.568bb722dd586p+1, + 0x1.0000000000034p+0, 0x1.5bf0a8b1457b0p+1, + 0x1.03fffffffffe2p+0, 0x1.616b5967376dfp+1, + 0x1.07fffffffff4bp+0, 0x1.66fc20f0337a9p+1, + 0x1.0bffffffffffdp+0, 0x1.6ca35859290f5p+1, + -0x1.fffffffffffe4p-7, 0x1.f80feabfeefa5p-1, + -0x1.ffffffffffb0bp-6, 0x1.f03f56a88b5fep-1, + -0x1.7ffffffffffa7p-5, 0x1.e88dc6afecfc5p-1, + -0x1.ffffffffffea8p-5, 0x1.e0fabfbc702b8p-1, + -0x1.3ffffffffffb3p-4, 0x1.d985c89d041acp-1, + -0x1.7ffffffffffe3p-4, 0x1.d22e6a0197c06p-1, + -0x1.bffffffffff9ap-4, 0x1.caf42e73a4c89p-1, + -0x1.fffffffffff98p-4, 0x1.c3d6a24ed822dp-1, + -0x1.1ffffffffffe9p-3, 0x1.bcd553b9d7b67p-1, + -0x1.3ffffffffffe0p-3, 0x1.b5efd29f24c2dp-1, + -0x1.5fffffffff553p-3, 0x1.af25b0a61a9f4p-1, + -0x1.7ffffffffff8bp-3, 0x1.a876812c08794p-1, + -0x1.9fffffffffe51p-3, 0x1.a1e1d93d68828p-1, + -0x1.bffffffffff6ep-3, 0x1.9b674f8f2f3f5p-1, + -0x1.dffffffffff7fp-3, 0x1.95067c7837a0cp-1, + -0x1.fffffffffff7ap-3, 0x1.8ebef9eac8225p-1, + -0x1.0fffffffffffep-2, 0x1.8890636e31f55p-1, + -0x1.1ffffffffff41p-2, 0x1.827a56188975ep-1, + -0x1.2ffffffffffbap-2, 0x1.7c7c708877656p-1, + -0x1.3fffffffffff8p-2, 0x1.769652df22f81p-1, + -0x1.4ffffffffff90p-2, 0x1.70c79eba33c2fp-1, + -0x1.5ffffffffffdbp-2, 0x1.6b0ff72deb8aap-1, + -0x1.6ffffffffff9ap-2, 0x1.656f00bf5798ep-1, + -0x1.7ffffffffff9fp-2, 0x1.5fe4615e98eb0p-1, + -0x1.8ffffffffffeep-2, 0x1.5a6fc061433cep-1, + -0x1.9fffffffffc4ap-2, 0x1.5510c67cd26cdp-1, + -0x1.affffffffff30p-2, 0x1.4fc71dc13566bp-1, + -0x1.bfffffffffff0p-2, 0x1.4a9271936fd0ep-1, + -0x1.cfffffffffff3p-2, 0x1.45726ea84fb8cp-1, + -0x1.dfffffffffff3p-2, 0x1.4066c2ff3912bp-1, + -0x1.effffffffff80p-2, 0x1.3b6f1ddd05ab9p-1, + -0x1.fffffffffffdfp-2, 0x1.368b2fc6f9614p-1, + -0x1.0800000000000p-1, 0x1.31baaa7dca843p-1, + -0x1.0ffffffffffa4p-1, 0x1.2cfd40f8bdce4p-1, + -0x1.17fffffffff0ap-1, 0x1.2852a760d5ce7p-1, + -0x1.2000000000000p-1, 0x1.23ba930c1568bp-1, + -0x1.27fffffffffbbp-1, 0x1.1f34ba78d568dp-1, + -0x1.2fffffffffe32p-1, 0x1.1ac0d5492c1dbp-1, + -0x1.37ffffffff042p-1, 0x1.165e9c3e67ef2p-1, + -0x1.3ffffffffff77p-1, 0x1.120dc93499431p-1, + -0x1.47fffffffff6bp-1, 0x1.0dce171e34ecep-1, + -0x1.4fffffffffff1p-1, 0x1.099f41ffbe588p-1, + -0x1.57ffffffffe02p-1, 0x1.058106eb8a7aep-1, + -0x1.5ffffffffffe5p-1, 0x1.017323fd9002ep-1, + -0x1.67fffffffffb0p-1, 0x1.faeab0ae9386cp-2, + -0x1.6ffffffffffb2p-1, 0x1.f30ec837503d7p-2, + -0x1.77fffffffff7fp-1, 0x1.eb5210d627133p-2, + -0x1.7ffffffffffe8p-1, 0x1.e3b40ebefcd95p-2, + -0x1.87fffffffffc8p-1, 0x1.dc3448110dae2p-2, + -0x1.8fffffffffb30p-1, 0x1.d4d244cf4ef06p-2, + -0x1.97fffffffffefp-1, 0x1.cd8d8ed8ee395p-2, + -0x1.9ffffffffffa7p-1, 0x1.c665b1e1f1e5cp-2, + -0x1.a7fffffffffdcp-1, 0x1.bf5a3b6bf18d6p-2, + -0x1.affffffffff95p-1, 0x1.b86ababeef93bp-2, + -0x1.b7fffffffffcbp-1, 0x1.b196c0e24d256p-2, + -0x1.bffffffffff32p-1, 0x1.aadde095dadf7p-2, + -0x1.c7fffffffff6ap-1, 0x1.a43fae4b047c9p-2, + -0x1.cffffffffffb6p-1, 0x1.9dbbc01e182a4p-2, + -0x1.d7fffffffffcap-1, 0x1.9751adcfa81ecp-2, + -0x1.dffffffffffcdp-1, 0x1.910110be0699ep-2, + -0x1.e7ffffffffffbp-1, 0x1.8ac983dedbc69p-2, + -0x1.effffffffff88p-1, 0x1.84aaa3b8d51a9p-2, + -0x1.f7fffffffffbbp-1, 0x1.7ea40e5d6d92ep-2, + -0x1.fffffffffffdbp-1, 0x1.78b56362cef53p-2, + -0x1.03fffffffff00p+0, 0x1.72de43ddcb1f2p-2, + -0x1.07ffffffffe6fp+0, 0x1.6d1e525bed085p-2, + -0x1.0bfffffffffd6p+0, 0x1.677532dda1c57p-2}; + +static const double + half = 0.5, + one = 1.0, +/* t2-t5 terms used for polynomial computation. */ + t2 = 0x1.5555555555555p-3, /* 1.6666666666666665741e-1 */ + t3 = 0x1.5555555555555p-5, /* 4.1666666666666664354e-2 */ + t4 = 0x1.1111111111111p-7, /* 8.3333333333333332177e-3 */ + t5 = 0x1.6c16c16c16c17p-10; /* 1.3888888888888889419e-3 */ diff --git a/sysdeps/ieee754/dbl-64/gamma_product.c b/sysdeps/ieee754/dbl-64/gamma_product.c index 7ae144aeb3..3300b5139c 100644 --- a/sysdeps/ieee754/dbl-64/gamma_product.c +++ b/sysdeps/ieee754/dbl-64/gamma_product.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,37 +18,7 @@ #include <math.h> #include <math_private.h> -#include <float.h> - -/* Calculate X * Y exactly and store the result in *HI + *LO. It is - given that the values are small enough that no overflow occurs and - large enough (or zero) that no underflow occurs. */ - -static void -mul_split (double *hi, double *lo, double x, double y) -{ -#ifdef __FP_FAST_FMA - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fma (x, y, -*hi); -#elif defined FP_FAST_FMA - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fma (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) - double x1 = x * C; - double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - double x2 = x - x1; - double y2 = y - y1; - *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; -#endif -} +#include <mul_split.h> /* Compute the product of X + X_EPS, X + X_EPS + 1, ..., X + X_EPS + N - 1, in the form R * (1 + *EPS) where the return value R is an diff --git a/sysdeps/ieee754/dbl-64/gamma_productf.c b/sysdeps/ieee754/dbl-64/gamma_productf.c index 58b4e761fe..011c6ff06e 100644 --- a/sysdeps/ieee754/dbl-64/gamma_productf.c +++ b/sysdeps/ieee754/dbl-64/gamma_productf.c @@ -1,5 +1,5 @@ /* Compute a product of X, X+1, ..., with an error estimate. - Copyright (C) 2013-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -17,6 +17,7 @@ <http://www.gnu.org/licenses/>. */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> #include <float.h> diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c deleted file mode 100644 index 5e3e731754..0000000000 --- a/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-2016 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; -} diff --git a/sysdeps/ieee754/dbl-64/k_rem_pio2.c b/sysdeps/ieee754/dbl-64/k_rem_pio2.c index e58c9e854c..d8403dc345 100644 --- a/sysdeps/ieee754/dbl-64/k_rem_pio2.c +++ b/sysdeps/ieee754/dbl-64/k_rem_pio2.c @@ -131,7 +131,9 @@ static char rcsid[] = "$NetBSD: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $ */ #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libc-diag.h> static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ @@ -251,8 +253,17 @@ recompute: j |= iq[i]; if (j == 0) /* need recomputation */ { + /* On s390x gcc 6.1 -O3 produces the warning "array subscript is below + array bounds [-Werror=array-bounds]". Only __ieee754_rem_pio2l + calls __kernel_rem_pio2 for normal numbers and |x| > pi/4 in case + of ldbl-96 and |x| > 3pi/4 in case of ldbl-128[ibm]. + Thus x can't be zero and ipio2 is not zero, too. Thus not all iq[] + values can't be zero. */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_NEEDS_COMMENT (6.1, "-Warray-bounds"); for (k = 1; iq[jk - k] == 0; k++) ; /* k = no. of terms needed */ + DIAG_POP_NEEDS_COMMENT; for (i = jz + 1; i <= jz + k; i++) /* add q[jz+1] to q[jz+k] */ { @@ -319,7 +330,16 @@ recompute: for (i = jz; i >= 0; i--) fv = math_narrow_eval (fv + fq[i]); y[0] = (ih == 0) ? fv : -fv; + /* GCC mainline (to be GCC 9), as of 2018-05-22 on i686, warns + that fq[0] may be used uninitialized. This is not possible + because jz is always nonnegative when the above loop + initializing fq is executed, because the result is never zero + to full precision (this function is not called for zero + arguments). */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_NEEDS_COMMENT (9, "-Wmaybe-uninitialized"); fv = math_narrow_eval (fq[0] - fv); + DIAG_POP_NEEDS_COMMENT; for (i = 1; i <= jz; i++) fv = math_narrow_eval (fv + fq[i]); y[1] = (ih == 0) ? fv : -fv; diff --git a/sysdeps/ieee754/dbl-64/lgamma_neg.c b/sysdeps/ieee754/dbl-64/lgamma_neg.c index fbf992203c..5bb2f10c71 100644 --- a/sysdeps/ieee754/dbl-64/lgamma_neg.c +++ b/sysdeps/ieee754/dbl-64/lgamma_neg.c @@ -1,5 +1,5 @@ /* lgamma expanding around zeros. - Copyright (C) 2015-2016 Free Software Foundation, Inc. + Copyright (C) 2015-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,6 +18,7 @@ #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> static const double lgamma_zeros[][2] = diff --git a/sysdeps/ieee754/dbl-64/lgamma_product.c b/sysdeps/ieee754/dbl-64/lgamma_product.c index d956575bc7..b37ecee73f 100644 --- a/sysdeps/ieee754/dbl-64/lgamma_product.c +++ b/sysdeps/ieee754/dbl-64/lgamma_product.c @@ -1,5 +1,5 @@ /* Compute a product of 1 + (T/X), 1 + (T/(X+1)), .... - Copyright (C) 2015-2016 Free Software Foundation, Inc. + Copyright (C) 2015-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,37 +18,7 @@ #include <math.h> #include <math_private.h> -#include <float.h> - -/* Calculate X * Y exactly and store the result in *HI + *LO. It is - given that the values are small enough that no overflow occurs and - large enough (or zero) that no underflow occurs. */ - -static void -mul_split (double *hi, double *lo, double x, double y) -{ -#ifdef __FP_FAST_FMA - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fma (x, y, -*hi); -#elif defined FP_FAST_FMA - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fma (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) - double x1 = x * C; - double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - double x2 = x - x1; - double y2 = y - y1; - *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; -#endif -} +#include <mul_split.h> /* Compute the product of 1 + (T / (X + X_EPS)), 1 + (T / (X + X_EPS + 1)), ..., 1 + (T / (X + X_EPS + N - 1)), minus 1. X is such that diff --git a/sysdeps/ieee754/dbl-64/mpa-arch.h b/sysdeps/ieee754/dbl-64/mpa-arch.h index ea33a14219..d1a2bf8ac4 100644 --- a/sysdeps/ieee754/dbl-64/mpa-arch.h +++ b/sysdeps/ieee754/dbl-64/mpa-arch.h @@ -1,5 +1,5 @@ /* Overridable constants and operations. - Copyright (C) 2013-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 Free Software Foundation, Inc. This program is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpa.c b/sysdeps/ieee754/dbl-64/mpa.c index 4b21d1d2e6..f59ec6a9d8 100644 --- a/sysdeps/ieee754/dbl-64/mpa.c +++ b/sysdeps/ieee754/dbl-64/mpa.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpa.h b/sysdeps/ieee754/dbl-64/mpa.h index b8ca297eb7..1e188de4d1 100644 --- a/sysdeps/ieee754/dbl-64/mpa.h +++ b/sysdeps/ieee754/dbl-64/mpa.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -119,36 +119,5 @@ void __dvd (const mp_no *, const mp_no *, mp_no *, int); extern void __mpatan (mp_no *, mp_no *, int); extern void __mpatan2 (mp_no *, mp_no *, mp_no *, int); extern void __mpsqrt (mp_no *, mp_no *, int); -extern void __mpexp (mp_no *, mp_no *, int); extern void __c32 (mp_no *, mp_no *, mp_no *, int); extern int __mpranred (double, mp_no *, int); - -/* Given a power POW, build a multiprecision number 2^POW. */ -static inline void -__pow_mp (int pow, mp_no *y, int p) -{ - int i, rem; - - /* The exponent is E such that E is a factor of 2^24. The remainder (of the - form 2^x) goes entirely into the first digit of the mantissa as it is - always less than 2^24. */ - EY = pow / 24; - rem = pow - EY * 24; - EY++; - - /* If the remainder is negative, it means that POW was negative since - |EY * 24| <= |pow|. Adjust so that REM is positive and still less than - 24 because of which, the mantissa digit is less than 2^24. */ - if (rem < 0) - { - EY--; - rem += 24; - } - /* The sign of any 2^x is always positive. */ - Y[0] = 1; - Y[1] = 1 << rem; - - /* Everything else is 0. */ - for (i = 2; i <= p; i++) - Y[i] = 0; -} diff --git a/sysdeps/ieee754/dbl-64/mpatan.c b/sysdeps/ieee754/dbl-64/mpatan.c index dbd49b2bee..656dc763aa 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.c +++ b/sysdeps/ieee754/dbl-64/mpatan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpatan.h b/sysdeps/ieee754/dbl-64/mpatan.h index 36c02602ac..2e34ce8a4d 100644 --- a/sysdeps/ieee754/dbl-64/mpatan.h +++ b/sysdeps/ieee754/dbl-64/mpatan.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpatan2.c b/sysdeps/ieee754/dbl-64/mpatan2.c index e7de6bc158..6eda7d39eb 100644 --- a/sysdeps/ieee754/dbl-64/mpatan2.c +++ b/sysdeps/ieee754/dbl-64/mpatan2.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpexp.c b/sysdeps/ieee754/dbl-64/mpexp.c deleted file mode 100644 index f17baf2139..0000000000 --- a/sysdeps/ieee754/dbl-64/mpexp.c +++ /dev/null @@ -1,163 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation; either version 2.1 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program; if not, see <http://www.gnu.org/licenses/>. - */ -/*************************************************************************/ -/* MODULE_NAME:mpexp.c */ -/* */ -/* FUNCTIONS: mpexp */ -/* */ -/* FILES NEEDED: mpa.h endian.h mpexp.h */ -/* mpa.c */ -/* */ -/* Multi-Precision exponential function subroutine */ -/* ( for p >= 4, 2**(-55) <= abs(x) <= 1024 ). */ -/*************************************************************************/ - -#include "endian.h" -#include "mpa.h" -#include <assert.h> - -#ifndef SECTION -# define SECTION -#endif - -/* Multi-Precision exponential function subroutine (for p >= 4, - 2**(-55) <= abs(x) <= 1024). */ -void -SECTION -__mpexp (mp_no *x, mp_no *y, int p) -{ - int i, j, k, m, m1, m2, n; - mantissa_t b; - static const int np[33] = - { - 0, 0, 0, 0, 3, 3, 4, 4, 5, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, - 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8 - }; - - static const int m1p[33] = - { - 0, 0, 0, 0, - 17, 23, 23, 28, - 27, 38, 42, 39, - 43, 47, 43, 47, - 50, 54, 57, 60, - 64, 67, 71, 74, - 68, 71, 74, 77, - 70, 73, 76, 78, - 81 - }; - static const int m1np[7][18] = - { - {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {0, 0, 0, 0, 36, 48, 60, 72, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, - {0, 0, 0, 0, 24, 32, 40, 48, 56, 64, 72, 0, 0, 0, 0, 0, 0, 0}, - {0, 0, 0, 0, 17, 23, 29, 35, 41, 47, 53, 59, 65, 0, 0, 0, 0, 0}, - {0, 0, 0, 0, 0, 0, 23, 28, 33, 38, 42, 47, 52, 57, 62, 66, 0, 0}, - {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, 0, 39, 43, 47, 51, 55, 59, 63}, - {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 43, 47, 50, 54} - }; - mp_no mps, mpk, mpt1, mpt2; - - /* Choose m,n and compute a=2**(-m). */ - n = np[p]; - m1 = m1p[p]; - b = X[1]; - m2 = 24 * EX; - for (; b < HALFRAD; m2--) - b *= 2; - if (b == HALFRAD) - { - for (i = 2; i <= p; i++) - { - if (X[i] != 0) - break; - } - if (i == p + 1) - m2--; - } - - m = m1 + m2; - if (__glibc_unlikely (m <= 0)) - { - /* The m1np array which is used to determine if we can reduce the - polynomial expansion iterations, has only 18 elements. Besides, - numbers smaller than those required by p >= 18 should not come here - at all since the fast phase of exp returns 1.0 for anything less - than 2^-55. */ - assert (p < 18); - m = 0; - for (i = n - 1; i > 0; i--, n--) - if (m1np[i][p] + m2 > 0) - break; - } - - /* Compute s=x*2**(-m). Put result in mps. This is the range-reduced input - that we will use to compute e^s. For the final result, simply raise it - to 2^m. */ - __pow_mp (-m, &mpt1, p); - __mul (x, &mpt1, &mps, p); - - /* Compute the Taylor series for e^s: - - 1 + x/1! + x^2/2! + x^3/3! ... - - for N iterations. We compute this as: - - e^x = 1 + (x * n!/1! + x^2 * n!/2! + x^3 * n!/3!) / n! - = 1 + (x * (n!/1! + x * (n!/2! + x * (n!/3! + x ...)))) / n! - - k! is computed on the fly as KF and at the end of the polynomial loop, KF - is n!, which can be used directly. */ - __cpy (&mps, &mpt2, p); - - double kf = 1.0; - - /* Evaluate the rest. The result will be in mpt2. */ - for (k = n - 1; k > 0; k--) - { - /* n! / k! = n * (n - 1) ... * (n - k + 1) */ - kf *= k + 1; - - __dbl_mp (kf, &mpk, p); - __add (&mpt2, &mpk, &mpt1, p); - __mul (&mps, &mpt1, &mpt2, p); - } - __dbl_mp (kf, &mpk, p); - __dvd (&mpt2, &mpk, &mpt1, p); - __add (&__mpone, &mpt1, &mpt2, p); - - /* Raise polynomial value to the power of 2**m. Put result in y. */ - for (k = 0, j = 0; k < m;) - { - __sqr (&mpt2, &mpt1, p); - k++; - if (k == m) - { - j = 1; - break; - } - __sqr (&mpt1, &mpt2, p); - k++; - } - if (j) - __cpy (&mpt1, y, p); - else - __cpy (&mpt2, y, p); - return; -} diff --git a/sysdeps/ieee754/dbl-64/mplog.c b/sysdeps/ieee754/dbl-64/mplog.c deleted file mode 100644 index b297153b10..0000000000 --- a/sysdeps/ieee754/dbl-64/mplog.c +++ /dev/null @@ -1,65 +0,0 @@ -/* - * IBM Accurate Mathematical Library - * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation; either version 2.1 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program; if not, see <http://www.gnu.org/licenses/>. - */ -/************************************************************************/ -/* */ -/* MODULE_NAME:mplog.c */ -/* */ -/* FUNCTIONS: mplog */ -/* */ -/* FILES NEEDED: endian.h mpa.h mplog.h */ -/* mpexp.c */ -/* */ -/* Multi-Precision logarithm function subroutine (for precision p >= 4, */ -/* 2**(-1024) < x < 2**1024) and x is outside of the interval */ -/* [1-2**(-54),1+2**(-54)]. Upon entry, x should be set to the */ -/* multi-precision value of the input and y should be set into a multi- */ -/* precision value of an approximation of log(x) with relative error */ -/* bound of at most 2**(-52). The routine improves the accuracy of y. */ -/* */ -/************************************************************************/ -#include "endian.h" -#include "mpa.h" - -void -__mplog (mp_no *x, mp_no *y, int p) -{ - int i, m; - static const int mp[33] = - { - 0, 0, 0, 0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 - }; - mp_no mpt1, mpt2; - - /* Choose m. */ - m = mp[p]; - - /* Perform m newton iterations to solve for y: exp(y) - x = 0. The - iterations formula is: y(n + 1) = y(n) + (x * exp(-y(n)) - 1). */ - __cpy (y, &mpt1, p); - for (i = 0; i < m; i++) - { - mpt1.d[0] = -mpt1.d[0]; - __mpexp (&mpt1, &mpt2, p); - __mul (x, &mpt2, &mpt1, p); - __sub (&mpt1, &__mpone, &mpt2, p); - __add (y, &mpt2, &mpt1, p); - __cpy (&mpt1, y, p); - } -} diff --git a/sysdeps/ieee754/dbl-64/mpn2dbl.c b/sysdeps/ieee754/dbl-64/mpn2dbl.c index 1050ec3b0e..693448fc8d 100644 --- a/sysdeps/ieee754/dbl-64/mpn2dbl.c +++ b/sysdeps/ieee754/dbl-64/mpn2dbl.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1995-2016 Free Software Foundation, Inc. +/* Copyright (C) 1995-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.c b/sysdeps/ieee754/dbl-64/mpsqrt.c index a77d3938a9..3ce4b9cda7 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.c +++ b/sysdeps/ieee754/dbl-64/mpsqrt.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mpsqrt.h b/sysdeps/ieee754/dbl-64/mpsqrt.h index 0f1cf58fb2..05a907bedf 100644 --- a/sysdeps/ieee754/dbl-64/mpsqrt.h +++ b/sysdeps/ieee754/dbl-64/mpsqrt.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mptan.c b/sysdeps/ieee754/dbl-64/mptan.c index 280c4909d6..26dbd5f601 100644 --- a/sysdeps/ieee754/dbl-64/mptan.c +++ b/sysdeps/ieee754/dbl-64/mptan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/mydefs.h b/sysdeps/ieee754/dbl-64/mydefs.h index b1d67d4445..3e15d5bdd6 100644 --- a/sysdeps/ieee754/dbl-64/mydefs.h +++ b/sysdeps/ieee754/dbl-64/mydefs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/powtwo.tbl b/sysdeps/ieee754/dbl-64/powtwo.tbl index 2d828f6022..518fd6cd4b 100644 --- a/sysdeps/ieee754/dbl-64/powtwo.tbl +++ b/sysdeps/ieee754/dbl-64/powtwo.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/root.tbl b/sysdeps/ieee754/dbl-64/root.tbl index 61c859248c..a458e25ee2 100644 --- a/sysdeps/ieee754/dbl-64/root.tbl +++ b/sysdeps/ieee754/dbl-64/root.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/s_asinh.c b/sysdeps/ieee754/dbl-64/s_asinh.c index 9193301b5e..192ff8594d 100644 --- a/sysdeps/ieee754/dbl-64/s_asinh.c +++ b/sysdeps/ieee754/dbl-64/s_asinh.c @@ -24,6 +24,8 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> static const double one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ @@ -54,19 +56,15 @@ __asinh (double x) double xa = fabs (x); if (ix > 0x40000000) /* 2**28 > |x| > 2.0 */ { - w = __ieee754_log (2.0 * xa + one / (__ieee754_sqrt (xa * xa + one) + + w = __ieee754_log (2.0 * xa + one / (sqrt (xa * xa + one) + xa)); } else /* 2.0 > |x| > 2**-28 */ { double t = xa * xa; - w = __log1p (xa + t / (one + __ieee754_sqrt (one + t))); + w = __log1p (xa + t / (one + sqrt (one + t))); } } return __copysign (w, x); } -weak_alias (__asinh, asinh) -#ifdef NO_LONG_DOUBLE -strong_alias (__asinh, __asinhl) -weak_alias (__asinh, asinhl) -#endif +libm_alias_double (__asinh, asinh) diff --git a/sysdeps/ieee754/dbl-64/s_atan.c b/sysdeps/ieee754/dbl-64/s_atan.c index 780c3ff17a..38db092d04 100644 --- a/sysdeps/ieee754/dbl-64/s_atan.c +++ b/sysdeps/ieee754/dbl-64/s_atan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -43,8 +43,10 @@ #include "atnat.h" #include <fenv.h> #include <float.h> +#include <libm-alias-double.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> #include <stap-probe.h> void __mpatan (mp_no *, mp_no *, int); /* see definition in mpatan.c */ @@ -61,7 +63,7 @@ __signArctan (double x, double y) /* An ultimate atan() routine. Given an IEEE double machine number x, */ /* routine computes the correctly rounded (to nearest) value of atan(x). */ double -atan (double x) +__atan (double x) { double cor, s1, ss1, s2, ss2, t1, t2, t3, t7, t8, t9, t10, u, u2, u3, v, vv, w, ww, y, yy, z, zz; @@ -323,6 +325,6 @@ atanMp (double x, const int pr[]) return y1; /*if impossible to do exact computing */ } -#ifdef NO_LONG_DOUBLE -weak_alias (atan, atanl) +#ifndef __atan +libm_alias_double (__atan, atan) #endif diff --git a/sysdeps/ieee754/dbl-64/s_cbrt.c b/sysdeps/ieee754/dbl-64/s_cbrt.c index 647f30bc71..6cd55dc0d6 100644 --- a/sysdeps/ieee754/dbl-64/s_cbrt.c +++ b/sysdeps/ieee754/dbl-64/s_cbrt.c @@ -1,5 +1,5 @@ /* Compute cubic root of double value. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Dirk Alboth <dirka@uni-paderborn.de> and Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #define CBRT2 1.2599210498948731648 /* 2^(1/3) */ @@ -69,8 +70,4 @@ __cbrt (double x) return __ldexp (x > 0.0 ? ym : -ym, xe / 3); } -weak_alias (__cbrt, cbrt) -#ifdef NO_LONG_DOUBLE -strong_alias (__cbrt, __cbrtl) -weak_alias (__cbrt, cbrtl) -#endif +libm_alias_double (__cbrt, cbrt) diff --git a/sysdeps/ieee754/dbl-64/s_ceil.c b/sysdeps/ieee754/dbl-64/s_ceil.c index b2154b407d..5a7434c737 100644 --- a/sysdeps/ieee754/dbl-64/s_ceil.c +++ b/sysdeps/ieee754/dbl-64/s_ceil.c @@ -15,27 +15,23 @@ * Return x rounded toward -inf to integral value * Method: * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). */ #include <math.h> #include <math_private.h> - -static const double huge = 1.0e300; +#include <libm-alias-double.h> double __ceil (double x) { int32_t i0, i1, j0; - u_int32_t i, j; + uint32_t i, j; EXTRACT_WORDS (i0, i1, x); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; if (j0 < 20) { - if (j0 < 0) /* raise inexact if x != 0 */ + if (j0 < 0) { - math_force_eval (huge + x); /* return 0*sign(x) if |x|<1 */ if (i0 < 0) { @@ -51,7 +47,6 @@ __ceil (double x) i = (0x000fffff) >> j0; if (((i0 & i) | i1) == 0) return x; /* x is integral */ - math_force_eval (huge + x); /* raise inexact flag */ if (i0 > 0) i0 += (0x00100000) >> j0; i0 &= (~i); i1 = 0; @@ -66,10 +61,9 @@ __ceil (double x) } else { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + i = ((uint32_t) (0xffffffff)) >> (j0 - 20); if ((i1 & i) == 0) return x; /* x is integral */ - math_force_eval (huge + x); /* raise inexact flag */ if (i0 > 0) { if (j0 == 20) @@ -88,9 +82,5 @@ __ceil (double x) return x; } #ifndef __ceil -weak_alias (__ceil, ceil) -# ifdef NO_LONG_DOUBLE -strong_alias (__ceil, __ceill) -weak_alias (__ceil, ceill) -# endif +libm_alias_double (__ceil, ceil) #endif diff --git a/sysdeps/ieee754/dbl-64/s_copysign.c b/sysdeps/ieee754/dbl-64/s_copysign.c index 9caf24e8f2..ab81d732ab 100644 --- a/sysdeps/ieee754/dbl-64/s_copysign.c +++ b/sysdeps/ieee754/dbl-64/s_copysign.c @@ -22,18 +22,15 @@ static char rcsid[] = "$NetBSD: s_copysign.c,v 1.8 1995/05/10 20:46:57 jtc Exp $ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> double __copysign (double x, double y) { - u_int32_t hx, hy; + uint32_t hx, hy; GET_HIGH_WORD (hx, x); GET_HIGH_WORD (hy, y); SET_HIGH_WORD (x, (hx & 0x7fffffff) | (hy & 0x80000000)); return x; } -weak_alias (__copysign, copysign) -#ifdef NO_LONG_DOUBLE -strong_alias (__copysign, __copysignl) -weak_alias (__copysign, copysignl) -#endif +libm_alias_double (__copysign, copysign) diff --git a/sysdeps/ieee754/dbl-64/s_erf.c b/sysdeps/ieee754/dbl-64/s_erf.c index b4975a8af8..5f820604f8 100644 --- a/sysdeps/ieee754/dbl-64/s_erf.c +++ b/sysdeps/ieee754/dbl-64/s_erf.c @@ -115,7 +115,10 @@ static char rcsid[] = "$NetBSD: s_erf.c,v 1.8 1995/05/10 20:47:05 jtc Exp $"; #include <errno.h> #include <float.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> #include <fix-int-fp-convert-zero.h> static const double @@ -201,7 +204,7 @@ __erf (double x) ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) /* erf(nan)=nan */ { - i = ((u_int32_t) hx >> 31) << 1; + i = ((uint32_t) hx >> 31) << 1; return (double) (1 - i) + one / x; /* erf(+-inf)=+-1 */ } @@ -294,11 +297,7 @@ __erf (double x) else return r / x - one; } -weak_alias (__erf, erf) -#ifdef NO_LONG_DOUBLE -strong_alias (__erf, __erfl) -weak_alias (__erf, erfl) -#endif +libm_alias_double (__erf, erf) double __erfc (double x) @@ -309,7 +308,7 @@ __erfc (double x) ix = hx & 0x7fffffff; if (ix >= 0x7ff00000) /* erfc(nan)=nan */ { /* erfc(+-inf)=0,2 */ - double ret = (double) (((u_int32_t) hx >> 31) << 1) + one / x; + double ret = (double) (((uint32_t) hx >> 31) << 1) + one / x; if (FIX_INT_FP_CONVERT_ZERO && ret == 0.0) return 0.0; return ret; @@ -421,8 +420,4 @@ __erfc (double x) return two - tiny; } } -weak_alias (__erfc, erfc) -#ifdef NO_LONG_DOUBLE -strong_alias (__erfc, __erfcl) -weak_alias (__erfc, erfcl) -#endif +libm_alias_double (__erfc, erfc) diff --git a/sysdeps/ieee754/dbl-64/s_expm1.c b/sysdeps/ieee754/dbl-64/s_expm1.c index 54d771007a..8f1c95bd04 100644 --- a/sysdeps/ieee754/dbl-64/s_expm1.c +++ b/sysdeps/ieee754/dbl-64/s_expm1.c @@ -111,7 +111,10 @@ #include <errno.h> #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> #define one Q[0] static const double huge = 1.0e+300, @@ -132,7 +135,7 @@ __expm1 (double x) { double y, hi, lo, c, t, e, hxs, hfx, r1, h2, h4, R1, R2, R3; int32_t k, xsb; - u_int32_t hx; + uint32_t hx; GET_HIGH_WORD (hx, x); xsb = hx & 0x80000000; /* sign bit of x */ @@ -149,7 +152,7 @@ __expm1 (double x) { if (hx >= 0x7ff00000) { - u_int32_t low; + uint32_t low; GET_LOW_WORD (low, x); if (((hx & 0xfffff) | low) != 0) return x + x; /* NaN */ @@ -228,7 +231,7 @@ __expm1 (double x) } if (k <= -2 || k > 56) /* suffice to return exp(x)-1 */ { - u_int32_t high; + uint32_t high; y = one - (e - x); GET_HIGH_WORD (high, y); SET_HIGH_WORD (y, high + (k << 20)); /* add k to y's exponent */ @@ -237,7 +240,7 @@ __expm1 (double x) t = one; if (k < 20) { - u_int32_t high; + uint32_t high; SET_HIGH_WORD (t, 0x3ff00000 - (0x200000 >> k)); /* t=1-2^-k */ y = t - (e - x); GET_HIGH_WORD (high, y); @@ -245,7 +248,7 @@ __expm1 (double x) } else { - u_int32_t high; + uint32_t high; SET_HIGH_WORD (t, ((0x3ff - k) << 20)); /* 2^-k */ y = x - (e + t); y += one; @@ -255,8 +258,4 @@ __expm1 (double x) } return y; } -weak_alias (__expm1, expm1) -#ifdef NO_LONG_DOUBLE -strong_alias (__expm1, __expm1l) -weak_alias (__expm1, expm1l) -#endif +libm_alias_double (__expm1, expm1) diff --git a/sysdeps/ieee754/dbl-64/s_f32xaddf64.c b/sysdeps/ieee754/dbl-64/s_f32xaddf64.c new file mode 100644 index 0000000000..6e20a9e939 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_f32xaddf64.c @@ -0,0 +1,30 @@ +/* Add _Float64 values, converting the result to _Float32x. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define daddl __hide_daddl +#include <math.h> +#undef daddl + +#include <math-narrow.h> + +_Float32x +__f32xaddf64 (_Float64 x, _Float64 y) +{ + NARROW_ADD_TRIVIAL (x, y, _Float32x); +} +libm_alias_float32x_float64 (add) diff --git a/sysdeps/ieee754/dbl-64/s_f32xdivf64.c b/sysdeps/ieee754/dbl-64/s_f32xdivf64.c new file mode 100644 index 0000000000..24dc25dc0c --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_f32xdivf64.c @@ -0,0 +1,30 @@ +/* Divide _Float64 values, converting the result to _Float32x. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define ddivl __hide_ddivl +#include <math.h> +#undef ddivl + +#include <math-narrow.h> + +_Float32x +__f32xdivf64 (_Float64 x, _Float64 y) +{ + NARROW_DIV_TRIVIAL (x, y, _Float32x); +} +libm_alias_float32x_float64 (div) diff --git a/sysdeps/ieee754/dbl-64/s_f32xmulf64.c b/sysdeps/ieee754/dbl-64/s_f32xmulf64.c new file mode 100644 index 0000000000..f899c84567 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_f32xmulf64.c @@ -0,0 +1,30 @@ +/* Multiply _Float64 values, converting the result to _Float32x. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define dmull __hide_dmull +#include <math.h> +#undef dmull + +#include <math-narrow.h> + +_Float32x +__f32xmulf64 (_Float64 x, _Float64 y) +{ + NARROW_MUL_TRIVIAL (x, y, _Float32x); +} +libm_alias_float32x_float64 (mul) diff --git a/sysdeps/ieee754/dbl-64/s_f32xsubf64.c b/sysdeps/ieee754/dbl-64/s_f32xsubf64.c new file mode 100644 index 0000000000..568f43a880 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_f32xsubf64.c @@ -0,0 +1,30 @@ +/* Subtract _Float64 values, converting the result to _Float32x. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define dsubl __hide_dsubl +#include <math.h> +#undef dsubl + +#include <math-narrow.h> + +_Float32x +__f32xsubf64 (_Float64 x, _Float64 y) +{ + NARROW_SUB_TRIVIAL (x, y, _Float32x); +} +libm_alias_float32x_float64 (sub) diff --git a/sysdeps/ieee754/dbl-64/s_fabs.c b/sysdeps/ieee754/dbl-64/s_fabs.c index 73c09a269e..8232183324 100644 --- a/sysdeps/ieee754/dbl-64/s_fabs.c +++ b/sysdeps/ieee754/dbl-64/s_fabs.c @@ -19,14 +19,11 @@ static char rcsid[] = "$NetBSD: s_fabs.c,v 1.7 1995/05/10 20:47:13 jtc Exp $"; */ #include <math.h> +#include <libm-alias-double.h> double __fabs (double x) { return __builtin_fabs (x); } -weak_alias (__fabs, fabs) -#ifdef NO_LONG_DOUBLE -strong_alias (__fabs, __fabsl) -weak_alias (__fabs, fabsl) -#endif +libm_alias_double (__fabs, fabs) diff --git a/sysdeps/ieee754/dbl-64/s_fadd.c b/sysdeps/ieee754/dbl-64/s_fadd.c new file mode 100644 index 0000000000..5ecb435497 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fadd.c @@ -0,0 +1,34 @@ +/* Add double values, narrowing the result to float. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define f32addf64 __hide_f32addf64 +#define f32addf32x __hide_f32addf32x +#define faddl __hide_faddl +#include <math.h> +#undef f32addf64 +#undef f32addf32x +#undef faddl + +#include <math-narrow.h> + +float +__fadd (double x, double y) +{ + NARROW_ADD_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1); +} +libm_alias_float_double (add) diff --git a/sysdeps/ieee754/dbl-64/s_fdiv.c b/sysdeps/ieee754/dbl-64/s_fdiv.c new file mode 100644 index 0000000000..cdd2649df0 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fdiv.c @@ -0,0 +1,34 @@ +/* Divide double values, narrowing the result to float. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define f32divf64 __hide_f32divf64 +#define f32divf32x __hide_f32divf32x +#define fdivl __hide_fdivl +#include <math.h> +#undef f32divf64 +#undef f32divf32x +#undef fdivl + +#include <math-narrow.h> + +float +__fdiv (double x, double y) +{ + NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1); +} +libm_alias_float_double (div) diff --git a/sysdeps/ieee754/dbl-64/s_finite.c b/sysdeps/ieee754/dbl-64/s_finite.c index 69141db75d..da1519b1d0 100644 --- a/sysdeps/ieee754/dbl-64/s_finite.c +++ b/sysdeps/ieee754/dbl-64/s_finite.c @@ -21,6 +21,7 @@ static char rcsid[] = "$NetBSD: s_finite.c,v 1.8 1995/05/10 20:47:17 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> #undef __finite @@ -33,12 +34,12 @@ int FINITE(double x) { int32_t hx; GET_HIGH_WORD (hx, x); - return (int) ((u_int32_t) ((hx & 0x7ff00000) - 0x7ff00000) >> 31); + return (int) ((uint32_t) ((hx & 0x7ff00000) - 0x7ff00000) >> 31); } hidden_def (__finite) weak_alias (__finite, finite) #ifdef NO_LONG_DOUBLE -# ifdef LDBL_CLASSIFY_COMPAT +# if LDBL_CLASSIFY_COMPAT # if SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __finite, __finitel, GLIBC_2_0); # endif diff --git a/sysdeps/ieee754/dbl-64/s_floor.c b/sysdeps/ieee754/dbl-64/s_floor.c index bd6afa72e8..f27c6f3ad2 100644 --- a/sysdeps/ieee754/dbl-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/s_floor.c @@ -15,27 +15,24 @@ * Return x rounded toward -inf to integral value * Method: * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). */ #include <math.h> #include <math_private.h> - -static const double huge = 1.0e300; +#include <libm-alias-double.h> double __floor (double x) { int32_t i0, i1, j0; - u_int32_t i, j; + uint32_t i, j; EXTRACT_WORDS (i0, i1, x); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; if (j0 < 20) { - if (j0 < 0) /* raise inexact if x != 0 */ + if (j0 < 0) { - math_force_eval (huge + x); /* return 0*sign(x) if |x|<1 */ + /* return 0*sign(x) if |x|<1 */ if (i0 >= 0) { i0 = i1 = 0; @@ -50,7 +47,6 @@ __floor (double x) i = (0x000fffff) >> j0; if (((i0 & i) | i1) == 0) return x; /* x is integral */ - math_force_eval (huge + x); /* raise inexact flag */ if (i0 < 0) i0 += (0x00100000) >> j0; i0 &= (~i); i1 = 0; @@ -65,10 +61,9 @@ __floor (double x) } else { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + i = ((uint32_t) (0xffffffff)) >> (j0 - 20); if ((i1 & i) == 0) return x; /* x is integral */ - math_force_eval (huge + x); /* raise inexact flag */ if (i0 < 0) { if (j0 == 20) @@ -87,9 +82,5 @@ __floor (double x) return x; } #ifndef __floor -weak_alias (__floor, floor) -# ifdef NO_LONG_DOUBLE -strong_alias (__floor, __floorl) -weak_alias (__floor, floorl) -# endif +libm_alias_double (__floor, floor) #endif diff --git a/sysdeps/ieee754/dbl-64/s_fma.c b/sysdeps/ieee754/dbl-64/s_fma.c index a3492434e4..57c7b5dfc2 100644 --- a/sysdeps/ieee754/dbl-64/s_fma.c +++ b/sysdeps/ieee754/dbl-64/s_fma.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2016 Free Software Foundation, Inc. + Copyright (C) 2010-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -21,7 +21,9 @@ #include <math.h> #include <fenv.h> #include <ieee754.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-double.h> #include <tininess.h> /* This implementation uses rounding to odd to avoid problems with @@ -292,10 +294,5 @@ __fma (double x, double y, double z) } } #ifndef __fma -weak_alias (__fma, fma) -#endif - -#ifdef NO_LONG_DOUBLE -strong_alias (__fma, __fmal) -weak_alias (__fmal, fmal) +libm_alias_double (__fma, fma) #endif diff --git a/sysdeps/ieee754/dbl-64/s_fmaf.c b/sysdeps/ieee754/dbl-64/s_fmaf.c index 582c205572..5c8b22ac16 100644 --- a/sysdeps/ieee754/dbl-64/s_fmaf.c +++ b/sysdeps/ieee754/dbl-64/s_fmaf.c @@ -1,5 +1,5 @@ /* Compute x * y + z as ternary operation. - Copyright (C) 2010-2016 Free Software Foundation, Inc. + Copyright (C) 2010-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Jakub Jelinek <jakub@redhat.com>, 2010. @@ -20,7 +20,9 @@ #include <math.h> #include <fenv.h> #include <ieee754.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-float.h> /* This implementation relies on double being more than twice as precise as float and uses rounding to odd in order to avoid problems @@ -60,5 +62,5 @@ __fmaf (float x, float y, float z) return (float) u.d; } #ifndef __fmaf -weak_alias (__fmaf, fmaf) +libm_alias_float (__fma, fma) #endif diff --git a/sysdeps/ieee754/dbl-64/s_fmul.c b/sysdeps/ieee754/dbl-64/s_fmul.c new file mode 100644 index 0000000000..ad9ab7e94d --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fmul.c @@ -0,0 +1,34 @@ +/* Multiply double values, narrowing the result to float. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define f32mulf64 __hide_f32mulf64 +#define f32mulf32x __hide_f32mulf32x +#define fmull __hide_fmull +#include <math.h> +#undef f32mulf64 +#undef f32mulf32x +#undef fmull + +#include <math-narrow.h> + +float +__fmul (double x, double y) +{ + NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1); +} +libm_alias_float_double (mul) diff --git a/sysdeps/ieee754/dbl-64/s_fpclassify.c b/sysdeps/ieee754/dbl-64/s_fpclassify.c index 176e29fccb..1c2d1938f9 100644 --- a/sysdeps/ieee754/dbl-64/s_fpclassify.c +++ b/sysdeps/ieee754/dbl-64/s_fpclassify.c @@ -1,5 +1,5 @@ /* Return classification value corresponding to argument. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -25,7 +25,7 @@ int __fpclassify (double x) { - u_int32_t hx, lx; + uint32_t hx, lx; int retval = FP_NORMAL; EXTRACT_WORDS (hx, lx, x); diff --git a/sysdeps/ieee754/dbl-64/s_frexp.c b/sysdeps/ieee754/dbl-64/s_frexp.c index 1b8d8500ba..c96a869665 100644 --- a/sysdeps/ieee754/dbl-64/s_frexp.c +++ b/sysdeps/ieee754/dbl-64/s_frexp.c @@ -26,6 +26,7 @@ static char rcsid[] = "$NetBSD: s_frexp.c,v 1.9 1995/05/10 20:47:24 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> static const double two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */ @@ -38,7 +39,7 @@ __frexp (double x, int *eptr) ix = 0x7fffffff & hx; *eptr = 0; if (ix >= 0x7ff00000 || ((ix | lx) == 0)) - return x; /* 0,inf,nan */ + return x + x; /* 0,inf,nan */ if (ix < 0x00100000) /* subnormal */ { x *= two54; @@ -51,8 +52,4 @@ __frexp (double x, int *eptr) SET_HIGH_WORD (x, hx); return x; } -weak_alias (__frexp, frexp) -#ifdef NO_LONG_DOUBLE -strong_alias (__frexp, __frexpl) -weak_alias (__frexp, frexpl) -#endif +libm_alias_double (__frexp, frexp) diff --git a/sysdeps/ieee754/dbl-64/s_fromfp.c b/sysdeps/ieee754/dbl-64/s_fromfp.c new file mode 100644 index 0000000000..30572b2a9b --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fromfp.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 0 +#define FUNC __fromfp +#include <s_fromfp_main.c> +libm_alias_double (__fromfp, fromfp) diff --git a/sysdeps/ieee754/dbl-64/s_fromfp_main.c b/sysdeps/ieee754/dbl-64/s_fromfp_main.c new file mode 100644 index 0000000000..dfb70b434e --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fromfp_main.c @@ -0,0 +1,83 @@ +/* Round to integer type. dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <errno.h> +#include <fenv.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <stdbool.h> +#include <stdint.h> + +#define BIAS 0x3ff +#define MANT_DIG 53 + +#if UNSIGNED +# define RET_TYPE uintmax_t +#else +# define RET_TYPE intmax_t +#endif + +#include <fromfp.h> + +RET_TYPE +FUNC (double x, int round, unsigned int width) +{ + if (width > INTMAX_WIDTH) + width = INTMAX_WIDTH; + uint64_t ix; + EXTRACT_WORDS64 (ix, x); + bool negative = (ix & 0x8000000000000000ULL) != 0; + if (width == 0) + return fromfp_domain_error (negative, width); + ix &= 0x7fffffffffffffffULL; + if (ix == 0) + return 0; + int exponent = ix >> (MANT_DIG - 1); + exponent -= BIAS; + int max_exponent = fromfp_max_exponent (negative, width); + if (exponent > max_exponent) + return fromfp_domain_error (negative, width); + + ix &= ((1ULL << (MANT_DIG - 1)) - 1); + ix |= 1ULL << (MANT_DIG - 1); + uintmax_t uret; + bool half_bit, more_bits; + if (exponent >= MANT_DIG - 1) + { + uret = ix; + uret <<= exponent - (MANT_DIG - 1); + half_bit = false; + more_bits = false; + } + else if (exponent >= -1) + { + uint64_t h = 1ULL << (MANT_DIG - 2 - exponent); + half_bit = (ix & h) != 0; + more_bits = (ix & (h - 1)) != 0; + uret = ix >> (MANT_DIG - 1 - exponent); + } + else + { + uret = 0; + half_bit = false; + more_bits = true; + } + return fromfp_round_and_return (negative, uret, half_bit, more_bits, round, + exponent, max_exponent, width); +} diff --git a/sysdeps/ieee754/dbl-64/s_fromfpx.c b/sysdeps/ieee754/dbl-64/s_fromfpx.c new file mode 100644 index 0000000000..b7a0d59a05 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fromfpx.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 1 +#define FUNC __fromfpx +#include <s_fromfp_main.c> +libm_alias_double (__fromfpx, fromfpx) diff --git a/sysdeps/ieee754/dbl-64/s_fsub.c b/sysdeps/ieee754/dbl-64/s_fsub.c new file mode 100644 index 0000000000..c3b911e801 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_fsub.c @@ -0,0 +1,34 @@ +/* Subtract double values, narrowing the result to float. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#define f32subf64 __hide_f32subf64 +#define f32subf32x __hide_f32subf32x +#define fsubl __hide_fsubl +#include <math.h> +#undef f32subf64 +#undef f32subf32x +#undef fsubl + +#include <math-narrow.h> + +float +__fsub (double x, double y) +{ + NARROW_SUB_ROUND_TO_ODD (x, y, float, union ieee754_double, , mantissa1); +} +libm_alias_float_double (sub) diff --git a/sysdeps/ieee754/dbl-64/w_exp.c b/sysdeps/ieee754/dbl-64/s_getpayload.c index 595ed992d2..7b995e94ee 100644 --- a/sysdeps/ieee754/dbl-64/w_exp.c +++ b/sysdeps/ieee754/dbl-64/s_getpayload.c @@ -1,6 +1,6 @@ -/* Copyright (C) 2011-2016 Free Software Foundation, Inc. +/* Get NaN payload. dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. - Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. The GNU C Library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public @@ -16,23 +16,21 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <fix-int-fp-convert-zero.h> #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> +#include <stdint.h> -/* wrapper exp */ double -__exp (double x) +__getpayload (const double *x) { - double z = __ieee754_exp (x); - if (__builtin_expect (!isfinite (z) || z == 0, 0) - && isfinite (x) && _LIB_VERSION != _IEEE_) - return __kernel_standard (x, x, 6 + !!signbit (x)); - - return z; + uint32_t hx, lx; + EXTRACT_WORDS (hx, lx, *x); + hx &= 0x7ffff; + uint64_t ix = ((uint64_t) hx << 32) | lx; + if (FIX_INT_FP_CONVERT_ZERO && ix == 0) + return 0.0f; + return (double) ix; } -hidden_def (__exp) -weak_alias (__exp, exp) -#ifdef NO_LONG_DOUBLE -strong_alias (__exp, __expl) -weak_alias (__exp, expl) -#endif +libm_alias_double (__getpayload, getpayload) diff --git a/sysdeps/ieee754/dbl-64/s_isinf.c b/sysdeps/ieee754/dbl-64/s_isinf.c index c0ad54538a..93eb65c147 100644 --- a/sysdeps/ieee754/dbl-64/s_isinf.c +++ b/sysdeps/ieee754/dbl-64/s_isinf.c @@ -15,6 +15,7 @@ static char rcsid[] = "$NetBSD: s_isinf.c,v 1.3 1995/05/11 23:20:14 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> int @@ -29,7 +30,7 @@ __isinf (double x) hidden_def (__isinf) weak_alias (__isinf, isinf) #ifdef NO_LONG_DOUBLE -# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +# if LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __isinf, __isinfl, GLIBC_2_0); # endif weak_alias (__isinf, isinfl) diff --git a/sysdeps/ieee754/dbl-64/s_isnan.c b/sysdeps/ieee754/dbl-64/s_isnan.c index 2174d988d8..82723eeb8a 100644 --- a/sysdeps/ieee754/dbl-64/s_isnan.c +++ b/sysdeps/ieee754/dbl-64/s_isnan.c @@ -21,6 +21,7 @@ static char rcsid[] = "$NetBSD: s_isnan.c,v 1.8 1995/05/10 20:47:36 jtc Exp $"; #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> #undef __isnan @@ -30,14 +31,14 @@ __isnan (double x) int32_t hx, lx; EXTRACT_WORDS (hx, lx, x); hx &= 0x7fffffff; - hx |= (u_int32_t) (lx | (-lx)) >> 31; + hx |= (uint32_t) (lx | (-lx)) >> 31; hx = 0x7ff00000 - hx; - return (int) (((u_int32_t) hx) >> 31); + return (int) (((uint32_t) hx) >> 31); } hidden_def (__isnan) weak_alias (__isnan, isnan) #ifdef NO_LONG_DOUBLE -# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +# if LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __isnan, __isnanl, GLIBC_2_0); # endif weak_alias (__isnan, isnanl) diff --git a/sysdeps/ieee754/dbl-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/s_issignaling.c index 4b93d6ae60..05aa3f2fdd 100644 --- a/sysdeps/ieee754/dbl-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/s_issignaling.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,19 +18,20 @@ #include <math.h> #include <math_private.h> +#include <nan-high-order-bit.h> int __issignaling (double x) { -#ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN - u_int32_t hxi; +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + uint32_t hxi; GET_HIGH_WORD (hxi, x); /* We only have to care about the high-order bit of x's significand, because having it set (sNaN) already makes the significand different from that used to designate infinity. */ return (hxi & 0x7ff80000) == 0x7ff80000; #else - u_int32_t hxi, lxi; + uint32_t hxi, lxi; EXTRACT_WORDS (hxi, lxi, x); /* To keep the following comparison simple, toggle the quiet/signaling bit, so that it is set for sNaNs. This is inverse to IEEE 754-2008 (as well as diff --git a/sysdeps/ieee754/dbl-64/s_llrint.c b/sysdeps/ieee754/dbl-64/s_llrint.c index ff40ee626d..8159706f90 100644 --- a/sysdeps/ieee754/dbl-64/s_llrint.c +++ b/sysdeps/ieee754/dbl-64/s_llrint.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,7 +22,9 @@ #include <limits.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libm-alias-double.h> #include <fix-fp-int-convert-overflow.h> static const double two52[2] = @@ -36,7 +38,7 @@ long long int __llrint (double x) { int32_t j0; - u_int32_t i1, i0; + uint32_t i1, i0; long long int result; double w; double t; @@ -96,8 +98,4 @@ __llrint (double x) return sx ? -result : result; } -weak_alias (__llrint, llrint) -#ifdef NO_LONG_DOUBLE -strong_alias (__llrint, __llrintl) -weak_alias (__llrint, llrintl) -#endif +libm_alias_double (__llrint, llrint) diff --git a/sysdeps/ieee754/dbl-64/s_llround.c b/sysdeps/ieee754/dbl-64/s_llround.c index 1d6a5205f1..1067d19859 100644 --- a/sysdeps/ieee754/dbl-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/s_llround.c @@ -1,5 +1,5 @@ /* Round double value to long long int. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <fix-fp-int-convert-overflow.h> @@ -29,7 +30,7 @@ long long int __llround (double x) { int32_t j0; - u_int32_t i1, i0; + uint32_t i1, i0; long long int result; int sign; @@ -56,7 +57,7 @@ __llround (double x) result = (((long long int) i0 << 32) | i1) << (j0 - 52); else { - u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + uint32_t j = i1 + (0x80000000 >> (j0 - 20)); if (j < i1) ++i0; @@ -84,8 +85,4 @@ __llround (double x) return sign * result; } -weak_alias (__llround, llround) -#ifdef NO_LONG_DOUBLE -strong_alias (__llround, __llroundl) -weak_alias (__llround, llroundl) -#endif +libm_alias_double (__llround, llround) diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c index 340f6377f7..e6476a8260 100644 --- a/sysdeps/ieee754/dbl-64/s_log1p.c +++ b/sysdeps/ieee754/dbl-64/s_log1p.c @@ -80,7 +80,10 @@ #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> +#include <libc-diag.h> static const double ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ @@ -191,5 +194,14 @@ __log1p (double x) if (k == 0) return f - (hfsq - s * (hfsq + R)); else - return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); + { + /* With GCC 7 when compiling with -Os the compiler warns that c + might be used uninitialized. This can't be true because k + must be 0 for c to be uninitialized and we handled that + computation earlier without using c. */ + DIAG_PUSH_NEEDS_COMMENT; + DIAG_IGNORE_Os_NEEDS_COMMENT (7, "-Wmaybe-uninitialized"); + return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f); + DIAG_POP_NEEDS_COMMENT; + } } diff --git a/sysdeps/ieee754/dbl-64/s_logb.c b/sysdeps/ieee754/dbl-64/s_logb.c index 3a26b18f78..a6de1f6e49 100644 --- a/sysdeps/ieee754/dbl-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/s_logb.c @@ -18,6 +18,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <fix-int-fp-convert-zero.h> double @@ -46,7 +47,6 @@ __logb (double x) return 0.0; return (double) (rix - 1023); } -weak_alias (__logb, logb) -#ifdef NO_LONG_DOUBLE -strong_alias (__logb, __logbl) weak_alias (__logb, logbl) +#ifndef __logb +libm_alias_double (__logb, logb) #endif diff --git a/sysdeps/ieee754/dbl-64/s_lrint.c b/sysdeps/ieee754/dbl-64/s_lrint.c index 4e7cbef97b..0e64ae1260 100644 --- a/sysdeps/ieee754/dbl-64/s_lrint.c +++ b/sysdeps/ieee754/dbl-64/s_lrint.c @@ -1,6 +1,6 @@ /* Round argument to nearest integral value according to current rounding direction. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,7 +22,9 @@ #include <limits.h> #include <math.h> +#include <math-narrow-eval.h> #include <math_private.h> +#include <libm-alias-double.h> #include <fix-fp-int-convert-overflow.h> static const double two52[2] = @@ -36,7 +38,7 @@ long int __lrint (double x) { int32_t j0; - u_int32_t i0, i1; + uint32_t i0, i1; double w; double t; long int result; @@ -120,8 +122,4 @@ __lrint (double x) return sx ? -result : result; } -weak_alias (__lrint, lrint) -#ifdef NO_LONG_DOUBLE -strong_alias (__lrint, __lrintl) -weak_alias (__lrint, lrintl) -#endif +libm_alias_double (__lrint, lrint) diff --git a/sysdeps/ieee754/dbl-64/s_lround.c b/sysdeps/ieee754/dbl-64/s_lround.c index aa1275ddc3..70f624eea1 100644 --- a/sysdeps/ieee754/dbl-64/s_lround.c +++ b/sysdeps/ieee754/dbl-64/s_lround.c @@ -1,5 +1,5 @@ /* Round double value to long int. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <fix-fp-int-convert-overflow.h> @@ -29,7 +30,7 @@ long int __lround (double x) { int32_t j0; - u_int32_t i1, i0; + uint32_t i1, i0; long int result; int sign; @@ -56,7 +57,7 @@ __lround (double x) result = ((long int) i0 << (j0 - 20)) | ((long int) i1 << (j0 - 52)); else { - u_int32_t j = i1 + (0x80000000 >> (j0 - 20)); + uint32_t j = i1 + (0x80000000 >> (j0 - 20)); if (j < i1) ++i0; @@ -106,8 +107,4 @@ __lround (double x) return sign * result; } -weak_alias (__lround, lround) -#ifdef NO_LONG_DOUBLE -strong_alias (__lround, __lroundl) -weak_alias (__lround, lroundl) -#endif +libm_alias_double (__lround, lround) diff --git a/sysdeps/ieee754/dbl-64/s_modf.c b/sysdeps/ieee754/dbl-64/s_modf.c index 0a1e13008f..722511c64a 100644 --- a/sysdeps/ieee754/dbl-64/s_modf.c +++ b/sysdeps/ieee754/dbl-64/s_modf.c @@ -21,6 +21,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> static const double one = 1.0; @@ -28,7 +29,7 @@ double __modf (double x, double *iptr) { int32_t i0, i1, j0; - u_int32_t i; + uint32_t i; EXTRACT_WORDS (i0, i1, x); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; /* exponent of x */ if (j0 < 20) /* integer part in high x */ @@ -65,7 +66,7 @@ __modf (double x, double *iptr) } else /* fraction part in low x */ { - i = ((u_int32_t) (0xffffffff)) >> (j0 - 20); + i = ((uint32_t) (0xffffffff)) >> (j0 - 20); if ((i1 & i) == 0) /* x is integral */ { *iptr = x; @@ -79,8 +80,6 @@ __modf (double x, double *iptr) } } } -weak_alias (__modf, modf) -#ifdef NO_LONG_DOUBLE -strong_alias (__modf, __modfl) -weak_alias (__modf, modfl) +#ifndef __modf +libm_alias_double (__modf, modf) #endif diff --git a/sysdeps/ieee754/dbl-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/s_nearbyint.c index dec0c5d6ee..903121d456 100644 --- a/sysdeps/ieee754/dbl-64/s_nearbyint.c +++ b/sysdeps/ieee754/dbl-64/s_nearbyint.c @@ -26,7 +26,9 @@ static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $"; #include <fenv.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-double.h> static const double TWO52[2] = { @@ -48,7 +50,7 @@ __nearbyint (double x) if (j0 < 0) { libc_feholdexcept (&env); - w = TWO52[sx] + x; + w = TWO52[sx] + math_opt_barrier (x); t = w - TWO52[sx]; math_force_eval (t); libc_fesetenv (&env); @@ -65,14 +67,10 @@ __nearbyint (double x) return x; /* x is integral */ } libc_feholdexcept (&env); - w = TWO52[sx] + x; + w = TWO52[sx] + math_opt_barrier (x); t = w - TWO52[sx]; math_force_eval (t); libc_fesetenv (&env); return t; } -weak_alias (__nearbyint, nearbyint) -#ifdef NO_LONG_DOUBLE -strong_alias (__nearbyint, __nearbyintl) -weak_alias (__nearbyint, nearbyintl) -#endif +libm_alias_double (__nearbyint, nearbyint) diff --git a/sysdeps/ieee754/dbl-64/s_nextup.c b/sysdeps/ieee754/dbl-64/s_nextup.c new file mode 100644 index 0000000000..d37a2b6657 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_nextup.c @@ -0,0 +1,56 @@ +/* Return the least floating-point number greater than X. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <float.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> + +/* Return the least floating-point number greater than X. */ +double +__nextup (double x) +{ + int32_t hx, ix; + uint32_t lx; + + EXTRACT_WORDS (hx, lx, x); + ix = hx & 0x7fffffff; + + if (((ix >= 0x7ff00000) && ((ix - 0x7ff00000) | lx) != 0)) /* x is nan. */ + return x + x; + if ((ix | lx) == 0) + return DBL_TRUE_MIN; + if (hx >= 0) + { /* x > 0. */ + if (isinf (x)) + return x; + lx += 1; + if (lx == 0) + hx += 1; + } + else + { /* x < 0. */ + if (lx == 0) + hx -= 1; + lx -= 1; + } + INSERT_WORDS (x, hx, lx); + return x; +} + +libm_alias_double (__nextup, nextup) diff --git a/sysdeps/ieee754/dbl-64/s_remquo.c b/sysdeps/ieee754/dbl-64/s_remquo.c index 16215e1e29..9a3fbe0ec8 100644 --- a/sysdeps/ieee754/dbl-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/s_remquo.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> static const double zero = 0.0; @@ -29,7 +30,7 @@ double __remquo (double x, double y, int *quo) { int32_t hx, hy; - u_int32_t sx, lx, ly; + uint32_t sx, lx, ly; int cquo, qs; EXTRACT_WORDS (hx, lx, x); @@ -108,8 +109,4 @@ __remquo (double x, double y, int *quo) x = -x; return x; } -weak_alias (__remquo, remquo) -#ifdef NO_LONG_DOUBLE -strong_alias (__remquo, __remquol) -weak_alias (__remquo, remquol) -#endif +libm_alias_double (__remquo, remquo) diff --git a/sysdeps/ieee754/dbl-64/s_rint.c b/sysdeps/ieee754/dbl-64/s_rint.c index a9c0d27842..cb0f5ca298 100644 --- a/sysdeps/ieee754/dbl-64/s_rint.c +++ b/sysdeps/ieee754/dbl-64/s_rint.c @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> static const double TWO52[2] = { @@ -59,9 +60,5 @@ __rint (double x) return w - TWO52[sx]; } #ifndef __rint -weak_alias (__rint, rint) -# ifdef NO_LONG_DOUBLE -strong_alias (__rint, __rintl) -weak_alias (__rint, rintl) -# endif +libm_alias_double (__rint, rint) #endif diff --git a/sysdeps/ieee754/dbl-64/s_round.c b/sysdeps/ieee754/dbl-64/s_round.c index aeed6faf0e..fa9e83196e 100644 --- a/sysdeps/ieee754/dbl-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,16 +20,14 @@ #include <math.h> #include <math_private.h> - - -static const double huge = 1.0e300; +#include <libm-alias-double.h> double __round (double x) { int32_t i0, j0; - u_int32_t i1; + uint32_t i1; EXTRACT_WORDS (i0, i1, x); j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; @@ -37,8 +35,6 @@ __round (double x) { if (j0 < 0) { - math_force_eval (huge + x); - i0 &= 0x80000000; if (j0 == -1) i0 |= 0x3ff00000; @@ -46,13 +42,11 @@ __round (double x) } else { - u_int32_t i = 0x000fffff >> j0; + uint32_t i = 0x000fffff >> j0; if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - math_force_eval (huge + x); - /* Raise inexact if x != 0. */ i0 += 0x00080000 >> j0; i0 &= ~i; i1 = 0; @@ -68,15 +62,12 @@ __round (double x) } else { - u_int32_t i = 0xffffffff >> (j0 - 20); + uint32_t i = 0xffffffff >> (j0 - 20); if ((i1 & i) == 0) /* X is integral. */ return x; - math_force_eval (huge + x); - - /* Raise inexact if x != 0. */ - u_int32_t j = i1 + (1 << (51 - j0)); + uint32_t j = i1 + (1 << (51 - j0)); if (j < i1) i0 += 1; i1 = j; @@ -86,8 +77,4 @@ __round (double x) INSERT_WORDS (x, i0, i1); return x; } -weak_alias (__round, round) -#ifdef NO_LONG_DOUBLE -strong_alias (__round, __roundl) -weak_alias (__round, roundl) -#endif +libm_alias_double (__round, round) diff --git a/sysdeps/ieee754/dbl-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/s_roundeven.c new file mode 100644 index 0000000000..1438e81d45 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_roundeven.c @@ -0,0 +1,105 @@ +/* Round to nearest integer value, rounding halfway cases to even. + dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <stdint.h> + +#define BIAS 0x3ff +#define MANT_DIG 53 +#define MAX_EXP (2 * BIAS + 1) + +double +__roundeven (double x) +{ + uint32_t hx, lx, uhx; + EXTRACT_WORDS (hx, lx, x); + uhx = hx & 0x7fffffff; + int exponent = uhx >> (MANT_DIG - 1 - 32); + if (exponent >= BIAS + MANT_DIG - 1) + { + /* Integer, infinity or NaN. */ + if (exponent == MAX_EXP) + /* Infinity or NaN; quiet signaling NaNs. */ + return x + x; + else + return x; + } + else if (exponent >= BIAS + MANT_DIG - 32) + { + /* Not necessarily an integer; integer bit is in low word. + Locate the bits with exponents 0 and -1. */ + int int_pos = (BIAS + MANT_DIG - 1) - exponent; + int half_pos = int_pos - 1; + uint32_t half_bit = 1U << half_pos; + uint32_t int_bit = 1U << int_pos; + if ((lx & (int_bit | (half_bit - 1))) != 0) + { + /* Carry into the exponent works correctly. No need to test + whether HALF_BIT is set. */ + lx += half_bit; + hx += lx < half_bit; + } + lx &= ~(int_bit - 1); + } + else if (exponent == BIAS + MANT_DIG - 33) + { + /* Not necessarily an integer; integer bit is bottom of high + word, half bit is top of low word. */ + if (((hx & 1) | (lx & 0x7fffffff)) != 0) + { + lx += 0x80000000; + hx += lx < 0x80000000; + } + lx = 0; + } + else if (exponent >= BIAS) + { + /* At least 1; not necessarily an integer, integer bit and half + bit are in the high word. Locate the bits with exponents 0 + and -1 (when the unbiased exponent is 0, the bit with + exponent 0 is implicit, but as the bias is odd it is OK to + take it from the low bit of the exponent). */ + int int_pos = (BIAS + MANT_DIG - 33) - exponent; + int half_pos = int_pos - 1; + uint32_t half_bit = 1U << half_pos; + uint32_t int_bit = 1U << int_pos; + if (((hx & (int_bit | (half_bit - 1))) | lx) != 0) + hx += half_bit; + hx &= ~(int_bit - 1); + lx = 0; + } + else if (exponent == BIAS - 1 && (uhx > 0x3fe00000 || lx != 0)) + { + /* Interval (0.5, 1). */ + hx = (hx & 0x80000000) | 0x3ff00000; + lx = 0; + } + else + { + /* Rounds to 0. */ + hx &= 0x80000000; + lx = 0; + } + INSERT_WORDS (x, hx, lx); + return x; +} +hidden_def (__roundeven) +libm_alias_double (__roundeven, roundeven) diff --git a/sysdeps/ieee754/dbl-64/s_setpayload.c b/sysdeps/ieee754/dbl-64/s_setpayload.c new file mode 100644 index 0000000000..0f536956f4 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_setpayload.c @@ -0,0 +1,4 @@ +#define SIG 0 +#define FUNC __setpayload +#include <s_setpayload_main.c> +libm_alias_double (__setpayload, setpayload) diff --git a/sysdeps/ieee754/dbl-64/s_setpayload_main.c b/sysdeps/ieee754/dbl-64/s_setpayload_main.c new file mode 100644 index 0000000000..0e87c2a587 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_setpayload_main.c @@ -0,0 +1,70 @@ +/* Set NaN payload. dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +#define SET_HIGH_BIT (HIGH_ORDER_BIT_IS_SET_FOR_SNAN ? SIG : !SIG) +#define BIAS 0x3ff +#define PAYLOAD_DIG 51 +#define EXPLICIT_MANT_DIG 52 + +int +FUNC (double *x, double payload) +{ + uint32_t hx, lx; + EXTRACT_WORDS (hx, lx, payload); + int exponent = hx >> (EXPLICIT_MANT_DIG - 32); + /* Test if argument is (a) negative or too large; (b) too small, + except for 0 when allowed; (c) not an integer. */ + if (exponent >= BIAS + PAYLOAD_DIG + || (exponent < BIAS && !(SET_HIGH_BIT && hx == 0 && lx == 0))) + { + INSERT_WORDS (*x, 0, 0); + return 1; + } + int shift = BIAS + EXPLICIT_MANT_DIG - exponent; + if (shift < 32 + ? (lx & ((1U << shift) - 1)) != 0 + : (lx != 0 || (hx & ((1U << (shift - 32)) - 1)) != 0)) + { + INSERT_WORDS (*x, 0, 0); + return 1; + } + if (exponent != 0) + { + hx &= (1U << (EXPLICIT_MANT_DIG - 32)) - 1; + hx |= 1U << (EXPLICIT_MANT_DIG - 32); + if (shift >= 32) + { + lx = hx >> (shift - 32); + hx = 0; + } + else if (shift != 0) + { + lx = (lx >> shift) | (hx << (32 - shift)); + hx >>= shift; + } + } + hx |= 0x7ff00000 | (SET_HIGH_BIT ? 0x80000 : 0); + INSERT_WORDS (*x, hx, lx); + return 0; +} diff --git a/sysdeps/ieee754/dbl-64/s_setpayloadsig.c b/sysdeps/ieee754/dbl-64/s_setpayloadsig.c new file mode 100644 index 0000000000..96ea34b070 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_setpayloadsig.c @@ -0,0 +1,4 @@ +#define SIG 1 +#define FUNC __setpayloadsig +#include <s_setpayload_main.c> +libm_alias_double (__setpayloadsig, setpayloadsig) diff --git a/sysdeps/ieee754/dbl-64/s_signbit.c b/sysdeps/ieee754/dbl-64/s_signbit.c index f4ef4d238f..71500b0a2d 100644 --- a/sysdeps/ieee754/dbl-64/s_signbit.c +++ b/sysdeps/ieee754/dbl-64/s_signbit.c @@ -1,5 +1,5 @@ /* Return nonzero value if number is negative. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index ca2532fb63..b369ac9f5b 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -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. */ /* */ @@ -52,6 +41,8 @@ #include "MathLib.h" #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> #include <fenv.h> /* Helper macros to compute sin of the input values. */ @@ -65,35 +56,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,156 +90,69 @@ 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, int n); -static double sloww1 (double x, double dx, double orig, int m, int n); -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 U and X such that U is an index into the - sin/cos table, this macro 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. */ -static double -do_cos (mynumber u, double x, double *corp) +/* 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. */ +static inline double +__always_inline +do_cos (double x, double dx) { - double xx, s, sn, ssn, c, cs, ccs, res, cor; + mynumber u; + + if (x < 0) + dx = -dx; + + u.x = big + fabs (x); + x = fabs (x) - (u.x - big) + dx; + + 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; + return cs + cor; } -/* A more precise variant of DO_COS where the number is partitioned into U, X - and DX. EPS is the adjustment to the correction COR. */ -static double -do_cos_slow (mynumber u, double x, double dx, double eps, double *corp) +/* 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. */ +static inline double +__always_inline +do_sin (double x, double dx) { - 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; - if (cor > 0) - cor = 1.0005 * cor + eps; - else - cor = 1.0005 * cor - eps; - *corp = cor; - return res; -} + 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); -/* Given a number partitioned into U and X and DX such that U is an index into - the sin/cos table, this macro 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. */ -static double -do_sin (mynumber u, double x, double dx, double *corp) -{ - double xx, s, sn, ssn, c, cs, ccs, cor, res; + mynumber u; + + if (x <= 0) + dx = -dx; + u.x = big + fabs (x); + x = fabs (x) - (u.x - big); + + 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 res = do_sin where the number is partitioned into U, X - and DX. EPS is the adjustment to the correction COR. */ -static double -do_sin_slow (mynumber u, double x, double dx, double eps, double *corp) -{ - 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; - if (cor > 0) - cor = 1.0005 * cor + eps; - else - cor = 1.0005 * cor - eps; - *corp = cor; - return res; -} - -/* Reduce range of X and compute sin of a + da. K is the amount by which to - rotate the quadrants. This allows us to use the same routine to compute cos - by simply rotating the quadrants by 1. */ -static inline double -__always_inline -reduce_and_compute (double x, unsigned int k) -{ - double retval = 0, a, da; - unsigned int n = __branred (x, &a, &da); - k = (n + k) % 4; - switch (k) - { - case 0: - if (a * a < 0.01588) - retval = bsloww (a, da, x, n); - else - retval = bsloww1 (a, da, x, n); - break; - case 2: - 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; @@ -281,198 +161,54 @@ reduce_sincos_1 (double x, double *a, double *da) 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 shifting the quadrant N - clockwise. */ -static double -__always_inline -do_sincos_1 (double a, double da, double x, int4 n, int4 k) -{ - double xx, retval, res, cor, y; - mynumber u; - int m; - double eps = fabs (x) * 1.2e-30; - - int k1 = (n + k) & 3; - switch (k1) - { /* quarter of unit circle */ - case 2: - a = -a; - da = -da; - case 0: - xx = a * a; - if (xx < 0.01588) - { - /* Taylor series. */ - res = TAYLOR_SIN (xx, a, da, cor); - cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : sloww (a, da, x, k); - } - else - { - if (a > 0) - m = 1; - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m, k)); - } - break; - - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - retval = ((res == res + cor) ? ((k1 & 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 b, db, t1, t2; + t1 = xn * pp3; + t2 = y - t1; + db = (y - t2) - t1; - 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; + 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 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, int4 k) +do_sincos (double a, double da, int4 n) { - double res, retval, cor, xx; - mynumber u; - - double eps = 1.0e-24; - - k = (n + k) & 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 = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps; - retval = (res == res + cor) ? res : bsloww (a, da, x, n); - } - else - { - double t, db, y; - int m; - if (a > 0) - { - m = 1; - t = a; - db = da; - } - else - { - m = 0; - t = -a; - db = -da; - } - u.x = big + t; - y = t - (u.x - big); - res = do_sin (u, y, db, &cor); - cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps; - retval = ((res == res + cor) ? ((m) ? res : -res) - : bsloww1 (a, da, x, n)); - } - break; + double retval; - case 1: - case 3: - if (a < 0) - { - a = -a; - da = -da; - } - u.x = big + a; - double y = a - (u.x - big) + da; - res = do_cos (u, y, &cor); - cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps; - 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, y, s, c, sn, ssn, cs, ccs; + 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]; @@ -482,77 +218,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) { - u.x = (m > 0) ? big + x : big - x; - y = (m > 0) ? x - (u.x - big) : x + (u.x - big); - xx = y * y; - s = y + y * xx * (sn3 + xx * sn5); - c = xx * (cs2 + xx * (cs4 + xx * cs6)); - SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); - if (m <= 0) - { - sn = -sn; - ssn = -ssn; - } - cor = (ssn + s * ccs - sn * c) + cs * s; - res = sn + cor; - cor = (sn - res) + cor; - retval = (res == res + 1.096 * cor) ? res : slow1 (x); + /* Max ULP is 0.548. */ + retval = do_sin (x, 0); } /* else if (k < 0x3feb6000) */ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ else if (k < 0x400368fd) { - - y = (m > 0) ? hp0 - x : hp0 + x; - if (y >= 0) - { - u.x = big + y; - y = (y - (u.x - big)) + hp1; - } - else - { - u.x = big - y; - y = (-hp1) - (y + (u.x - big)); - } - res = do_cos (u, y, &cor); - retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x); + t = hp0 - fabs (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, 0); + 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, 0); - } /* else if (k < 0x42F00000 ) */ - -/* -----------------281474976710656 <|x| <2^1024----------------------------*/ +/* --------------------105414350 <|x| <2^1024------------------------------*/ else if (k < 0x7ff00000) - retval = reduce_and_compute (x, 0); - + { + n = __branred (x, &a, &da); + retval = do_sincos (a, da, n); + } /*--------------------- |x| > 2^1024 ----------------------------------*/ else { @@ -560,7 +253,6 @@ __sin (double x) __set_errno (EDOM); retval = x / x; } -#endif return retval; } @@ -571,23 +263,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]; @@ -599,11 +285,8 @@ __cos (double x) else if (k < 0x3feb6000) { /* 2^-27 < |x| < 0.855469 */ - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_cos (u, y, &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) @@ -611,55 +294,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 = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31; - retval = (res == res + cor) ? res : sloww (a, da, x, 1); - } - else - { - if (a > 0) - { - m = 1; - } - else - { - m = 0; - a = -a; - da = -da; - } - u.x = big + a; - y = a - (u.x - big); - res = do_sin (u, y, da, &cor); - cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31; - retval = ((res == res + cor) ? ((m) ? res : -res) - : sloww1 (a, da, x, m, 1)); - } - + /* 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, 1); + 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, 1); - } /* else if (k < 0x42F00000 ) */ - - /* 281474976710656 <|x| <2^1024 */ + /* 105414350 <|x| <2^1024 */ else if (k < 0x7ff00000) - retval = reduce_and_compute (x, 1); + { + n = __branred (x, &a, &da); + retval = do_sincos (a, da, n + 1); + } else { @@ -667,366 +318,15 @@ __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 double -SECTION -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 (x > 0) ? w[0] : -w[0]; - - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); -} - -/*******************************************************************************/ -/* 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 double -SECTION -slow1 (double x) -{ - mynumber u; - double w[2], y, cor, res; - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_sin_slow (u, y, 0, 0, &cor); - if (res == res + cor) - return (x > 0) ? res : -res; - - __dubsin (fabs (x), 0, w); - if (w[0] == w[0] + 1.000000005 * w[1]) - return (x > 0) ? w[0] : -w[0]; - - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-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 double -SECTION -slow2 (double x) -{ - mynumber u; - double w[2], y, y1, y2, cor, res, del; - - y = fabs (x); - y = hp0 - y; - if (y >= 0) - { - u.x = big + y; - y = y - (u.x - big); - del = hp1; - } - else - { - u.x = big - y; - y = -(y + (u.x - big)); - del = -hp1; - } - res = do_cos_slow (u, y, del, 0, &cor); - if (res == res + cor) - return (x > 0) ? res : -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 (x > 0) ? w[0] : -w[0]; - - return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); -} - -/***************************************************************************/ -/* Routine compute sin(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 mpsin1 or dubsin */ -/***************************************************************************/ - -static double -SECTION -sloww (double x, double dx, double orig, int k) -{ - double y, t, res, cor, w[2], a, da, xn; - mynumber v; - int4 n; - res = TAYLOR_SLOW (x, dx, cor); - - if (cor > 0) - cor = 1.0005 * cor + fabs (orig) * 3.1e-30; - else - cor = 1.0005 * cor - fabs (orig) * 3.1e-30; - - if (res == res + cor) - return res; - - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-30; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-30; - - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - - t = (orig * hpinv + toint); - xn = t - toint; - v.x = t; - y = (orig - xn * mp1) - xn * mp2; - n = (v.i[LOW_HALF] + k) & 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 || n == 1) - { - a = -a; - da = -da; - } - (a > 0) ? __dubsin (a, da, w) : __dubsin (-a, -da, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + fabs (orig) * 1.1e-40; - else - cor = 1.000000001 * w[1] - fabs (orig) * 1.1e-40; - - if (w[0] == w[0] + cor) - return (a > 0) ? w[0] : -w[0]; - - return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); -} - -/***************************************************************************/ -/* Routine compute sin(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 mpsin1 or dubsin */ -/***************************************************************************/ - -static double -SECTION -sloww1 (double x, double dx, double orig, int m, int k) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); - - if (res == res + cor) - return (m > 0) ? res : -res; - - __dubsin (x, dx, w); - - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - - if (w[0] == w[0] + cor) - return (m > 0) ? w[0] : -w[0]; - - return (k == 1) ? __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 double -SECTION -sloww2 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - u.x = big + x; - y = x - (u.x - big); - res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor); - - if (res == res + cor) - return (n & 2) ? -res : res; - - __docos (x, dx, w); - - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-30 * fabs (orig); - else - cor = 1.000000005 * w[1] - 1.1e-30 * fabs (orig); - - 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 double -SECTION -bsloww (double x, double dx, double orig, int n) -{ - double res, cor, w[2]; - - res = TAYLOR_SLOW (x, dx, cor); - cor = (cor > 0) ? 1.0005 * cor + 1.1e-24 : 1.0005 * cor - 1.1e-24; - if (res == res + cor) - return res; - - (x > 0) ? __dubsin (x, dx, w) : __dubsin (-x, -dx, w); - if (w[1] > 0) - cor = 1.000000001 * w[1] + 1.1e-24; - else - cor = 1.000000001 * w[1] - 1.1e-24; - - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - - 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 double -SECTION -bsloww1 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - dx = (x > 0) ? dx : -dx; - res = do_sin_slow (u, y, dx, 1.1e-24, &cor); - if (res == res + cor) - return (x > 0) ? res : -res; - - __dubsin (fabs (x), dx, w); - - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; - - if (w[0] == w[0] + cor) - return (x > 0) ? w[0] : -w[0]; - - 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 double -SECTION -bsloww2 (double x, double dx, double orig, int n) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - dx = (x > 0) ? dx : -dx; - res = do_cos_slow (u, y, dx, 1.1e-24, &cor); - if (res == res + cor) - return (n & 2) ? -res : res; - - __docos (fabs (x), dx, w); - - if (w[1] > 0) - cor = 1.000000005 * w[1] + 1.1e-24; - else - cor = 1.000000005 * w[1] - 1.1e-24; - - 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 double -SECTION -cslow2 (double x) -{ - mynumber u; - double w[2], y, cor, res; - - y = fabs (x); - u.x = big + y; - y = y - (u.x - big); - res = do_cos_slow (u, y, 0, 0, &cor); - if (res == res + cor) - return res; - - y = fabs (x); - __docos (y, 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 -strong_alias (__cos, __cosl) -weak_alias (__cos, cosl) -# endif +libm_alias_double (__cos, cos) #endif #ifndef __sin -weak_alias (__sin, sin) -# ifdef NO_LONG_DOUBLE -strong_alias (__sin, __sinl) -weak_alias (__sin, sinl) -# endif +libm_alias_double (__sin, sin) +#endif + #endif diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c index c389226b04..1d8d44befe 100644 --- a/sysdeps/ieee754/dbl-64/s_sincos.c +++ b/sysdeps/ieee754/dbl-64/s_sincos.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -21,43 +21,12 @@ #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.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 +36,62 @@ __sincos (double x, double *sinx, double *cosx) 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, 0); - *cosx = do_sincos_1 (a, da, x, n, 1); - - return; - } - if (k < 0x42F00000) - { - double a, da; - int4 n = reduce_sincos_2 (x, &a, &da); - - *sinx = do_sincos_2 (a, da, x, n, 0); - *cosx = do_sincos_2 (a, da, x, n, 1); - + 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; } @@ -106,8 +100,4 @@ __sincos (double x, double *sinx, double *cosx) *sinx = *cosx = x / x; } -weak_alias (__sincos, sincos) -#ifdef NO_LONG_DOUBLE -strong_alias (__sincos, __sincosl) -weak_alias (__sincos, sincosl) -#endif +libm_alias_double (__sincos, sincos) diff --git a/sysdeps/ieee754/dbl-64/s_tan.c b/sysdeps/ieee754/dbl-64/s_tan.c index ccffdd62d7..04ff8b6717 100644 --- a/sysdeps/ieee754/dbl-64/s_tan.c +++ b/sysdeps/ieee754/dbl-64/s_tan.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -41,6 +41,8 @@ #include "MathLib.h" #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> #include <fenv.h> #include <stap-probe.h> @@ -53,7 +55,7 @@ void __mptan (double, mp_no *, int); double SECTION -tan (double x) +__tan (double x) { #include "utan.h" #include "utan.tbl" @@ -843,6 +845,6 @@ tanMp (double x) return y; } -#ifdef NO_LONG_DOUBLE -weak_alias (tan, tanl) +#ifndef __tan +libm_alias_double (__tan, tan) #endif diff --git a/sysdeps/ieee754/dbl-64/s_tanh.c b/sysdeps/ieee754/dbl-64/s_tanh.c index 344a2f0330..673a97102d 100644 --- a/sysdeps/ieee754/dbl-64/s_tanh.c +++ b/sysdeps/ieee754/dbl-64/s_tanh.c @@ -41,6 +41,8 @@ static char rcsid[] = "$NetBSD: s_tanh.c,v 1.7 1995/05/10 20:48:22 jtc Exp $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-double.h> static const double one = 1.0, two = 2.0, tiny = 1.0e-300; @@ -91,8 +93,4 @@ __tanh (double x) } return (jx >= 0) ? z : -z; } -weak_alias (__tanh, tanh) -#ifdef NO_LONG_DOUBLE -strong_alias (__tanh, __tanhl) -weak_alias (__tanh, tanhl) -#endif +libm_alias_double (__tanh, tanh) diff --git a/sysdeps/ieee754/dbl-64/s_totalorder.c b/sysdeps/ieee754/dbl-64/s_totalorder.c new file mode 100644 index 0000000000..59092dceda --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_totalorder.c @@ -0,0 +1,53 @@ +/* Total order operation. dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalorder (double x, double y) +{ + int32_t hx, hy; + uint32_t lx, ly; + EXTRACT_WORDS (hx, lx, x); + EXTRACT_WORDS (hy, ly, y); +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + uint32_t uhx = hx & 0x7fffffff, uhy = hy & 0x7fffffff; + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the arguments interpreted as + sign-magnitude integers. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if ((uhx > 0x7ff00000 || (uhx == 0x7ff00000 && lx != 0)) + && (uhy > 0x7ff00000 || (uhy == 0x7ff00000 && ly != 0))) + { + hx ^= 0x00080000; + hy ^= 0x00080000; + } +#endif + uint32_t hx_sign = hx >> 31; + uint32_t hy_sign = hy >> 31; + hx ^= hx_sign >> 1; + lx ^= hx_sign; + hy ^= hy_sign >> 1; + ly ^= hy_sign; + return hx < hy || (hx == hy && lx <= ly); +} +libm_alias_double (__totalorder, totalorder) diff --git a/sysdeps/ieee754/dbl-64/s_totalordermag.c b/sysdeps/ieee754/dbl-64/s_totalordermag.c new file mode 100644 index 0000000000..3b45113962 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_totalordermag.c @@ -0,0 +1,48 @@ +/* Total order operation on absolute values. dbl-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalordermag (double x, double y) +{ + uint32_t hx, hy; + uint32_t lx, ly; + EXTRACT_WORDS (hx, lx, x); + EXTRACT_WORDS (hy, ly, y); + hx &= 0x7fffffff; + hy &= 0x7fffffff; +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the absolute values of the + arguments. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if ((hx > 0x7ff00000 || (hx == 0x7ff00000 && lx != 0)) + && (hy > 0x7ff00000 || (hy == 0x7ff00000 && ly != 0))) + { + hx ^= 0x00080000; + hy ^= 0x00080000; + } +#endif + return hx < hy || (hx == hy && lx <= ly); +} +libm_alias_double (__totalordermag, totalordermag) diff --git a/sysdeps/ieee754/dbl-64/s_trunc.c b/sysdeps/ieee754/dbl-64/s_trunc.c index f64e097f16..6ffabb410a 100644 --- a/sysdeps/ieee754/dbl-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/s_trunc.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,13 +20,14 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> double __trunc (double x) { int32_t i0, j0; - u_int32_t i1; + uint32_t i1; int sx; EXTRACT_WORDS (i0, i1, x); @@ -53,8 +54,6 @@ __trunc (double x) return x; } -weak_alias (__trunc, trunc) -#ifdef NO_LONG_DOUBLE -strong_alias (__trunc, __truncl) -weak_alias (__trunc, truncl) +#ifndef __trunc +libm_alias_double (__trunc, trunc) #endif diff --git a/sysdeps/ieee754/dbl-64/s_ufromfp.c b/sysdeps/ieee754/dbl-64/s_ufromfp.c new file mode 100644 index 0000000000..2532215981 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_ufromfp.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 0 +#define FUNC __ufromfp +#include <s_fromfp_main.c> +libm_alias_double (__ufromfp, ufromfp) diff --git a/sysdeps/ieee754/dbl-64/s_ufromfpx.c b/sysdeps/ieee754/dbl-64/s_ufromfpx.c new file mode 100644 index 0000000000..0945dfce08 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/s_ufromfpx.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 1 +#define FUNC __ufromfpx +#include <s_fromfp_main.c> +libm_alias_double (__ufromfpx, ufromfpx) diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index ad43465f64..a8572910d2 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/sincos32.h b/sysdeps/ieee754/dbl-64/sincos32.h index 5a5f2e223a..d4e72bd82a 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.h +++ b/sysdeps/ieee754/dbl-64/sincos32.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -25,7 +25,7 @@ /******************************************************************/ #ifndef SINCOS32_H -#define SINCCOS32_H +#define SINCOS32_H #ifdef BIG_ENDI static const number diff --git a/sysdeps/ieee754/dbl-64/sincostab.c b/sysdeps/ieee754/dbl-64/sincostab.c index cd88af7331..082ecad7f6 100644 --- a/sysdeps/ieee754/dbl-64/sincostab.c +++ b/sysdeps/ieee754/dbl-64/sincostab.c @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/slowexp.c b/sysdeps/ieee754/dbl-64/slowexp.c deleted file mode 100644 index b602262783..0000000000 --- a/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-2016 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 -} diff --git a/sysdeps/ieee754/dbl-64/slowpow.c b/sysdeps/ieee754/dbl-64/slowpow.c deleted file mode 100644 index 266910c979..0000000000 --- a/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-2016 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; -} diff --git a/sysdeps/ieee754/dbl-64/t_exp.c b/sysdeps/ieee754/dbl-64/t_exp.c index 3c921f1d3d..555c4ff01b 100644 --- a/sysdeps/ieee754/dbl-64/t_exp.c +++ b/sysdeps/ieee754/dbl-64/t_exp.c @@ -1,5 +1,5 @@ /* Accurate tables for exp(). - Copyright (C) 1998-2016 Free Software Foundation, Inc. + Copyright (C) 1998-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Geoffrey Keating <geoffk@ozemail.com.au> diff --git a/sysdeps/ieee754/dbl-64/uasncs.h b/sysdeps/ieee754/dbl-64/uasncs.h index e88bd1bf3e..e980fc5c55 100644 --- a/sysdeps/ieee754/dbl-64/uasncs.h +++ b/sysdeps/ieee754/dbl-64/uasncs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uatan.tbl b/sysdeps/ieee754/dbl-64/uatan.tbl index 979d65d8b2..3f0b2bff90 100644 --- a/sysdeps/ieee754/dbl-64/uatan.tbl +++ b/sysdeps/ieee754/dbl-64/uatan.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/uexp.h b/sysdeps/ieee754/dbl-64/uexp.h index eb7aa5f5fb..64ab2c8fca 100644 --- a/sysdeps/ieee754/dbl-64/uexp.h +++ b/sysdeps/ieee754/dbl-64/uexp.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -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; diff --git a/sysdeps/ieee754/dbl-64/uexp.tbl b/sysdeps/ieee754/dbl-64/uexp.tbl index 9abfc7ef71..b1eff5253e 100644 --- a/sysdeps/ieee754/dbl-64/uexp.tbl +++ b/sysdeps/ieee754/dbl-64/uexp.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/ulog.h b/sysdeps/ieee754/dbl-64/ulog.h index d507c693d8..087b76e2ab 100644 --- a/sysdeps/ieee754/dbl-64/ulog.h +++ b/sysdeps/ieee754/dbl-64/ulog.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -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 diff --git a/sysdeps/ieee754/dbl-64/ulog.tbl b/sysdeps/ieee754/dbl-64/ulog.tbl index dd92cbb886..93fa1b36d3 100644 --- a/sysdeps/ieee754/dbl-64/ulog.tbl +++ b/sysdeps/ieee754/dbl-64/ulog.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/upow.h b/sysdeps/ieee754/dbl-64/upow.h index dd9ba32b71..c1b9d8e3cc 100644 --- a/sysdeps/ieee754/dbl-64/upow.h +++ b/sysdeps/ieee754/dbl-64/upow.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/upow.tbl b/sysdeps/ieee754/dbl-64/upow.tbl index 55d878edad..c1ad333583 100644 --- a/sysdeps/ieee754/dbl-64/upow.tbl +++ b/sysdeps/ieee754/dbl-64/upow.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/urem.h b/sysdeps/ieee754/dbl-64/urem.h index 42b26273fd..c810791a1c 100644 --- a/sysdeps/ieee754/dbl-64/urem.h +++ b/sysdeps/ieee754/dbl-64/urem.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/usncs.h b/sysdeps/ieee754/dbl-64/usncs.h index 1408695c85..97d1418a23 100644 --- a/sysdeps/ieee754/dbl-64/usncs.h +++ b/sysdeps/ieee754/dbl-64/usncs.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/utan.h b/sysdeps/ieee754/dbl-64/utan.h index 0a99812754..d619d0004b 100644 --- a/sysdeps/ieee754/dbl-64/utan.h +++ b/sysdeps/ieee754/dbl-64/utan.h @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/utan.tbl b/sysdeps/ieee754/dbl-64/utan.tbl index 8a8c2afa53..92163ec3e8 100644 --- a/sysdeps/ieee754/dbl-64/utan.tbl +++ b/sysdeps/ieee754/dbl-64/utan.tbl @@ -1,7 +1,7 @@ /* * IBM Accurate Mathematical Library * Written by International Business Machines Corp. - * Copyright (C) 2001-2016 Free Software Foundation, Inc. + * Copyright (C) 2001-2018 Free Software Foundation, Inc. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c index ccccdaf106..0af05a0222 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_acosh.c @@ -51,13 +51,13 @@ __ieee754_acosh (double x) /* 2**28 > x > 2 */ double t = x * x; - return __ieee754_log (2.0 * x - one / (x + __ieee754_sqrt (t - one))); + return __ieee754_log (2.0 * x - one / (x + sqrt (t - one))); } else if (__glibc_likely (hx > INT64_C (0x3ff0000000000000))) { /* 1<x<2 */ double t = x - one; - return __log1p (t + __ieee754_sqrt (2.0 * t + t * t)); + return __log1p (t + sqrt (2.0 * t + t * t)); } else if (__glibc_likely (hx == INT64_C (0x3ff0000000000000))) return 0.0; /* acosh(1) = 0 */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c index 4f5a81669e..cd5567182f 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_log10.c @@ -65,7 +65,7 @@ __ieee754_log10 (double x) if (hx < INT64_C(0x0010000000000000)) { /* x < 2**-1022 */ if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0)) - return -two54 / (x - x); /* log(+-0)=-inf */ + return -two54 / fabs (x); /* log(+-0)=-inf */ if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c b/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c index 5ccb78cf03..f08d5b337d 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/e_log2.c @@ -81,7 +81,7 @@ __ieee754_log2 (double x) if (hx < INT64_C(0x0010000000000000)) { /* x < 2**-1022 */ if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0)) - return -two54 / (x - x); /* log(+-0)=-inf */ + return -two54 / fabs (x); /* log(+-0)=-inf */ if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); /* log(-#) = NaN */ k -= 54; diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c index c687525d4b..b99829d2b0 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_ceil.c @@ -15,14 +15,11 @@ * Return x rounded toward -inf to integral value * Method: * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to ceil(x). */ #include <math.h> #include <math_private.h> - -static const double huge = 1.0e300; +#include <libm-alias-double.h> double __ceil(double x) @@ -32,14 +29,13 @@ __ceil(double x) EXTRACT_WORDS64(i0,x); j0 = ((i0>>52)&0x7ff)-0x3ff; if(j0<=51) { - if(j0<0) { /* raise inexact if x != 0 */ - math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ if(i0<0) {i0=INT64_C(0x8000000000000000);} else if(i0!=0) { i0=INT64_C(0x3ff0000000000000);} } else { i = INT64_C(0x000fffffffffffff)>>j0; if((i0&i)==0) return x; /* x is integral */ - math_force_eval(huge+x); /* raise inexact flag */ if(i0>0) i0 += UINT64_C(0x0010000000000000)>>j0; i0 &= (~i); } @@ -51,9 +47,5 @@ __ceil(double x) return x; } #ifndef __ceil -weak_alias (__ceil, ceil) -# ifdef NO_LONG_DOUBLE -strong_alias (__ceil, __ceill) -weak_alias (__ceil, ceill) -# endif +libm_alias_double (__ceil, ceil) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c index ef51608f6e..40676924fe 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_finite.c @@ -16,6 +16,7 @@ #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> #include <stdint.h> @@ -30,7 +31,7 @@ __finite(double x) hidden_def (__finite) weak_alias (__finite, finite) #ifdef NO_LONG_DOUBLE -# ifdef LDBL_CLASSIFY_COMPAT +# if LDBL_CLASSIFY_COMPAT # if SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __finite, __finitel, GLIBC_2_0); # endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c index b7ed14bfa2..f7e0a77ec3 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_floor.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 2011-2016 Free Software Foundation, Inc. + Copyright (C) 2011-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 2011. @@ -33,18 +33,15 @@ #include <math.h> #include <math_private.h> #include <stdint.h> +#include <libm-alias-double.h> /* * floor(x) * Return x rounded toward -inf to integral value * Method: * Bit twiddling. - * Exception: - * Inexact flag raised if x not equal to floor(x). */ -static const double huge = 1.0e300; - double __floor (double x) @@ -53,15 +50,14 @@ __floor (double x) EXTRACT_WORDS64(i0,x); int32_t j0 = ((i0>>52)&0x7ff)-0x3ff; if(__builtin_expect(j0<52, 1)) { - if(j0<0) { /* raise inexact if x != 0 */ - math_force_eval(huge+x);/* return 0*sign(x) if |x|<1 */ + if(j0<0) { + /* return 0*sign(x) if |x|<1 */ if(i0>=0) {i0=0;} else if((i0&0x7fffffffffffffffl)!=0) { i0=0xbff0000000000000l;} } else { uint64_t i = (0x000fffffffffffffl)>>j0; if((i0&i)==0) return x; /* x is integral */ - math_force_eval(huge+x); /* raise inexact flag */ if(i0<0) i0 += (0x0010000000000000l)>>j0; i0 &= (~i); } @@ -71,9 +67,5 @@ __floor (double x) return x; } #ifndef __floor -weak_alias (__floor, floor) -# ifdef NO_LONG_DOUBLE -strong_alias (__floor, __floorl) -weak_alias (__floor, floorl) -# endif +libm_alias_double (__floor, floor) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c index 42068f8187..c73434f5f3 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_frexp.c @@ -1,4 +1,4 @@ -/* Copyright (C) 2011-2016 Free Software Foundation, Inc. +/* Copyright (C) 2011-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -19,6 +19,7 @@ #include <inttypes.h> #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> /* * for non-zero, finite x @@ -55,12 +56,11 @@ __frexp (double x, int *eptr) ix = (ix & INT64_C (0x800fffffffffffff)) | INT64_C (0x3fe0000000000000); INSERT_WORDS64 (x, ix); } + else + /* Quiet signaling NaNs. */ + x += x; *eptr = e; return x; } -weak_alias (__frexp, frexp) -#ifdef NO_LONG_DOUBLE -strong_alias (__frexp, __frexpl) -weak_alias (__frexp, frexpl) -#endif +libm_alias_double (__frexp, frexp) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c new file mode 100644 index 0000000000..1103ef2ed9 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c @@ -0,0 +1,32 @@ +/* Get NaN payload. dbl-64/wordsize-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <stdint.h> + +double +__getpayload (const double *x) +{ + uint64_t ix; + EXTRACT_WORDS64 (ix, *x); + ix &= 0x7ffffffffffffULL; + return (double) ix; +} +libm_alias_double (__getpayload, getpayload) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c index 951fb73239..2b427a8b4c 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isinf.c @@ -11,6 +11,7 @@ #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> int @@ -26,7 +27,7 @@ __isinf (double x) hidden_def (__isinf) weak_alias (__isinf, isinf) #ifdef NO_LONG_DOUBLE -# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +# if LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __isinf, __isinfl, GLIBC_2_0); # endif weak_alias (__isinf, isinfl) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c index bcff9e3b67..cd805d157b 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_isnan.c @@ -17,6 +17,7 @@ #include <math.h> #include <math_private.h> +#include <ldbl-classify-compat.h> #include <shlib-compat.h> #include <stdint.h> @@ -32,7 +33,7 @@ int __isnan(double x) hidden_def (__isnan) weak_alias (__isnan, isnan) #ifdef NO_LONG_DOUBLE -# if defined LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) +# if LDBL_CLASSIFY_COMPAT && SHLIB_COMPAT (libc, GLIBC_2_0, GLIBC_2_23) compat_symbol (libc, __isnan, __isnanl, GLIBC_2_0); # endif weak_alias (__isnan, isnanl) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c index c22e608c6e..63ef5ca7d8 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_issignaling.c @@ -1,5 +1,5 @@ /* Test for signaling NaN. - Copyright (C) 2013-2016 Free Software Foundation, Inc. + Copyright (C) 2013-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,13 +18,14 @@ #include <math.h> #include <math_private.h> +#include <nan-high-order-bit.h> int __issignaling (double x) { - u_int64_t xi; + uint64_t xi; EXTRACT_WORDS64 (xi, x); -#ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN /* We only have to care about the high-order bit of x's significand, because having it set (sNaN) already makes the significand different from that used to designate infinity. */ diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c index 73ddd8dd9f..e3b18eb893 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_llround.c @@ -1,5 +1,5 @@ /* Round double value to long long int. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -24,6 +24,7 @@ #include <sysdep.h> #include <math_private.h> +#include <libm-alias-double.h> long long int @@ -63,20 +64,12 @@ __llround (double x) return sign * result; } -weak_alias (__llround, llround) -#ifdef NO_LONG_DOUBLE -strong_alias (__llround, __llroundl) -weak_alias (__llround, llroundl) -#endif +libm_alias_double (__llround, llround) /* long has the same width as long long on LP64 machines, so use an alias. */ #undef lround #undef __lround #ifdef _LP64 strong_alias (__llround, __lround) -weak_alias (__llround, lround) -# ifdef NO_LONG_DOUBLE -strong_alias (__llround, __lroundl) -weak_alias (__llround, lroundl) -# endif +libm_alias_double (__lround, lround) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c index c7846054a6..c7fa8169c5 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_logb.c @@ -1,5 +1,5 @@ /* Compute radix independent exponent. - Copyright (C) 2011-2016 Free Software Foundation, Inc. + Copyright (C) 2011-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@gmail.com>, 2011. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> double @@ -41,8 +42,6 @@ __logb (double x) } return (double) (ex - 1023); } -weak_alias (__logb, logb) -#ifdef NO_LONG_DOUBLE -strong_alias (__logb, __logbl) -weak_alias (__logb, logbl) +#ifndef __logb +libm_alias_double (__logb, logb) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c index 9e665b01ec..a88c6c8788 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_lround.c @@ -1,5 +1,5 @@ /* Round double value to long int. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -21,6 +21,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> /* For LP64, lround is an alias for llround. */ #ifndef _LP64 @@ -80,10 +81,6 @@ __lround (double x) return sign * result; } -weak_alias (__lround, lround) -# ifdef NO_LONG_DOUBLE -strong_alias (__lround, __lroundl) -weak_alias (__lround, lroundl) -# endif +libm_alias_double (__lround, lround) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c index c309e56272..8d14e78ef0 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_modf.c @@ -22,6 +22,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <stdint.h> static const double one = 1.0; @@ -59,8 +60,6 @@ __modf(double x, double *iptr) return x; } } -weak_alias (__modf, modf) -#ifdef NO_LONG_DOUBLE -strong_alias (__modf, __modfl) -weak_alias (__modf, modfl) +#ifndef __modf +libm_alias_double (__modf, modf) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c index 8293819981..a4a081724e 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c @@ -22,7 +22,9 @@ #include <fenv.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <libm-alias-double.h> static const double TWO52[2]={ @@ -42,9 +44,9 @@ __nearbyint(double x) if(__builtin_expect(j0<52, 1)) { if(j0<0) { libc_feholdexcept (&env); - double w = TWO52[sx]+x; + double w = TWO52[sx] + math_opt_barrier (x); double t = w-TWO52[sx]; - math_opt_barrier(t); + math_force_eval (t); libc_fesetenv (&env); return __copysign (t, x); } @@ -53,14 +55,10 @@ __nearbyint(double x) else return x; /* x is integral */ } libc_feholdexcept (&env); - double w = TWO52[sx]+x; + double w = TWO52[sx] + math_opt_barrier (x); double t = w-TWO52[sx]; - math_opt_barrier (t); + math_force_eval (t); libc_fesetenv (&env); return t; } -weak_alias (__nearbyint, nearbyint) -#ifdef NO_LONG_DOUBLE -strong_alias (__nearbyint, __nearbyintl) -weak_alias (__nearbyint, nearbyintl) -#endif +libm_alias_double (__nearbyint, nearbyint) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c index 716e07e54d..bbfb1195ac 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_remquo.c @@ -1,5 +1,5 @@ /* Compute remainder and a congruent to the quotient. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <stdint.h> static const double zero = 0.0; @@ -107,8 +108,4 @@ __remquo (double x, double y, int *quo) x = -x; return x; } -weak_alias (__remquo, remquo) -#ifdef NO_LONG_DOUBLE -strong_alias (__remquo, __remquol) -weak_alias (__remquo, remquol) -#endif +libm_alias_double (__remquo, remquo) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c index 87b2339d43..622e479c5f 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c @@ -21,6 +21,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> static const double TWO52[2]={ @@ -52,9 +53,5 @@ __rint(double x) return w-TWO52[sx]; } #ifndef __rint -weak_alias (__rint, rint) -# ifdef NO_LONG_DOUBLE -strong_alias (__rint, __rintl) -weak_alias (__rint, rintl) -# endif +libm_alias_double (__rint, rint) #endif diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c index f1d75bbe68..3323621ce3 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_round.c @@ -1,5 +1,5 @@ /* Round double to integer away from zero. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,10 +20,9 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> #include <stdint.h> -static const double huge = 1.0e300; - double __round (double x) @@ -36,8 +35,6 @@ __round (double x) { if (j0 < 0) { - math_force_eval (huge + x); - i0 &= UINT64_C(0x8000000000000000); if (j0 == -1) i0 |= UINT64_C(0x3ff0000000000000); @@ -48,9 +45,7 @@ __round (double x) if ((i0 & i) == 0) /* X is integral. */ return x; - math_force_eval (huge + x); - /* Raise inexact if x != 0. */ i0 += UINT64_C(0x0008000000000000) >> j0; i0 &= ~i; } @@ -67,8 +62,4 @@ __round (double x) INSERT_WORDS64 (x, i0); return x; } -weak_alias (__round, round) -#ifdef NO_LONG_DOUBLE -strong_alias (__round, __roundl) -weak_alias (__round, roundl) -#endif +libm_alias_double (__round, round) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c new file mode 100644 index 0000000000..7bbbb2dc20 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c @@ -0,0 +1,71 @@ +/* Round to nearest integer value, rounding halfway cases to even. + dbl-64/wordsize-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <stdint.h> + +#define BIAS 0x3ff +#define MANT_DIG 53 +#define MAX_EXP (2 * BIAS + 1) + +double +__roundeven (double x) +{ + uint64_t ix, ux; + EXTRACT_WORDS64 (ix, x); + ux = ix & 0x7fffffffffffffffULL; + int exponent = ux >> (MANT_DIG - 1); + if (exponent >= BIAS + MANT_DIG - 1) + { + /* Integer, infinity or NaN. */ + if (exponent == MAX_EXP) + /* Infinity or NaN; quiet signaling NaNs. */ + return x + x; + else + return x; + } + else if (exponent >= BIAS) + { + /* At least 1; not necessarily an integer. Locate the bits with + exponents 0 and -1 (when the unbiased exponent is 0, the bit + with exponent 0 is implicit, but as the bias is odd it is OK + to take it from the low bit of the exponent). */ + int int_pos = (BIAS + MANT_DIG - 1) - exponent; + int half_pos = int_pos - 1; + uint64_t half_bit = 1ULL << half_pos; + uint64_t int_bit = 1ULL << int_pos; + if ((ix & (int_bit | (half_bit - 1))) != 0) + /* Carry into the exponent works correctly. No need to test + whether HALF_BIT is set. */ + ix += half_bit; + ix &= ~(int_bit - 1); + } + else if (exponent == BIAS - 1 && ux > 0x3fe0000000000000ULL) + /* Interval (0.5, 1). */ + ix = (ix & 0x8000000000000000ULL) | 0x3ff0000000000000ULL; + else + /* Rounds to 0. */ + ix &= 0x8000000000000000ULL; + INSERT_WORDS64 (x, ix); + return x; +} +hidden_def (__roundeven) +libm_alias_double (__roundeven, roundeven) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c new file mode 100644 index 0000000000..6af92d7fe2 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c @@ -0,0 +1,54 @@ +/* Set NaN payload. dbl-64/wordsize-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <libm-alias-double.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +#define SET_HIGH_BIT (HIGH_ORDER_BIT_IS_SET_FOR_SNAN ? SIG : !SIG) +#define BIAS 0x3ff +#define PAYLOAD_DIG 51 +#define EXPLICIT_MANT_DIG 52 + +int +FUNC (double *x, double payload) +{ + uint64_t ix; + EXTRACT_WORDS64 (ix, payload); + int exponent = ix >> EXPLICIT_MANT_DIG; + /* Test if argument is (a) negative or too large; (b) too small, + except for 0 when allowed; (c) not an integer. */ + if (exponent >= BIAS + PAYLOAD_DIG + || (exponent < BIAS && !(SET_HIGH_BIT && ix == 0)) + || (ix & ((1ULL << (BIAS + EXPLICIT_MANT_DIG - exponent)) - 1)) != 0) + { + INSERT_WORDS64 (*x, 0); + return 1; + } + if (ix != 0) + { + ix &= (1ULL << EXPLICIT_MANT_DIG) - 1; + ix |= 1ULL << EXPLICIT_MANT_DIG; + ix >>= BIAS + EXPLICIT_MANT_DIG - exponent; + } + ix |= 0x7ff0000000000000ULL | (SET_HIGH_BIT ? 0x8000000000000ULL : 0); + INSERT_WORDS64 (*x, ix); + return 0; +} diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c new file mode 100644 index 0000000000..c235e55c20 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c @@ -0,0 +1,49 @@ +/* Total order operation. dbl-64/wordsize-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <nan-high-order-bit.h> +#include <libm-alias-double.h> +#include <stdint.h> + +int +__totalorder (double x, double y) +{ + int64_t ix, iy; + EXTRACT_WORDS64 (ix, x); + EXTRACT_WORDS64 (iy, y); +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the arguments interpreted as + sign-magnitude integers. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if ((ix & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL + && (iy & 0x7fffffffffffffffULL) > 0x7ff0000000000000ULL) + { + ix ^= 0x0008000000000000ULL; + iy ^= 0x0008000000000000ULL; + } +#endif + uint64_t ix_sign = ix >> 63; + uint64_t iy_sign = iy >> 63; + ix ^= ix_sign >> 1; + iy ^= iy_sign >> 1; + return ix <= iy; +} +libm_alias_double (__totalorder, totalorder) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c new file mode 100644 index 0000000000..6da4f118a3 --- /dev/null +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c @@ -0,0 +1,46 @@ +/* Total order operation on absolute values. dbl-64/wordsize-64 version. + Copyright (C) 2016-2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <math.h> +#include <math_private.h> +#include <nan-high-order-bit.h> +#include <libm-alias-double.h> +#include <stdint.h> + +int +__totalordermag (double x, double y) +{ + uint64_t ix, iy; + EXTRACT_WORDS64 (ix, x); + EXTRACT_WORDS64 (iy, y); + ix &= 0x7fffffffffffffffULL; + iy &= 0x7fffffffffffffffULL; +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN + /* For the preferred quiet NaN convention, this operation is a + comparison of the representations of the absolute values of the + arguments. If both arguments are NaNs, invert the + quiet/signaling bit so comparing that way works. */ + if (ix > 0x7ff0000000000000ULL && iy > 0x7ff0000000000000ULL) + { + ix ^= 0x0008000000000000ULL; + iy ^= 0x0008000000000000ULL; + } +#endif + return ix <= iy; +} +libm_alias_double (__totalordermag, totalordermag) diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c index 81ac55e2f6..19a09b894e 100644 --- a/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c +++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_trunc.c @@ -1,5 +1,5 @@ /* Truncate argument to nearest integral value not larger than the argument. - Copyright (C) 1997-2016 Free Software Foundation, Inc. + Copyright (C) 1997-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. @@ -20,6 +20,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-double.h> double @@ -48,8 +49,6 @@ __trunc (double x) return x; } -weak_alias (__trunc, trunc) -#ifdef NO_LONG_DOUBLE -strong_alias (__trunc, __truncl) -weak_alias (__trunc, truncl) +#ifndef __trunc +libm_alias_double (__trunc, trunc) #endif diff --git a/sysdeps/ieee754/dbl-64/x2y2m1.c b/sysdeps/ieee754/dbl-64/x2y2m1.c index 96078888a7..aa0d9703fa 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2016 Free Software Foundation, Inc. + Copyright (C) 2012-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -18,7 +18,7 @@ #include <math.h> #include <math_private.h> -#include <float.h> +#include <mul_split.h> #include <stdlib.h> /* Calculate X + Y exactly and store the result in *HI + *LO. It is @@ -33,36 +33,6 @@ add_split (double *hi, double *lo, double x, double y) *lo = (x - *hi) + y; } -/* Calculate X * Y exactly and store the result in *HI + *LO. It is - given that the values are small enough that no overflow occurs and - large enough (or zero) that no underflow occurs. */ - -static void -mul_split (double *hi, double *lo, double x, double y) -{ -#ifdef __FP_FAST_FMA - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fma (x, y, -*hi); -#elif defined FP_FAST_FMA - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fma (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1) - double x1 = x * C; - double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - double x2 = x - x1; - double y2 = y - y1; - *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; -#endif -} - /* Compare absolute values of floating-point values pointed to by P and Q for qsort. */ diff --git a/sysdeps/ieee754/dbl-64/x2y2m1f.c b/sysdeps/ieee754/dbl-64/x2y2m1f.c index 480b4c0c67..a61dc4185f 100644 --- a/sysdeps/ieee754/dbl-64/x2y2m1f.c +++ b/sysdeps/ieee754/dbl-64/x2y2m1f.c @@ -1,5 +1,5 @@ /* Compute x^2 + y^2 - 1, without large cancellation error. - Copyright (C) 2012-2016 Free Software Foundation, Inc. + Copyright (C) 2012-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or |