diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-96')
80 files changed, 1369 insertions, 262 deletions
diff --git a/sysdeps/ieee754/ldbl-96/Makefile b/sysdeps/ieee754/ldbl-96/Makefile new file mode 100644 index 0000000000..790f670e44 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/Makefile @@ -0,0 +1,21 @@ +# Makefile for sysdeps/ieee754/ldbl-96. +# 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/>. + +ifeq ($(subdir),math) +tests += test-canonical-ldbl-96 test-totalorderl-ldbl-96 +endif diff --git a/sysdeps/ieee754/ldbl-96/bits/iscanonical.h b/sysdeps/ieee754/ldbl-96/bits/iscanonical.h new file mode 100644 index 0000000000..e1ee1356b7 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/bits/iscanonical.h @@ -0,0 +1,54 @@ +/* Define iscanonical macro. ldbl-96 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/>. */ + +#ifndef _MATH_H +# error "Never use <bits/iscanonical.h> directly; include <math.h> instead." +#endif + +extern int __iscanonicall (long double __x) + __THROW __attribute__ ((__const__)); +#define __iscanonicalf(x) ((void) (__typeof (x)) (x), 1) +#define __iscanonical(x) ((void) (__typeof (x)) (x), 1) +#if __HAVE_DISTINCT_FLOAT128 +# define __iscanonicalf128(x) ((void) (__typeof (x)) (x), 1) +#endif + +/* Return nonzero value if X is canonical. In IEEE interchange binary + formats, all values are canonical, but the argument must still be + converted to its semantic type for any exceptions arising from the + conversion, before being discarded; in extended precision, there + are encodings that are not consistently handled as corresponding to + any particular value of the type, and we return 0 for those. */ +#ifndef __cplusplus +# define iscanonical(x) __MATH_TG ((x), __iscanonical, (x)) +#else +/* In C++ mode, __MATH_TG cannot be used, because it relies on + __builtin_types_compatible_p, which is a C-only builtin. On the + other hand, overloading provides the means to distinguish between + the floating-point types. The overloading resolution will match + the correct parameter (regardless of type qualifiers (i.e.: const + and volatile)). */ +extern "C++" { +inline int iscanonical (float __val) { return __iscanonicalf (__val); } +inline int iscanonical (double __val) { return __iscanonical (__val); } +inline int iscanonical (long double __val) { return __iscanonicall (__val); } +# if __HAVE_DISTINCT_FLOAT128 +inline int iscanonical (_Float128 __val) { return __iscanonicalf128 (__val); } +# endif +} +#endif /* __cplusplus */ diff --git a/sysdeps/ieee754/ldbl-96/bits/long-double.h b/sysdeps/ieee754/ldbl-96/bits/long-double.h new file mode 100644 index 0000000000..28488e0b05 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/bits/long-double.h @@ -0,0 +1,20 @@ +/* Properties of long double type. ldbl-96 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 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/>. */ + +/* long double is distinct from double, so there is nothing to + define here. */ diff --git a/sysdeps/ieee754/ldbl-96/e_acoshl.c b/sysdeps/ieee754/ldbl-96/e_acoshl.c index cf9a6db0ef..56b04d4cc4 100644 --- a/sysdeps/ieee754/ldbl-96/e_acoshl.c +++ b/sysdeps/ieee754/ldbl-96/e_acoshl.c @@ -39,7 +39,7 @@ long double __ieee754_acoshl(long double x) { long double t; - u_int32_t se,i0,i1; + uint32_t se,i0,i1; GET_LDOUBLE_WORDS(se,i0,i1,x); if(se<0x3fff || se & 0x8000) { /* x < 1 */ return (x-x)/(x-x); @@ -52,10 +52,10 @@ __ieee754_acoshl(long double x) return 0.0; /* acosh(1) = 0 */ } else if (se > 0x4000) { /* 2**28 > x > 2 */ t=x*x; - return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one))); + return __ieee754_logl(2.0*x-one/(x+sqrtl(t-one))); } else { /* 1<x<2 */ t = x-one; - return __log1pl(t+__ieee754_sqrtl(2.0*t+t*t)); + return __log1pl(t+sqrtl(2.0*t+t*t)); } } strong_alias (__ieee754_acoshl, __acoshl_finite) diff --git a/sysdeps/ieee754/ldbl-96/e_asinl.c b/sysdeps/ieee754/ldbl-96/e_asinl.c index f52b931459..806906a58a 100644 --- a/sysdeps/ieee754/ldbl-96/e_asinl.c +++ b/sysdeps/ieee754/ldbl-96/e_asinl.c @@ -61,6 +61,7 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const long double one = 1.0L, @@ -96,7 +97,7 @@ __ieee754_asinl (long double x) { long double t, w, p, q, c, r, s; int32_t ix; - u_int32_t se, i0, i1, k; + uint32_t se, i0, i1, k; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -132,7 +133,7 @@ __ieee754_asinl (long double x) t = w * 0.5; p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5))))); q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t)))); - s = __ieee754_sqrtl (t); + s = sqrtl (t); if (ix >= 0x3ffef999) { /* if |x| > 0.975 */ w = p / q; diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c index b99a83c6ee..7312f84329 100644 --- a/sysdeps/ieee754/ldbl-96/e_atanhl.c +++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c @@ -34,7 +34,9 @@ #include <float.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> +#include <math-underflow.h> static const long double one = 1.0, huge = 1e4900L; @@ -45,7 +47,7 @@ __ieee754_atanhl(long double x) { long double t; int32_t ix; - u_int32_t se,i0,i1; + uint32_t se,i0,i1; GET_LDOUBLE_WORDS(se,i0,i1,x); ix = se&0x7fff; if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff) diff --git a/sysdeps/ieee754/ldbl-96/e_coshl.c b/sysdeps/ieee754/ldbl-96/e_coshl.c index dd22cae363..1edf2c1542 100644 --- a/sysdeps/ieee754/ldbl-96/e_coshl.c +++ b/sysdeps/ieee754/ldbl-96/e_coshl.c @@ -44,7 +44,7 @@ __ieee754_coshl (long double x) { long double t,w; int32_t ex; - u_int32_t mx,lx; + uint32_t mx,lx; /* High word of |x|. */ GET_LDOUBLE_WORDS(ex,mx,lx,x); diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c index 8dd7a03918..fc7a5c55dc 100644 --- a/sysdeps/ieee754/ldbl-96/e_gammal_r.c +++ b/sysdeps/ieee754/ldbl-96/e_gammal_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. @@ -19,6 +19,7 @@ #include <math.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 @@ -100,7 +101,7 @@ gammal_positive (long double x, int *exp2_adj) long double ret = (__ieee754_powl (x_adj_mant, x_adj) * __ieee754_exp2l (x_adj_log2 * x_adj_frac) * __ieee754_expl (-x_adj) - * __ieee754_sqrtl (2 * M_PIl / x_adj) + * sqrtl (2 * M_PIl / x_adj) / prod); exp_adj += x_eps * __ieee754_logl (x_adj); long double bsum = gamma_coeff[NCOEFF - 1]; @@ -115,7 +116,7 @@ gammal_positive (long double x, int *exp2_adj) long double __ieee754_gammal_r (long double x, int *signgamp) { - u_int32_t es, hx, lx; + uint32_t es, hx, lx; long double ret; GET_LDOUBLE_WORDS (es, hx, lx, x); diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c index ee3a07055b..f664e30c98 100644 --- a/sysdeps/ieee754/ldbl-96/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c @@ -48,11 +48,12 @@ #include <math.h> #include <math_private.h> +#include <math-underflow.h> long double __ieee754_hypotl(long double x, long double y) { long double a,b,t1,t2,y1,y2,w; - u_int32_t j,k,ea,eb; + uint32_t j,k,ea,eb; GET_LDOUBLE_EXP(ea,x); ea &= 0x7fff; @@ -65,9 +66,11 @@ long double __ieee754_hypotl(long double x, long double y) k=0; if(__builtin_expect(ea > 0x5f3f,0)) { /* a>2**8000 */ if(ea == 0x7fff) { /* Inf or NaN */ - u_int32_t exp __attribute__ ((unused)); - u_int32_t high,low; + uint32_t exp __attribute__ ((unused)); + uint32_t high,low; w = a+b; /* for sNaN */ + if (issignaling (a) || issignaling (b)) + return w; GET_LDOUBLE_WORDS(exp,high,low,a); if(((high&0x7fffffff)|low)==0) w = a; GET_LDOUBLE_WORDS(exp,high,low,b); @@ -81,8 +84,8 @@ long double __ieee754_hypotl(long double x, long double y) } if(__builtin_expect(eb < 0x20bf, 0)) { /* b < 2**-8000 */ if(eb == 0) { /* subnormal b or 0 */ - u_int32_t exp __attribute__ ((unused)); - u_int32_t high,low; + uint32_t exp __attribute__ ((unused)); + uint32_t high,low; GET_LDOUBLE_WORDS(exp,high,low,b); if((high|low)==0) return a; SET_LDOUBLE_WORDS(t1, 0x7ffd, 0x80000000, 0); /* t1=2^16382 */ @@ -111,13 +114,13 @@ long double __ieee754_hypotl(long double x, long double y) /* medium size a and b */ w = a-b; if (w>b) { - u_int32_t high; + uint32_t high; GET_LDOUBLE_MSW(high,a); SET_LDOUBLE_WORDS(t1,ea,high,0); t2 = a-t1; - w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); + w = sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); } else { - u_int32_t high; + uint32_t high; GET_LDOUBLE_MSW(high,b); a = a+a; SET_LDOUBLE_WORDS(y1,eb,high,0); @@ -125,10 +128,10 @@ long double __ieee754_hypotl(long double x, long double y) GET_LDOUBLE_MSW(high,a); SET_LDOUBLE_WORDS(t1,ea+1,high,0); t2 = a - t1; - w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); + w = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { - u_int32_t exp; + uint32_t exp; t1 = 1.0; GET_LDOUBLE_EXP(exp,t1); SET_LDOUBLE_EXP(t1,exp+k); diff --git a/sysdeps/ieee754/ldbl-96/e_j0l.c b/sysdeps/ieee754/ldbl-96/e_j0l.c index a536054cde..e720ae9558 100644 --- a/sysdeps/ieee754/ldbl-96/e_j0l.c +++ b/sysdeps/ieee754/ldbl-96/e_j0l.c @@ -72,6 +72,7 @@ */ #include <math.h> +#include <math-barriers.h> #include <math_private.h> static long double pzero (long double), qzero (long double); @@ -108,7 +109,7 @@ __ieee754_j0l (long double x) { long double z, s, c, ss, cc, r, u, v; int32_t ix; - u_int32_t se; + uint32_t se; GET_LDOUBLE_EXP (se, x); ix = se & 0x7fff; @@ -133,12 +134,12 @@ __ieee754_j0l (long double x) * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) */ if (__glibc_unlikely (ix > 0x4080)) /* 2^129 */ - z = (invsqrtpi * cc) / __ieee754_sqrtl (x); + z = (invsqrtpi * cc) / sqrtl (x); else { u = pzero (x); v = qzero (x); - z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x); + z = invsqrtpi * (u * cc - v * ss) / sqrtl (x); } return z; } @@ -194,7 +195,7 @@ __ieee754_y0l (long double x) { long double z, s, c, ss, cc, u, v; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -235,12 +236,12 @@ __ieee754_y0l (long double x) ss = z / cc; } if (__glibc_unlikely (ix > 0x4080)) /* 1e39 */ - z = (invsqrtpi * ss) / __ieee754_sqrtl (x); + z = (invsqrtpi * ss) / sqrtl (x); else { u = pzero (x); v = qzero (x); - z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x); + z = invsqrtpi * (u * ss + v * cc) / sqrtl (x); } return z; } @@ -352,7 +353,7 @@ pzero (long double x) const long double *p, *q; long double z, r, s; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -490,7 +491,7 @@ qzero (long double x) const long double *p, *q; long double s, r, z; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; diff --git a/sysdeps/ieee754/ldbl-96/e_j1l.c b/sysdeps/ieee754/ldbl-96/e_j1l.c index e8a7349cf4..581615d563 100644 --- a/sysdeps/ieee754/ldbl-96/e_j1l.c +++ b/sysdeps/ieee754/ldbl-96/e_j1l.c @@ -75,6 +75,7 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static long double pone (long double), qone (long double); @@ -112,7 +113,7 @@ __ieee754_j1l (long double x) { long double z, c, r, s, ss, cc, u, v, y; int32_t ix; - u_int32_t se; + uint32_t se; GET_LDOUBLE_EXP (se, x); ix = se & 0x7fff; @@ -137,12 +138,12 @@ __ieee754_j1l (long double x) * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) */ if (__glibc_unlikely (ix > 0x4080)) - z = (invsqrtpi * cc) / __ieee754_sqrtl (y); + z = (invsqrtpi * cc) / sqrtl (y); else { u = pone (y); v = qone (y); - z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y); + z = invsqrtpi * (u * cc - v * ss) / sqrtl (y); } if (se & 0x8000) return -z; @@ -195,7 +196,7 @@ __ieee754_y1l (long double x) { long double z, s, c, ss, cc, u, v; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -231,12 +232,12 @@ __ieee754_y1l (long double x) * to compute the worse one. */ if (__glibc_unlikely (ix > 0x4080)) - z = (invsqrtpi * ss) / __ieee754_sqrtl (x); + z = (invsqrtpi * ss) / sqrtl (x); else { u = pone (x); v = qone (x); - z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x); + z = invsqrtpi * (u * ss + v * cc) / sqrtl (x); } return z; } @@ -362,7 +363,7 @@ pone (long double x) const long double *p, *q; long double z, r, s; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -507,9 +508,9 @@ static long double qone (long double x) { const long double *p, *q; - static long double s, r, z; + long double s, r, z; int32_t ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; diff --git a/sysdeps/ieee754/ldbl-96/e_jnl.c b/sysdeps/ieee754/ldbl-96/e_jnl.c index 92f96921a7..394921f564 100644 --- a/sysdeps/ieee754/ldbl-96/e_jnl.c +++ b/sysdeps/ieee754/ldbl-96/e_jnl.c @@ -60,6 +60,7 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const long double invsqrtpi = 5.64189583547756286948079e-1L, two = 2.0e0L, one = 1.0e0L; @@ -69,7 +70,7 @@ static const long double zero = 0.0L; long double __ieee754_jnl (int n, long double x) { - u_int32_t se, i0, i1; + uint32_t se, i0, i1; int32_t i, ix, sgn; long double a, b, temp, di, ret; long double z, w; @@ -142,7 +143,7 @@ __ieee754_jnl (int n, long double x) temp = c - s; break; } - b = invsqrtpi * temp / __ieee754_sqrtl (x); + b = invsqrtpi * temp / sqrtl (x); } else { @@ -302,7 +303,7 @@ strong_alias (__ieee754_jnl, __jnl_finite) long double __ieee754_ynl (int n, long double x) { - u_int32_t se, i0, i1; + uint32_t se, i0, i1; int32_t i, ix; int32_t sign; long double a, b, temp, ret; @@ -371,7 +372,7 @@ __ieee754_ynl (int n, long double x) temp = s + c; break; } - b = invsqrtpi * temp / __ieee754_sqrtl (x); + b = invsqrtpi * temp / sqrtl (x); } else { diff --git a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c index 9862fe8d5c..200421f5cc 100644 --- a/sysdeps/ieee754/ldbl-96/e_lgammal_r.c +++ b/sysdeps/ieee754/ldbl-96/e_lgammal_r.c @@ -91,9 +91,9 @@ * */ -#include <libc-internal.h> #include <math.h> #include <math_private.h> +#include <libc-diag.h> static const long double half = 0.5L, @@ -208,7 +208,7 @@ sin_pi (long double x) { long double y, z; int n, ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -275,7 +275,7 @@ __ieee754_lgammal_r (long double x, int *signgamp) { long double t, y, z, nadj, p, p1, p2, q, r, w; int i, ix; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; *signgamp = 1; GET_LDOUBLE_WORDS (se, i0, i1, x); diff --git a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c index 0d8e64675c..f67805f2d3 100644 --- a/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c +++ b/sysdeps/ieee754/ldbl-96/e_rem_pio2l.c @@ -1,5 +1,5 @@ /* Extended-precision floating point argument reduction. - Copyright (C) 1999-2016 Free Software Foundation, Inc. + Copyright (C) 1999-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision code by Jakub Jelinek <jj@ultra.linux.cz> @@ -186,7 +186,7 @@ __ieee754_rem_pio2l (long double x, long double *y) { double tx[3], ty[3]; int32_t se, j0; - u_int32_t i0, i1; + uint32_t i0, i1; int sx; int n, exp; diff --git a/sysdeps/ieee754/ldbl-96/e_sinhl.c b/sysdeps/ieee754/ldbl-96/e_sinhl.c index 095b142621..a4b39783bc 100644 --- a/sysdeps/ieee754/ldbl-96/e_sinhl.c +++ b/sysdeps/ieee754/ldbl-96/e_sinhl.c @@ -39,6 +39,7 @@ static char rcsid[] = "$NetBSD: $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> static const long double one = 1.0, shuge = 1.0e4931L; @@ -46,7 +47,7 @@ long double __ieee754_sinhl(long double x) { long double t,w,h; - u_int32_t jx,ix,i0,i1; + uint32_t jx,ix,i0,i1; /* Words of |x|. */ GET_LDOUBLE_WORDS(jx,i0,i1,x); diff --git a/sysdeps/ieee754/ldbl-96/gamma_product.c b/sysdeps/ieee754/ldbl-96/gamma_product.c index 419d11598f..f1b65e12e2 100644 --- a/sysdeps/ieee754/ldbl-96/gamma_product.c +++ b/sysdeps/ieee754/ldbl-96/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 @@ -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/ldbl-96/gamma_productl.c b/sysdeps/ieee754/ldbl-96/gamma_productl.c index 849b57d95d..ed0c166d78 100644 --- a/sysdeps/ieee754/ldbl-96/gamma_productl.c +++ b/sysdeps/ieee754/ldbl-96/gamma_productl.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 inline void -mul_split (long double *hi, long double *lo, long double x, long double y) -{ -#ifdef __FP_FAST_FMAL - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fmal (x, y, -*hi); -#elif defined FP_FAST_FMAL - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fmal (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) - long double x1 = x * C; - long double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - long double x2 = x - x1; - long double y2 = y - y1; - *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; -#endif -} +#include <mul_splitl.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 @@ -68,7 +38,7 @@ __gamma_productl (long double x, long double x_eps, int n, long double *eps) { *eps += x_eps / (x + i); long double lo; - mul_split (&ret, &lo, ret, x + i); + mul_splitl (&ret, &lo, ret, x + i); *eps += lo / ret; } return ret; diff --git a/sysdeps/ieee754/ldbl-96/include/bits/iscanonical.h b/sysdeps/ieee754/ldbl-96/include/bits/iscanonical.h new file mode 100644 index 0000000000..bee080bd29 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/include/bits/iscanonical.h @@ -0,0 +1,5 @@ +#include_next <bits/iscanonical.h> + +#ifndef _ISOMAC +libm_hidden_proto (__iscanonicall) +#endif diff --git a/sysdeps/ieee754/ldbl-96/k_cosl.c b/sysdeps/ieee754/ldbl-96/k_cosl.c index 08b11b3733..da20385210 100644 --- a/sysdeps/ieee754/ldbl-96/k_cosl.c +++ b/sysdeps/ieee754/ldbl-96/k_cosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point cosine on <-pi/4,pi/4>. - Copyright (C) 1999-2016 Free Software Foundation, Inc. + Copyright (C) 1999-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision cosine by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/k_sinl.c b/sysdeps/ieee754/ldbl-96/k_sinl.c index 6ba7ceddc0..2549f71d19 100644 --- a/sysdeps/ieee754/ldbl-96/k_sinl.c +++ b/sysdeps/ieee754/ldbl-96/k_sinl.c @@ -1,5 +1,5 @@ /* Quad-precision floating point sine on <-pi/4,pi/4>. - Copyright (C) 1999-2016 Free Software Foundation, Inc. + Copyright (C) 1999-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision sine by Jakub Jelinek <jj@ultra.linux.cz> @@ -23,6 +23,7 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> /* The polynomials have not been optimized for extended-precision and may contain more terms than needed. */ diff --git a/sysdeps/ieee754/ldbl-96/k_tanl.c b/sysdeps/ieee754/ldbl-96/k_tanl.c index 0c050c112f..9b5151baa2 100644 --- a/sysdeps/ieee754/ldbl-96/k_tanl.c +++ b/sysdeps/ieee754/ldbl-96/k_tanl.c @@ -57,9 +57,11 @@ */ #include <float.h> -#include <libc-internal.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libc-diag.h> + static const long double one = 1.0L, pio4hi = 0xc.90fdaa22168c235p-4L, diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c index fe7002f640..0805051723 100644 --- a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c +++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.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/ldbl-96/lgamma_negl.c b/sysdeps/ieee754/ldbl-96/lgamma_negl.c index 99ecf7e85f..6d2e0b7165 100644 --- a/sysdeps/ieee754/ldbl-96/lgamma_negl.c +++ b/sysdeps/ieee754/ldbl-96/lgamma_negl.c @@ -1,5 +1,5 @@ /* lgammal 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 diff --git a/sysdeps/ieee754/ldbl-96/lgamma_product.c b/sysdeps/ieee754/ldbl-96/lgamma_product.c index e3ba72d8e4..1eb99c40a5 100644 --- a/sysdeps/ieee754/ldbl-96/lgamma_product.c +++ b/sysdeps/ieee754/ldbl-96/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 diff --git a/sysdeps/ieee754/ldbl-96/lgamma_productl.c b/sysdeps/ieee754/ldbl-96/lgamma_productl.c index de67cbe665..9141a3177a 100644 --- a/sysdeps/ieee754/ldbl-96/lgamma_productl.c +++ b/sysdeps/ieee754/ldbl-96/lgamma_productl.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 (long double *hi, long double *lo, long double x, long double y) -{ -#ifdef __FP_FAST_FMAL - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fmal (x, y, -*hi); -#elif defined FP_FAST_FMAL - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fmal (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) - long double x1 = x * C; - long double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - long double x2 = x - x1; - long double y2 = y - y1; - *lo = (((x1 * y1 - *hi) + x1 * y2) + x2 * y1) + x2 * y2; -#endif -} +#include <mul_splitl.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 @@ -65,11 +35,11 @@ __lgamma_productl (long double t, long double x, long double x_eps, int n) long double xi = x + i; long double quot = t / xi; long double mhi, mlo; - mul_split (&mhi, &mlo, quot, xi); + mul_splitl (&mhi, &mlo, quot, xi); long double quot_lo = (t - mhi - mlo) / xi - t * x_eps / (xi * xi); /* We want (1 + RET + RET_EPS) * (1 + QUOT + QUOT_LO) - 1. */ long double rhi, rlo; - mul_split (&rhi, &rlo, ret, quot); + mul_splitl (&rhi, &rlo, ret, quot); long double rpq = ret + quot; long double rpq_eps = (ret - rpq) + quot; long double nret = rpq + rhi; diff --git a/sysdeps/ieee754/ldbl-96/strtod_nan_ldouble.h b/sysdeps/ieee754/ldbl-96/math-nan-payload-ldouble.h index 2694b5ee34..ab2542c097 100644 --- a/sysdeps/ieee754/ldbl-96/strtod_nan_ldouble.h +++ b/sysdeps/ieee754/ldbl-96/math-nan-payload-ldouble.h @@ -1,5 +1,5 @@ -/* Convert string for NaN payload to corresponding NaN. For ldbl-96. - Copyright (C) 1997-2016 Free Software Foundation, Inc. +/* NaN payload handling for ldbl-96. + 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 @@ -16,8 +16,7 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ -#define FLOAT long double -#define SET_MANTISSA(flt, mant) \ +#define SET_NAN_PAYLOAD(flt, mant) \ do \ { \ union ieee854_long_double u; \ diff --git a/sysdeps/ieee754/ldbl-96/math_ldbl.h b/sysdeps/ieee754/ldbl-96/math_ldbl.h index cca30657ce..99428f6eeb 100644 --- a/sysdeps/ieee754/ldbl-96/math_ldbl.h +++ b/sysdeps/ieee754/ldbl-96/math_ldbl.h @@ -1,11 +1,31 @@ -#ifndef _MATH_PRIVATE_H_ -#error "Never use <math_ldbl.h> directly; include <math_private.h> instead." -#endif +/* Manipulation of the bit representation of 'long double' quantities. + Copyright (C) 1999-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/>. */ + +#ifndef _MATH_LDBL_H_ +#define _MATH_LDBL_H_ 1 + +#include <stdint.h> +#include <endian.h> /* A union which permits us to convert between a long double and three 32 bit ints. */ -#if __FLOAT_WORD_ORDER == BIG_ENDIAN +#if __FLOAT_WORD_ORDER == __BIG_ENDIAN typedef union { @@ -14,22 +34,22 @@ typedef union { int sign_exponent:16; unsigned int empty:16; - u_int32_t msw; - u_int32_t lsw; + uint32_t msw; + uint32_t lsw; } parts; } ieee_long_double_shape_type; #endif -#if __FLOAT_WORD_ORDER == LITTLE_ENDIAN +#if __FLOAT_WORD_ORDER == __LITTLE_ENDIAN typedef union { long double value; struct { - u_int32_t lsw; - u_int32_t msw; + uint32_t lsw; + uint32_t msw; int sign_exponent:16; unsigned int empty:16; } parts; @@ -96,3 +116,5 @@ do { \ se_u.parts.sign_exponent = (exp); \ (d) = se_u.value; \ } while (0) + +#endif /* math_ldbl.h */ diff --git a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c index 6159d7dabf..8cbf28256b 100644 --- a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c +++ b/sysdeps/ieee754/ldbl-96/mpn2ldbl.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/ldbl-96/printf_fphex.c b/sysdeps/ieee754/ldbl-96/printf_fphex.c index 65bf538590..45d97447a7 100644 --- a/sysdeps/ieee754/ldbl-96/printf_fphex.c +++ b/sysdeps/ieee754/ldbl-96/printf_fphex.c @@ -1,5 +1,5 @@ /* Print floating point number in hexadecimal notation according to ISO C99. - 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 diff --git a/sysdeps/ieee754/ldbl-96/s_asinhl.c b/sysdeps/ieee754/ldbl-96/s_asinhl.c index da49ea5988..2b9ae1f677 100644 --- a/sysdeps/ieee754/ldbl-96/s_asinhl.c +++ b/sysdeps/ieee754/ldbl-96/s_asinhl.c @@ -32,6 +32,8 @@ static char rcsid[] = "$NetBSD: $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-ldouble.h> static const long double one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */ @@ -54,12 +56,12 @@ long double __asinhl(long double x) } else { long double xa = fabsl(x); if (ix>0x4000) { /* 2**34 > |x| > 2.0 */ - w = __ieee754_logl(2.0*xa+one/(__ieee754_sqrtl(xa*xa+one)+xa)); + w = __ieee754_logl(2.0*xa+one/(sqrtl(xa*xa+one)+xa)); } else { /* 2.0 > |x| > 2**-28 */ t = xa*xa; - w =__log1pl(xa+t/(one+__ieee754_sqrtl(one+t))); + w =__log1pl(xa+t/(one+sqrtl(one+t))); } } return __copysignl(w, x); } -weak_alias (__asinhl, asinhl) +libm_alias_ldouble (__asinh, asinh) diff --git a/sysdeps/ieee754/ldbl-96/s_cbrtl.c b/sysdeps/ieee754/ldbl-96/s_cbrtl.c index 42849ab517..67cf86dd7a 100644 --- a/sysdeps/ieee754/ldbl-96/s_cbrtl.c +++ b/sysdeps/ieee754/ldbl-96/s_cbrtl.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-ldouble.h> #define CBRT2 1.2599210498948731648 /* 2^(1/3) */ @@ -67,4 +68,4 @@ __cbrtl (long double x) u -= (u - (x / (u * u))) * third; return u; } -weak_alias (__cbrtl, cbrtl) +libm_alias_ldouble (__cbrt, cbrt) diff --git a/sysdeps/ieee754/ldbl-96/s_copysignl.c b/sysdeps/ieee754/ldbl-96/s_copysignl.c index b1c442452f..3c16d54783 100644 --- a/sysdeps/ieee754/ldbl-96/s_copysignl.c +++ b/sysdeps/ieee754/ldbl-96/s_copysignl.c @@ -26,13 +26,14 @@ static char rcsid[] = "$NetBSD: $"; #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long double __copysignl(long double x, long double y) { - u_int32_t es1,es2; + uint32_t es1,es2; GET_LDOUBLE_EXP(es1,x); GET_LDOUBLE_EXP(es2,y); SET_LDOUBLE_EXP(x,(es1&0x7fff)|(es2&0x8000)); return x; } -weak_alias (__copysignl, copysignl) +libm_alias_ldouble (__copysign, copysign) diff --git a/sysdeps/ieee754/ldbl-96/s_cosl.c b/sysdeps/ieee754/ldbl-96/s_cosl.c index 8b0b7d3cc2..324e5b9663 100644 --- a/sysdeps/ieee754/ldbl-96/s_cosl.c +++ b/sysdeps/ieee754/ldbl-96/s_cosl.c @@ -52,6 +52,7 @@ static char rcsid[] = "$NetBSD: $"; #include <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long double __cosl(long double x) { @@ -85,4 +86,4 @@ long double __cosl(long double x) } } } -weak_alias (__cosl, cosl) +libm_alias_ldouble (__cos, cos) diff --git a/sysdeps/ieee754/ldbl-96/s_daddl.c b/sysdeps/ieee754/ldbl-96/s_daddl.c new file mode 100644 index 0000000000..d1e3c17bc7 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_daddl.c @@ -0,0 +1,33 @@ +/* Add long double (ldbl-96) values, narrowing the result to double. + 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 f32xaddf64x __hide_f32xaddf64x +#define f64addf64x __hide_f64addf64x +#include <math.h> +#undef f32xaddf64x +#undef f64addf64x + +#include <math-narrow.h> + +double +__daddl (long double x, long double y) +{ + NARROW_ADD_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l, + mantissa1); +} +libm_alias_double_ldouble (add) diff --git a/sysdeps/ieee754/ldbl-96/s_ddivl.c b/sysdeps/ieee754/ldbl-96/s_ddivl.c new file mode 100644 index 0000000000..9c266d1ff3 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_ddivl.c @@ -0,0 +1,33 @@ +/* Divide long double (ldbl-96) values, narrowing the result to double. + 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 f32xdivf64x __hide_f32xdivf64x +#define f64divf64x __hide_f64divf64x +#include <math.h> +#undef f32xdivf64x +#undef f64divf64x + +#include <math-narrow.h> + +double +__ddivl (long double x, long double y) +{ + NARROW_DIV_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l, + mantissa1); +} +libm_alias_double_ldouble (div) diff --git a/sysdeps/ieee754/ldbl-96/s_dmull.c b/sysdeps/ieee754/ldbl-96/s_dmull.c new file mode 100644 index 0000000000..a717b0aa07 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_dmull.c @@ -0,0 +1,33 @@ +/* Multiply long double (ldbl-96) values, narrowing the result to double. + 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 f32xmulf64x __hide_f32xmulf64x +#define f64mulf64x __hide_f64mulf64x +#include <math.h> +#undef f32xmulf64x +#undef f64mulf64x + +#include <math-narrow.h> + +double +__dmull (long double x, long double y) +{ + NARROW_MUL_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l, + mantissa1); +} +libm_alias_double_ldouble (mul) diff --git a/sysdeps/ieee754/ldbl-96/s_dsubl.c b/sysdeps/ieee754/ldbl-96/s_dsubl.c new file mode 100644 index 0000000000..5a855790f6 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_dsubl.c @@ -0,0 +1,33 @@ +/* Subtract long double (ldbl-96) values, narrowing the result to double. + 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 f32xsubf64x __hide_f32xsubf64x +#define f64subf64x __hide_f64subf64x +#include <math.h> +#undef f32xsubf64x +#undef f64subf64x + +#include <math-narrow.h> + +double +__dsubl (long double x, long double y) +{ + NARROW_SUB_ROUND_TO_ODD (x, y, double, union ieee854_long_double, l, + mantissa1); +} +libm_alias_double_ldouble (sub) diff --git a/sysdeps/ieee754/ldbl-96/s_erfl.c b/sysdeps/ieee754/ldbl-96/s_erfl.c index d00adb1000..1e42df70a7 100644 --- a/sysdeps/ieee754/ldbl-96/s_erfl.c +++ b/sysdeps/ieee754/ldbl-96/s_erfl.c @@ -108,6 +108,8 @@ #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-ldouble.h> static const long double tiny = 1e-4931L, @@ -254,7 +256,7 @@ __erfl (long double x) { long double R, S, P, Q, s, y, z, r; int32_t ix, i; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -335,13 +337,13 @@ __erfl (long double x) return r / x - one; } -weak_alias (__erfl, erfl) +libm_alias_ldouble (__erf, erf) long double __erfcl (long double x) { int32_t hx, ix; long double R, S, P, Q, s, y, z, r; - u_int32_t se, i0, i1; + uint32_t se, i0, i1; GET_LDOUBLE_WORDS (se, i0, i1, x); ix = se & 0x7fff; @@ -448,4 +450,4 @@ __erfcl (long double x) } } -weak_alias (__erfcl, erfcl) +libm_alias_ldouble (__erfc, erfc) diff --git a/sysdeps/ieee754/ldbl-96/s_faddl.c b/sysdeps/ieee754/ldbl-96/s_faddl.c new file mode 100644 index 0000000000..4164774cd4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_faddl.c @@ -0,0 +1,31 @@ +/* Add long double (ldbl-96) 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 f32addf64x __hide_f32addf64x +#include <math.h> +#undef f32addf64x + +#include <math-narrow.h> + +float +__faddl (long double x, long double y) +{ + NARROW_ADD_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l, + mantissa1); +} +libm_alias_float_ldouble (add) diff --git a/sysdeps/ieee754/ldbl-96/s_fdivl.c b/sysdeps/ieee754/ldbl-96/s_fdivl.c new file mode 100644 index 0000000000..ccb87ccd15 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fdivl.c @@ -0,0 +1,31 @@ +/* Divide long double (ldbl-96) 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 f32divf64x __hide_f32divf64x +#include <math.h> +#undef f32divf64x + +#include <math-narrow.h> + +float +__fdivl (long double x, long double y) +{ + NARROW_DIV_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l, + mantissa1); +} +libm_alias_float_ldouble (div) diff --git a/sysdeps/ieee754/ldbl-96/s_fma.c b/sysdeps/ieee754/ldbl-96/s_fma.c index ab45bcfce2..986879cda5 100644 --- a/sysdeps/ieee754/ldbl-96/s_fma.c +++ b/sysdeps/ieee754/ldbl-96/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> /* This implementation uses rounding to odd to avoid problems with double rounding. See a paper by Boldo and Melquiond: @@ -30,14 +32,12 @@ double __fma (double x, double y, double z) { - if (__glibc_unlikely (isinf (z))) - { - /* If z is Inf, but x and y are finite, the result should be - z rather than NaN. */ - if (isfinite (x) && isfinite (y)) - return (z + x) + y; - return (x * y) + z; - } + if (__glibc_unlikely (!isfinite (x) || !isfinite (y))) + return x * y + z; + else if (__glibc_unlikely (!isfinite (z))) + /* If z is Inf, but x and y are finite, the result should be z + rather than NaN. */ + return (z + x) + y; /* Ensure correct sign of exact 0 + 0. */ if (__glibc_unlikely ((x == 0 || y == 0) && z == 0)) @@ -97,5 +97,5 @@ __fma (double x, double y, double z) return u.d; } #ifndef __fma -weak_alias (__fma, fma) +libm_alias_double (__fma, fma) #endif diff --git a/sysdeps/ieee754/ldbl-96/s_fmal.c b/sysdeps/ieee754/ldbl-96/s_fmal.c index f1467fda3d..0b261fd17a 100644 --- a/sysdeps/ieee754/ldbl-96/s_fmal.c +++ b/sysdeps/ieee754/ldbl-96/s_fmal.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-ldouble.h> #include <tininess.h> /* This implementation uses rounding to odd to avoid problems with @@ -293,4 +295,4 @@ __fmal (long double x, long double y, long double z) return v.d * 0x1p-130L; } } -weak_alias (__fmal, fmal) +libm_alias_ldouble (__fma, fma) diff --git a/sysdeps/ieee754/ldbl-96/s_fmull.c b/sysdeps/ieee754/ldbl-96/s_fmull.c new file mode 100644 index 0000000000..b7582526a6 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fmull.c @@ -0,0 +1,31 @@ +/* Multiply long double (ldbl-96) 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 f32mulf64x __hide_f32mulf64x +#include <math.h> +#undef f32mulf64x + +#include <math-narrow.h> + +float +__fmull (long double x, long double y) +{ + NARROW_MUL_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l, + mantissa1); +} +libm_alias_float_ldouble (mul) diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c index ab217a659b..7c31ed9936 100644 --- a/sysdeps/ieee754/ldbl-96/s_frexpl.c +++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c @@ -31,6 +31,7 @@ static char rcsid[] = "$NetBSD: $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> static const long double #if LDBL_MANT_DIG == 64 @@ -42,11 +43,11 @@ two65 = 3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */ long double __frexpl(long double x, int *eptr) { - u_int32_t se, hx, ix, lx; + uint32_t se, hx, ix, lx; GET_LDOUBLE_WORDS(se,hx,lx,x); ix = 0x7fff&se; *eptr = 0; - if(ix==0x7fff||((ix|hx|lx)==0)) return x; /* 0,inf,nan */ + if(ix==0x7fff||((ix|hx|lx)==0)) return x + x; /* 0,inf,nan */ if (ix==0x0000) { /* subnormal */ x *= two65; GET_LDOUBLE_EXP(se,x); @@ -58,4 +59,4 @@ long double __frexpl(long double x, int *eptr) SET_LDOUBLE_EXP(x,se); return x; } -weak_alias (__frexpl, frexpl) +libm_alias_ldouble (__frexp, frexp) diff --git a/sysdeps/ieee754/ldbl-96/s_fromfpl.c b/sysdeps/ieee754/ldbl-96/s_fromfpl.c new file mode 100644 index 0000000000..bcedceea8e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fromfpl.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 0 +#define FUNC __fromfpl +#include <s_fromfpl_main.c> +libm_alias_ldouble (__fromfp, fromfp) diff --git a/sysdeps/ieee754/ldbl-96/s_fromfpl_main.c b/sysdeps/ieee754/ldbl-96/s_fromfpl_main.c new file mode 100644 index 0000000000..6f24ccf488 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fromfpl_main.c @@ -0,0 +1,85 @@ +/* Round to integer type. ldbl-96 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-ldouble.h> +#include <stdbool.h> +#include <stdint.h> + +#define BIAS 0x3fff +#define MANT_DIG 64 + +#if UNSIGNED +# define RET_TYPE uintmax_t +#else +# define RET_TYPE intmax_t +#endif + +#include <fromfp.h> + +RET_TYPE +FUNC (long double x, int round, unsigned int width) +{ + if (width > INTMAX_WIDTH) + width = INTMAX_WIDTH; + uint16_t se; + uint32_t hx, lx; + GET_LDOUBLE_WORDS (se, hx, lx, x); + bool negative = (se & 0x8000) != 0; + if (width == 0) + return fromfp_domain_error (negative, width); + if ((hx | lx) == 0) + return 0; + int exponent = se & 0x7fff; + exponent -= BIAS; + int max_exponent = fromfp_max_exponent (negative, width); + if (exponent > max_exponent) + return fromfp_domain_error (negative, width); + + uint64_t ix = (((uint64_t) hx) << 32) | lx; + uintmax_t uret; + bool half_bit, more_bits; + if (exponent >= MANT_DIG - 1) + { + uret = ix; + /* Exponent 63; no shifting required. */ + 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; + if (exponent == -1) + uret = 0; + else + 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/ldbl-96/s_fromfpxl.c b/sysdeps/ieee754/ldbl-96/s_fromfpxl.c new file mode 100644 index 0000000000..0a342a22d1 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fromfpxl.c @@ -0,0 +1,5 @@ +#define UNSIGNED 0 +#define INEXACT 1 +#define FUNC __fromfpxl +#include <s_fromfpl_main.c> +libm_alias_ldouble (__fromfpx, fromfpx) diff --git a/sysdeps/ieee754/ldbl-96/s_fsubl.c b/sysdeps/ieee754/ldbl-96/s_fsubl.c new file mode 100644 index 0000000000..54aaf68b74 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fsubl.c @@ -0,0 +1,31 @@ +/* Subtract long double (ldbl-96) 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 f32subf64x __hide_f32subf64x +#include <math.h> +#undef f32subf64x + +#include <math-narrow.h> + +float +__fsubl (long double x, long double y) +{ + NARROW_SUB_ROUND_TO_ODD (x, y, float, union ieee854_long_double, l, + mantissa1); +} +libm_alias_float_ldouble (sub) diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/s_getpayloadl.c index 5c7a1812bc..4b7b734f3d 100644 --- a/sysdeps/ieee754/ldbl-96/w_expl.c +++ b/sysdeps/ieee754/ldbl-96/s_getpayloadl.c @@ -1,6 +1,6 @@ -/* Copyright (C) 2011-2016 Free Software Foundation, Inc. +/* Get NaN payload. ldbl-96 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 @@ -18,17 +18,17 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> +#include <stdint.h> -/* wrapper expl */ long double -__expl (long double x) +__getpayloadl (const long double *x) { - long double z = __ieee754_expl (x); - if (__builtin_expect (!isfinite (z) || z == 0, 0) - && isfinite (x) && _LIB_VERSION != _IEEE_) - return __kernel_standard_l (x, x, 206 + !!signbit (x)); - - return z; + uint16_t se __attribute__ ((unused)); + uint32_t hx, lx; + GET_LDOUBLE_WORDS (se, hx, lx, *x); + hx &= 0x3fffffff; + uint64_t ix = ((uint64_t) hx << 32) | lx; + return (long double) ix; } -hidden_def (__expl) -weak_alias (__expl, expl) +libm_alias_ldouble (__getpayload, getpayload) diff --git a/sysdeps/ieee754/ldbl-96/s_iscanonicall.c b/sysdeps/ieee754/ldbl-96/s_iscanonicall.c new file mode 100644 index 0000000000..413c6bd42c --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_iscanonicall.c @@ -0,0 +1,44 @@ +/* Test whether long double value is canonical. ldbl-96 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 <float.h> +#include <math.h> +#include <math_private.h> +#include <stdbool.h> +#include <stdint.h> + +int +__iscanonicall (long double x) +{ + uint32_t se, i0, i1 __attribute__ ((unused)); + + GET_LDOUBLE_WORDS (se, i0, i1, x); + int32_t ix = se & 0x7fff; + bool mant_high = (i0 & 0x80000000) != 0; + + if (LDBL_MIN_EXP == -16381) + /* Intel variant: the high mantissa bit should have a value + determined by the exponent. */ + return ix > 0 ? mant_high : !mant_high; + else + /* M68K variant: both values of the high bit are valid for the + greatest and smallest exponents, while other exponents require + the high bit to be set. */ + return ix == 0 || ix == 0x7fff || mant_high; +} +libm_hidden_def (__iscanonicall) diff --git a/sysdeps/ieee754/ldbl-96/s_issignalingl.c b/sysdeps/ieee754/ldbl-96/s_issignalingl.c index 73646cac0c..b1dd41ecfd 100644 --- a/sysdeps/ieee754/ldbl-96/s_issignalingl.c +++ b/sysdeps/ieee754/ldbl-96/s_issignalingl.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 __issignalingl (long double x) { - u_int32_t exi, hxi, lxi; + uint32_t exi, hxi, lxi; GET_LDOUBLE_WORDS (exi, hxi, lxi, x); -#ifdef HIGH_ORDER_BIT_IS_SET_FOR_SNAN +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN # error not implemented #else /* To keep the following comparison simple, toggle the quiet/signaling bit, diff --git a/sysdeps/ieee754/ldbl-96/s_llrintl.c b/sysdeps/ieee754/ldbl-96/s_llrintl.c index 592d51c607..d45a69a1f7 100644 --- a/sysdeps/ieee754/ldbl-96/s_llrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_llrintl.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. @@ -23,6 +23,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> static const long double two63[2] = { @@ -35,7 +36,7 @@ long long int __llrintl (long double x) { int32_t se,j0; - u_int32_t i0, i1; + uint32_t i0, i1; long long int result; long double w; long double t; @@ -88,4 +89,4 @@ __llrintl (long double x) return sx ? -result : result; } -weak_alias (__llrintl, llrintl) +libm_alias_ldouble (__llrint, llrint) diff --git a/sysdeps/ieee754/ldbl-96/s_llroundl.c b/sysdeps/ieee754/ldbl-96/s_llroundl.c index 483199d442..601fd0e644 100644 --- a/sysdeps/ieee754/ldbl-96/s_llroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_llroundl.c @@ -1,5 +1,5 @@ /* Round long 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,13 +22,14 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long long int __llroundl (long double x) { int32_t j0; - u_int32_t se, i1, i0; + uint32_t se, i1, i0; long long int result; int sign; @@ -42,7 +43,7 @@ __llroundl (long double x) return j0 < -1 ? 0 : sign; else { - u_int32_t j = i0 + (0x40000000 >> j0); + uint32_t j = i0 + (0x40000000 >> j0); if (j < i0) { j >>= 1; @@ -59,7 +60,7 @@ __llroundl (long double x) result = (((long long int) i0 << 32) | i1) << (j0 - 63); else { - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + uint32_t j = i1 + (0x80000000 >> (j0 - 31)); result = (long long int) i0; if (j < i1) @@ -86,4 +87,4 @@ __llroundl (long double x) return sign * result; } -weak_alias (__llroundl, llroundl) +libm_alias_ldouble (__llround, llround) diff --git a/sysdeps/ieee754/ldbl-96/s_lrintl.c b/sysdeps/ieee754/ldbl-96/s_lrintl.c index bd902deb47..df3222c7f2 100644 --- a/sysdeps/ieee754/ldbl-96/s_lrintl.c +++ b/sysdeps/ieee754/ldbl-96/s_lrintl.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. @@ -23,6 +23,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> static const long double two63[2] = { @@ -35,7 +36,7 @@ long int __lrintl (long double x) { int32_t se,j0; - u_int32_t i0, i1; + uint32_t i0, i1; long int result; long double w; long double t; @@ -123,4 +124,4 @@ __lrintl (long double x) return sx ? -result : result; } -weak_alias (__lrintl, lrintl) +libm_alias_ldouble (__lrint, lrint) diff --git a/sysdeps/ieee754/ldbl-96/s_lroundl.c b/sysdeps/ieee754/ldbl-96/s_lroundl.c index 3b43d77e72..0cc9f9c5d6 100644 --- a/sysdeps/ieee754/ldbl-96/s_lroundl.c +++ b/sysdeps/ieee754/ldbl-96/s_lroundl.c @@ -1,5 +1,5 @@ /* Round long 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,13 +22,14 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long int __lroundl (long double x) { int32_t j0; - u_int32_t se, i1, i0; + uint32_t se, i1, i0; long int result; int sign; @@ -42,7 +43,7 @@ __lroundl (long double x) return j0 < -1 ? 0 : sign; else { - u_int32_t j = i0 + (0x40000000 >> j0); + uint32_t j = i0 + (0x40000000 >> j0); if (j < i0) { j >>= 1; @@ -66,7 +67,7 @@ __lroundl (long double x) result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); else { - u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + uint32_t j = i1 + (0x80000000 >> (j0 - 31)); unsigned long int ures = i0; if (j < i1) @@ -108,4 +109,4 @@ __lroundl (long double x) return sign * result; } -weak_alias (__lroundl, lroundl) +libm_alias_ldouble (__lround, lround) diff --git a/sysdeps/ieee754/ldbl-96/s_modfl.c b/sysdeps/ieee754/ldbl-96/s_modfl.c index e9401d0f5d..380b6f0389 100644 --- a/sysdeps/ieee754/ldbl-96/s_modfl.c +++ b/sysdeps/ieee754/ldbl-96/s_modfl.c @@ -26,6 +26,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> static const long double one = 1.0; @@ -33,7 +34,7 @@ long double __modfl(long double x, long double *iptr) { int32_t i0,i1,j0; - u_int32_t i,se; + uint32_t i,se; GET_LDOUBLE_WORDS(se,i0,i1,x); j0 = (se&0x7fff)-0x3fff; /* exponent of x */ if(j0<32) { /* integer part in high x */ @@ -59,7 +60,7 @@ __modfl(long double x, long double *iptr) SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */ return x; } else { /* fraction part in low x */ - i = ((u_int32_t)(0x7fffffff))>>(j0-32); + i = ((uint32_t)(0x7fffffff))>>(j0-32); if((i1&i)==0) { /* x is integral */ *iptr = x; SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */ @@ -70,4 +71,4 @@ __modfl(long double x, long double *iptr) } } } -weak_alias (__modfl, modfl) +libm_alias_ldouble (__modf, modf) diff --git a/sysdeps/ieee754/ldbl-96/s_nexttoward.c b/sysdeps/ieee754/ldbl-96/s_nexttoward.c index 3d0382eac9..1d8d9c7f91 100644 --- a/sysdeps/ieee754/ldbl-96/s_nexttoward.c +++ b/sysdeps/ieee754/ldbl-96/s_nexttoward.c @@ -27,13 +27,14 @@ static char rcsid[] = "$NetBSD: $"; #include <errno.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> #include <float.h> double __nexttoward(double x, long double y) { int32_t hx,ix,iy; - u_int32_t lx,hy,ly,esy; + uint32_t lx,hy,ly,esy; EXTRACT_WORDS(hx,lx,x); GET_LDOUBLE_WORDS(esy,hy,ly,y); diff --git a/sysdeps/ieee754/ldbl-96/s_nexttowardf.c b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c index ae7538942f..9a08e1c8ff 100644 --- a/sysdeps/ieee754/ldbl-96/s_nexttowardf.c +++ b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c @@ -19,13 +19,14 @@ static char rcsid[] = "$NetBSD: $"; #include <errno.h> #include <math.h> +#include <math-barriers.h> #include <math_private.h> #include <float.h> float __nexttowardf(float x, long double y) { int32_t hx,ix,iy; - u_int32_t hy,ly,esy; + uint32_t hy,ly,esy; GET_FLOAT_WORD(hx,x); GET_LDOUBLE_WORDS(esy,hy,ly,y); diff --git a/sysdeps/ieee754/ldbl-96/s_nextupl.c b/sysdeps/ieee754/ldbl-96/s_nextupl.c new file mode 100644 index 0000000000..5dff32ce73 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_nextupl.c @@ -0,0 +1,86 @@ +/* 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-ldouble.h> + +/* Return the least floating-point number greater than X. */ +long double +__nextupl (long double x) +{ + uint32_t hx, ix; + uint32_t lx; + int32_t esx; + + GET_LDOUBLE_WORDS (esx, hx, lx, x); + ix = esx & 0x7fff; + + if (((ix == 0x7fff) && (((hx & 0x7fffffff) | lx) != 0))) /* x is nan. */ + return x + x; + if ((ix | hx | lx) == 0) + return LDBL_TRUE_MIN; + if (esx >= 0) + { /* x > 0. */ + if (isinf (x)) + return x; + lx += 1; + if (lx == 0) + { + hx += 1; +#if LDBL_MIN_EXP == -16381 + if (hx == 0 || (esx == 0 && hx == 0x80000000)) +#else + if (hx == 0) +#endif + { + esx += 1; + hx |= 0x80000000; + } + } + } + else + { /* x < 0. */ + if (lx == 0) + { +#if LDBL_MIN_EXP == -16381 + if (hx <= 0x80000000 && esx != 0xffff8000) + { + esx -= 1; + hx = hx - 1; + if ((esx & 0x7fff) > 0) + hx |= 0x80000000; + } + else + hx -= 1; +#else + if (ix != 0 && hx == 0x80000000) + hx = 0; + if (hx == 0) + esx -= 1; + hx -= 1; +#endif + } + lx -= 1; + } + SET_LDOUBLE_WORDS (x, esx, hx, lx); + return x; +} + +libm_alias_ldouble (__nextup, nextup) diff --git a/sysdeps/ieee754/ldbl-96/s_remquol.c b/sysdeps/ieee754/ldbl-96/s_remquol.c index 89b2630d46..88c5ea2084 100644 --- a/sysdeps/ieee754/ldbl-96/s_remquol.c +++ b/sysdeps/ieee754/ldbl-96/s_remquol.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-ldouble.h> static const long double zero = 0.0; @@ -29,7 +30,7 @@ long double __remquol (long double x, long double p, int *quo) { int32_t ex,ep,hx,hp; - u_int32_t sx,lx,lp; + uint32_t sx,lx,lp; int cquo,qs; GET_LDOUBLE_WORDS (ex, hx, lx, x); @@ -108,4 +109,4 @@ __remquol (long double x, long double p, int *quo) x = -x; return x; } -weak_alias (__remquol, remquol) +libm_alias_ldouble (__remquo, remquo) diff --git a/sysdeps/ieee754/ldbl-96/s_roundevenl.c b/sysdeps/ieee754/ldbl-96/s_roundevenl.c new file mode 100644 index 0000000000..be2e4fa49e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_roundevenl.c @@ -0,0 +1,126 @@ +/* Round to nearest integer value, rounding halfway cases to even. + ldbl-96 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-ldouble.h> +#include <stdint.h> + +#define BIAS 0x3fff +#define MANT_DIG 64 +#define MAX_EXP (2 * BIAS + 1) + +long double +__roundevenl (long double x) +{ + uint16_t se; + uint32_t hx, lx; + GET_LDOUBLE_WORDS (se, hx, lx, x); + int exponent = se & 0x7fff; + 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) + { + /* No need to test whether HALF_BIT is set. */ + lx += half_bit; + if (lx < half_bit) + { + hx++; + if (hx == 0) + { + hx = 0x80000000; + se++; + } + } + } + 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; + if (lx < 0x80000000) + { + hx++; + if (hx == 0) + { + hx = 0x80000000; + se++; + } + } + } + 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. */ + 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; + if (hx < half_bit) + { + hx = 0x80000000; + se++; + } + } + hx &= ~(int_bit - 1); + lx = 0; + } + else if (exponent == BIAS - 1 && (hx > 0x80000000 || lx != 0)) + { + /* Interval (0.5, 1). */ + se = (se & 0x8000) | 0x3fff; + hx = 0x80000000; + lx = 0; + } + else + { + /* Rounds to 0. */ + se &= 0x8000; + hx = 0; + lx = 0; + } + SET_LDOUBLE_WORDS (x, se, hx, lx); + return x; +} +libm_alias_ldouble (__roundeven, roundeven) diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c index 4f35c4847b..c5c304cb2e 100644 --- a/sysdeps/ieee754/ldbl-96/s_roundl.c +++ b/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -1,5 +1,5 @@ /* Round long 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 long double huge = 1.0e4930L; +#include <libm-alias-ldouble.h> long double __roundl (long double x) { int32_t j0; - u_int32_t se, i1, i0; + uint32_t se, i1, i0; GET_LDOUBLE_WORDS (se, i0, i1, x); j0 = (se & 0x7fff) - 0x3fff; @@ -37,7 +35,6 @@ __roundl (long double x) { if (j0 < 0) { - math_force_eval (huge + x); se &= 0x8000; i0 = i1 = 0; if (j0 == -1) @@ -48,14 +45,12 @@ __roundl (long double x) } else { - u_int32_t i = 0x7fffffff >> j0; + uint32_t i = 0x7fffffff >> j0; if (((i0 & i) | i1) == 0) /* X is integral. */ return x; - /* Raise inexact if x != 0. */ - math_force_eval (huge + x); - u_int32_t j = i0 + (0x40000000 >> j0); + uint32_t j = i0 + (0x40000000 >> j0); if (j < i0) se += 1; i0 = (j & ~i) | 0x80000000; @@ -72,17 +67,15 @@ __roundl (long double x) } else { - u_int32_t i = 0xffffffff >> (j0 - 31); + uint32_t i = 0xffffffff >> (j0 - 31); 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 << (62 - j0)); + uint32_t j = i1 + (1 << (62 - j0)); if (j < i1) { - u_int32_t k = i0 + 1; + uint32_t k = i0 + 1; if (k < i0) { se += 1; @@ -97,4 +90,4 @@ __roundl (long double x) SET_LDOUBLE_WORDS (x, se, i0, i1); return x; } -weak_alias (__roundl, roundl) +libm_alias_ldouble (__round, round) diff --git a/sysdeps/ieee754/ldbl-96/s_setpayloadl.c b/sysdeps/ieee754/ldbl-96/s_setpayloadl.c new file mode 100644 index 0000000000..9f43c259ec --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_setpayloadl.c @@ -0,0 +1,4 @@ +#define SIG 0 +#define FUNC __setpayloadl +#include <s_setpayloadl_main.c> +libm_alias_ldouble (__setpayload, setpayload) diff --git a/sysdeps/ieee754/ldbl-96/s_setpayloadl_main.c b/sysdeps/ieee754/ldbl-96/s_setpayloadl_main.c new file mode 100644 index 0000000000..6b85bdf1a9 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_setpayloadl_main.c @@ -0,0 +1,69 @@ +/* Set NaN payload. ldbl-96 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-ldouble.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 0x3fff +#define PAYLOAD_DIG 62 +#define EXPLICIT_MANT_DIG 63 + +int +FUNC (long double *x, long double payload) +{ + uint32_t hx, lx; + uint16_t exponent; + GET_LDOUBLE_WORDS (exponent, hx, lx, payload); + /* 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 + && exponent == 0 && hx == 0 && lx == 0))) + { + SET_LDOUBLE_WORDS (*x, 0, 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)) + { + SET_LDOUBLE_WORDS (*x, 0, 0, 0); + return 1; + } + if (exponent != 0) + { + if (shift >= 32) + { + lx = hx >> (shift - 32); + hx = 0; + } + else if (shift != 0) + { + lx = (lx >> shift) | (hx << (32 - shift)); + hx >>= shift; + } + } + hx |= 0x80000000 | (SET_HIGH_BIT ? 0x40000000 : 0); + SET_LDOUBLE_WORDS (*x, 0x7fff, hx, lx); + return 0; +} diff --git a/sysdeps/ieee754/ldbl-96/s_setpayloadsigl.c b/sysdeps/ieee754/ldbl-96/s_setpayloadsigl.c new file mode 100644 index 0000000000..cd82f295aa --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_setpayloadsigl.c @@ -0,0 +1,4 @@ +#define SIG 1 +#define FUNC __setpayloadsigl +#include <s_setpayloadl_main.c> +libm_alias_ldouble (__setpayloadsig, setpayloadsig) diff --git a/sysdeps/ieee754/ldbl-96/s_signbitl.c b/sysdeps/ieee754/ldbl-96/s_signbitl.c index ee5d77e27a..19953c180a 100644 --- a/sysdeps/ieee754/ldbl-96/s_signbitl.c +++ b/sysdeps/ieee754/ldbl-96/s_signbitl.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/ldbl-96/s_sincosl.c b/sysdeps/ieee754/ldbl-96/s_sincosl.c index ab32b73e7d..355c25dba9 100644 --- a/sysdeps/ieee754/ldbl-96/s_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/s_sincosl.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,6 +21,7 @@ #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> void @@ -73,4 +74,4 @@ __sincosl (long double x, long double *sinx, long double *cosx) } } } -weak_alias (__sincosl, sincosl) +libm_alias_ldouble (__sincos, sincos) diff --git a/sysdeps/ieee754/ldbl-96/s_sinl.c b/sysdeps/ieee754/ldbl-96/s_sinl.c index 11e1899822..cfbe9bf153 100644 --- a/sysdeps/ieee754/ldbl-96/s_sinl.c +++ b/sysdeps/ieee754/ldbl-96/s_sinl.c @@ -52,6 +52,7 @@ static char rcsid[] = "$NetBSD: $"; #include <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long double __sinl(long double x) { @@ -85,4 +86,4 @@ long double __sinl(long double x) } } } -weak_alias (__sinl, sinl) +libm_alias_ldouble (__sin, sin) diff --git a/sysdeps/ieee754/ldbl-96/s_tanhl.c b/sysdeps/ieee754/ldbl-96/s_tanhl.c index 38edf9f75e..b1b3e0637b 100644 --- a/sysdeps/ieee754/ldbl-96/s_tanhl.c +++ b/sysdeps/ieee754/ldbl-96/s_tanhl.c @@ -45,6 +45,8 @@ static char rcsid[] = "$NetBSD: $"; #include <float.h> #include <math.h> #include <math_private.h> +#include <math-underflow.h> +#include <libm-alias-ldouble.h> static const long double one=1.0, two=2.0, tiny = 1.0e-4900L; @@ -52,7 +54,7 @@ long double __tanhl(long double x) { long double t,z; int32_t se; - u_int32_t j0,j1,ix; + uint32_t j0,j1,ix; /* High word of |x|. */ GET_LDOUBLE_WORDS(se,j0,j1,x); @@ -87,4 +89,4 @@ long double __tanhl(long double x) } return (se&0x8000)? -z: z; } -weak_alias (__tanhl, tanhl) +libm_alias_ldouble (__tanh, tanh) diff --git a/sysdeps/ieee754/ldbl-96/s_tanl.c b/sysdeps/ieee754/ldbl-96/s_tanl.c index 3fbe4a8f6b..b4163792c5 100644 --- a/sysdeps/ieee754/ldbl-96/s_tanl.c +++ b/sysdeps/ieee754/ldbl-96/s_tanl.c @@ -51,6 +51,7 @@ static char rcsid[] = "$NetBSD: $"; #include <errno.h> #include <math.h> #include <math_private.h> +#include <libm-alias-ldouble.h> long double __tanl(long double x) { @@ -78,4 +79,4 @@ long double __tanl(long double x) -1 -- n odd */ } } -weak_alias (__tanl, tanl) +libm_alias_ldouble (__tan, tan) diff --git a/sysdeps/ieee754/ldbl-96/s_totalorderl.c b/sysdeps/ieee754/ldbl-96/s_totalorderl.c new file mode 100644 index 0000000000..f5be9bf28e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_totalorderl.c @@ -0,0 +1,59 @@ +/* Total order operation. ldbl-96 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 <float.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-ldouble.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalorderl (long double x, long double y) +{ + int16_t expx, expy; + uint32_t hx, hy; + uint32_t lx, ly; + GET_LDOUBLE_WORDS (expx, hx, lx, x); + GET_LDOUBLE_WORDS (expy, hy, ly, y); + if (LDBL_MIN_EXP == -16382) + { + /* M68K variant: for the greatest exponent, the high mantissa + bit is not significant and both values of it are valid, so + set it before comparing. For the Intel variant, only one + value of the high mantissa bit is valid for each exponent, so + this is not necessary. */ + if ((expx & 0x7fff) == 0x7fff) + hx |= 0x80000000; + if ((expy & 0x7fff) == 0x7fff) + hy |= 0x80000000; + } +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN +# error not implemented +#endif + uint32_t x_sign = expx >> 15; + uint32_t y_sign = expy >> 15; + expx ^= x_sign >> 17; + hx ^= x_sign; + lx ^= x_sign; + expy ^= y_sign >> 17; + hy ^= y_sign; + ly ^= y_sign; + return expx < expy || (expx == expy && (hx < hy || (hx == hy && lx <= ly))); +} +libm_alias_ldouble (__totalorder, totalorder) diff --git a/sysdeps/ieee754/ldbl-96/s_totalordermagl.c b/sysdeps/ieee754/ldbl-96/s_totalordermagl.c new file mode 100644 index 0000000000..18efefaee1 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_totalordermagl.c @@ -0,0 +1,53 @@ +/* Total order operation on absolute values. ldbl-96 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 <float.h> +#include <math.h> +#include <math_private.h> +#include <libm-alias-ldouble.h> +#include <nan-high-order-bit.h> +#include <stdint.h> + +int +__totalordermagl (long double x, long double y) +{ + uint16_t expx, expy; + uint32_t hx, hy; + uint32_t lx, ly; + GET_LDOUBLE_WORDS (expx, hx, lx, x); + GET_LDOUBLE_WORDS (expy, hy, ly, y); + expx &= 0x7fff; + expy &= 0x7fff; + if (LDBL_MIN_EXP == -16382) + { + /* M68K variant: for the greatest exponent, the high mantissa + bit is not significant and both values of it are valid, so + set it before comparing. For the Intel variant, only one + value of the high mantissa bit is valid for each exponent, so + this is not necessary. */ + if (expx == 0x7fff) + hx |= 0x80000000; + if (expy == 0x7fff) + hy |= 0x80000000; + } +#if HIGH_ORDER_BIT_IS_SET_FOR_SNAN +# error not implemented +#endif + return expx < expy || (expx == expy && (hx < hy || (hx == hy && lx <= ly))); +} +libm_alias_ldouble (__totalordermag, totalordermag) diff --git a/sysdeps/ieee754/ldbl-96/s_ufromfpl.c b/sysdeps/ieee754/ldbl-96/s_ufromfpl.c new file mode 100644 index 0000000000..22935e6ef7 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_ufromfpl.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 0 +#define FUNC __ufromfpl +#include <s_fromfpl_main.c> +libm_alias_ldouble (__ufromfp, ufromfp) diff --git a/sysdeps/ieee754/ldbl-96/s_ufromfpxl.c b/sysdeps/ieee754/ldbl-96/s_ufromfpxl.c new file mode 100644 index 0000000000..77a5423de8 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_ufromfpxl.c @@ -0,0 +1,5 @@ +#define UNSIGNED 1 +#define INEXACT 1 +#define FUNC __ufromfpxl +#include <s_fromfpl_main.c> +libm_alias_ldouble (__ufromfpx, ufromfpx) diff --git a/sysdeps/ieee754/ldbl-96/strtold_l.c b/sysdeps/ieee754/ldbl-96/strtold_l.c index b51560e18a..145e64f766 100644 --- a/sysdeps/ieee754/ldbl-96/strtold_l.c +++ b/sysdeps/ieee754/ldbl-96/strtold_l.c @@ -1,4 +1,4 @@ -/* Copyright (C) 1999-2016 Free Software Foundation, Inc. +/* Copyright (C) 1999-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 @@ -34,4 +34,19 @@ #define MPN2FLOAT __mpn_construct_long_double #define FLOAT_HUGE_VAL HUGE_VALL +#if __HAVE_FLOAT64X_LONG_DOUBLE +# define strtof64x_l __hide_strtof64x_l +# define wcstof64x_l __hide_wcstof64x_l +#endif + #include <stdlib/strtod_l.c> + +#if __HAVE_FLOAT64X_LONG_DOUBLE +# undef strtof64x_l +# undef wcstof64x_l +# ifdef USE_WIDE_CHAR +weak_alias (wcstold_l, wcstof64x_l) +# else +weak_alias (strtold_l, strtof64x_l) +# endif +#endif diff --git a/sysdeps/ieee754/ldbl-96/t_sincosl.c b/sysdeps/ieee754/ldbl-96/t_sincosl.c index f00a182061..5d3531f677 100644 --- a/sysdeps/ieee754/ldbl-96/t_sincosl.c +++ b/sysdeps/ieee754/ldbl-96/t_sincosl.c @@ -1,5 +1,5 @@ /* Extended-precision floating point sine and cosine tables. - Copyright (C) 1999-2016 Free Software Foundation, Inc. + Copyright (C) 1999-2018 Free Software Foundation, Inc. This file is part of the GNU C Library. Based on quad-precision tables by Jakub Jelinek <jj@ultra.linux.cz> diff --git a/sysdeps/ieee754/ldbl-96/test-canonical-ldbl-96.c b/sysdeps/ieee754/ldbl-96/test-canonical-ldbl-96.c new file mode 100644 index 0000000000..cd63e01140 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/test-canonical-ldbl-96.c @@ -0,0 +1,141 @@ +/* Test iscanonical and canonicalizel for ldbl-96. + 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_ldbl.h> +#include <stdbool.h> +#include <stdint.h> +#include <stdio.h> + +struct test +{ + bool sign; + uint16_t exponent; + bool high; + uint64_t mantissa; + bool canonical; +}; + +#define M68K_VARIANT (LDBL_MIN_EXP == -16382) + +static const struct test tests[] = + { + { false, 0, true, 0, M68K_VARIANT }, + { true, 0, true, 0, M68K_VARIANT }, + { false, 0, true, 1, M68K_VARIANT }, + { true, 0, true, 1, M68K_VARIANT }, + { false, 0, true, 0x100000000ULL, M68K_VARIANT }, + { true, 0, true, 0x100000000ULL, M68K_VARIANT }, + { false, 0, false, 0, true }, + { true, 0, false, 0, true }, + { false, 0, false, 1, true }, + { true, 0, false, 1, true }, + { false, 0, false, 0x100000000ULL, true }, + { true, 0, false, 0x100000000ULL, true }, + { false, 1, true, 0, true }, + { true, 1, true, 0, true }, + { false, 1, true, 1, true }, + { true, 1, true, 1, true }, + { false, 1, true, 0x100000000ULL, true }, + { true, 1, true, 0x100000000ULL, true }, + { false, 1, false, 0, false }, + { true, 1, false, 0, false }, + { false, 1, false, 1, false }, + { true, 1, false, 1, false }, + { false, 1, false, 0x100000000ULL, false }, + { true, 1, false, 0x100000000ULL, false }, + { false, 0x7ffe, true, 0, true }, + { true, 0x7ffe, true, 0, true }, + { false, 0x7ffe, true, 1, true }, + { true, 0x7ffe, true, 1, true }, + { false, 0x7ffe, true, 0x100000000ULL, true }, + { true, 0x7ffe, true, 0x100000000ULL, true }, + { false, 0x7ffe, false, 0, false }, + { true, 0x7ffe, false, 0, false }, + { false, 0x7ffe, false, 1, false }, + { true, 0x7ffe, false, 1, false }, + { false, 0x7ffe, false, 0x100000000ULL, false }, + { true, 0x7ffe, false, 0x100000000ULL, false }, + { false, 0x7fff, true, 0, true }, + { true, 0x7fff, true, 0, true }, + { false, 0x7fff, true, 1, true }, + { true, 0x7fff, true, 1, true }, + { false, 0x7fff, true, 0x100000000ULL, true }, + { true, 0x7fff, true, 0x100000000ULL, true }, + { false, 0x7fff, false, 0, M68K_VARIANT }, + { true, 0x7fff, false, 0, M68K_VARIANT }, + { false, 0x7fff, false, 1, M68K_VARIANT }, + { true, 0x7fff, false, 1, M68K_VARIANT }, + { false, 0x7fff, false, 0x100000000ULL, M68K_VARIANT }, + { true, 0x7fff, false, 0x100000000ULL, M68K_VARIANT }, + }; + +static int +do_test (void) +{ + int result = 0; + + for (size_t i = 0; i < sizeof (tests) / sizeof (tests[0]); i++) + { + long double ld; + SET_LDOUBLE_WORDS (ld, tests[i].exponent | (tests[i].sign << 15), + (tests[i].mantissa >> 32) | (tests[i].high << 31), + tests[i].mantissa & 0xffffffffULL); + bool canonical = iscanonical (ld); + if (canonical == tests[i].canonical) + { + printf ("PASS: iscanonical test %zu\n", i); + long double ldc = 12345.0L; + bool canonicalize_ret = canonicalizel (&ldc, &ld); + if (canonicalize_ret == !canonical) + { + printf ("PASS: canonicalizel test %zu\n", i); + bool canon_ok; + if (!canonical) + canon_ok = ldc == 12345.0L; + else if (isnan (ld)) + canon_ok = isnan (ldc) && !issignaling (ldc); + else + canon_ok = ldc == ld; + if (canon_ok) + printf ("PASS: canonicalized value test %zu\n", i); + else + { + printf ("FAIL: canonicalized value test %zu\n", i); + result = 1; + } + } + else + { + printf ("FAIL: canonicalizel test %zu\n", i); + result = 1; + } + } + else + { + printf ("FAIL: iscanonical test %zu\n", i); + result = 1; + } + } + + return result; +} + +#define TEST_FUNCTION do_test () +#include "../test-skeleton.c" diff --git a/sysdeps/ieee754/ldbl-96/test-totalorderl-ldbl-96.c b/sysdeps/ieee754/ldbl-96/test-totalorderl-ldbl-96.c new file mode 100644 index 0000000000..552a97783a --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/test-totalorderl-ldbl-96.c @@ -0,0 +1,82 @@ +/* Test totalorderl and totalordermagl for ldbl-96. + 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_ldbl.h> +#include <stdbool.h> +#include <stdint.h> +#include <stdio.h> + +static const uint64_t tests[] = + { + 0, 1, 0x4000000000000000ULL, 0x4000000000000001ULL, + 0x7fffffffffffffffULL + }; + +static int +do_test (void) +{ + int result = 0; + + if (LDBL_MIN_EXP == -16382) + for (size_t i = 0; i < sizeof (tests) / sizeof (tests[0]); i++) + { + long double ldx, ldy, ldnx, ldny; + /* Verify that the high bit of the mantissa is ignored for + infinities and NaNs for the M68K variant of this + format. */ + SET_LDOUBLE_WORDS (ldx, 0x7fff, + tests[i] >> 32, tests[i] & 0xffffffffULL); + SET_LDOUBLE_WORDS (ldy, 0x7fff, + (tests[i] >> 32) | 0x80000000, + tests[i] & 0xffffffffULL); + SET_LDOUBLE_WORDS (ldnx, 0xffff, + tests[i] >> 32, tests[i] & 0xffffffffULL); + SET_LDOUBLE_WORDS (ldny, 0xffff, + (tests[i] >> 32) | 0x80000000, + tests[i] & 0xffffffffULL); + bool to1 = totalorderl (ldx, ldy); + bool to2 = totalorderl (ldy, ldx); + bool to3 = totalorderl (ldnx, ldny); + bool to4 = totalorderl (ldny, ldnx); + if (to1 && to2 && to3 && to4) + printf ("PASS: test %zu\n", i); + else + { + printf ("FAIL: test %zu\n", i); + result = 1; + } + to1 = totalordermagl (ldx, ldy); + to2 = totalordermagl (ldy, ldx); + to3 = totalordermagl (ldnx, ldny); + to4 = totalordermagl (ldny, ldnx); + if (to1 && to2 && to3 && to4) + printf ("PASS: test %zu (totalordermagl)\n", i); + else + { + printf ("FAIL: test %zu (totalordermagl)\n", i); + result = 1; + } + } + + return result; +} + +#define TEST_FUNCTION do_test () +#include "../test-skeleton.c" diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1.c b/sysdeps/ieee754/ldbl-96/x2y2m1.c index 4119f42acc..afe7ab161b 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1.c +++ b/sysdeps/ieee754/ldbl-96/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 diff --git a/sysdeps/ieee754/ldbl-96/x2y2m1l.c b/sysdeps/ieee754/ldbl-96/x2y2m1l.c index 733742da04..392830c1b0 100644 --- a/sysdeps/ieee754/ldbl-96/x2y2m1l.c +++ b/sysdeps/ieee754/ldbl-96/x2y2m1l.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_splitl.h> #include <stdlib.h> /* Calculate X + Y exactly and store the result in *HI + *LO. It is @@ -33,36 +33,6 @@ add_split (long double *hi, long double *lo, long double x, long 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 inline void -mul_split (long double *hi, long double *lo, long double x, long double y) -{ -#ifdef __FP_FAST_FMAL - /* Fast built-in fused multiply-add. */ - *hi = x * y; - *lo = __builtin_fmal (x, y, -*hi); -#elif defined FP_FAST_FMAL - /* Fast library fused multiply-add, compiler before GCC 4.6. */ - *hi = x * y; - *lo = __fmal (x, y, -*hi); -#else - /* Apply Dekker's algorithm. */ - *hi = x * y; -# define C ((1LL << (LDBL_MANT_DIG + 1) / 2) + 1) - long double x1 = x * C; - long double y1 = y * C; -# undef C - x1 = (x - x1) + x1; - y1 = (y - y1) + y1; - long double x2 = x - x1; - long 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. */ @@ -88,8 +58,8 @@ __x2y2m1l (long double x, long double y) { long double vals[5]; SET_RESTORE_ROUNDL (FE_TONEAREST); - mul_split (&vals[1], &vals[0], x, x); - mul_split (&vals[3], &vals[2], y, y); + mul_splitl (&vals[1], &vals[0], x, x); + mul_splitl (&vals[3], &vals[2], y, y); vals[4] = -1.0L; qsort (vals, 5, sizeof (long double), compare); /* Add up the values so that each element of VALS has absolute value |