From 5e902a4a9a25a9e151b5d3945c68517bc0c73139 Mon Sep 17 00:00:00 2001 From: wolfbeast Date: Fri, 15 Nov 2019 15:17:01 +0100 Subject: Issue #1291 - Part 3: Update fdlibm to Sept 2019 version --- modules/fdlibm/src/e_asin.cpp | 2 +- modules/fdlibm/src/e_atan2.cpp | 4 +- modules/fdlibm/src/e_exp.cpp | 4 +- modules/fdlibm/src/e_hypot.cpp | 2 +- modules/fdlibm/src/e_pow.cpp | 50 ++++++++------- modules/fdlibm/src/k_exp.cpp | 2 + modules/fdlibm/src/math_private.h | 128 ++++++++++++++++++++++++------------- modules/fdlibm/src/moz.build | 17 +++-- modules/fdlibm/src/s_cbrt.cpp | 1 + modules/fdlibm/src/s_expm1.cpp | 2 +- modules/fdlibm/src/s_fabs.cpp | 7 +- modules/fdlibm/src/s_nearbyint.cpp | 2 + modules/fdlibm/src/s_scalbn.cpp | 12 ++-- 13 files changed, 144 insertions(+), 89 deletions(-) (limited to 'modules/fdlibm/src') diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp index 98607937d..e896bde9e 100644 --- a/modules/fdlibm/src/e_asin.cpp +++ b/modules/fdlibm/src/e_asin.cpp @@ -6,7 +6,7 @@ * * Developed at SunSoft, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp index 9990072cf..f45ad187f 100644 --- a/modules/fdlibm/src/e_atan2.cpp +++ b/modules/fdlibm/src/e_atan2.cpp @@ -69,8 +69,8 @@ __ieee754_atan2(double y, double x) iy = hy&0x7fffffff; if(((ix|((lx|-lx)>>31))>0x7ff00000)|| ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */ - return x+y; - if((hx-0x3ff00000|lx)==0) return atan(y); /* x=1.0 */ + return nan_mix(x, y); + if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */ m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ /* when y = 0 */ diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp index c4e530d4d..92af819ce 100644 --- a/modules/fdlibm/src/e_exp.cpp +++ b/modules/fdlibm/src/e_exp.cpp @@ -147,9 +147,9 @@ __ieee754_exp(double x) /* default IEEE double exp */ /* x is now in primary range */ t = x*x; if(k >= -1021) - INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0); + INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0); else - INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0); + INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0); c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); if(k==0) return one-((x*c)/(c-2.0)-x); else y = one-((lo-(x*c)/(2.0-c))-hi); diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp index b34a1f575..a23571150 100644 --- a/modules/fdlibm/src/e_hypot.cpp +++ b/modules/fdlibm/src/e_hypot.cpp @@ -70,7 +70,7 @@ __ieee754_hypot(double x, double y) if(ha >= 0x7ff00000) { /* Inf or NaN */ u_int32_t low; /* Use original arg order iff result is NaN; quieten sNaNs. */ - w = fabs(x+0.0)-fabs(y+0.0); + w = fabsl(x+0.0L)-fabs(y+0); GET_LOW_WORD(low,a); if(((ha&0xfffff)|low)==0) w = a; GET_LOW_WORD(low,b); diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp index 7f22afce8..c18226b8a 100644 --- a/modules/fdlibm/src/e_pow.cpp +++ b/modules/fdlibm/src/e_pow.cpp @@ -4,7 +4,7 @@ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. * * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -19,7 +19,7 @@ * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. - * 2. Perform y*log2(x) = n+y' by simulating multi-precision + * 2. Perform y*log2(x) = n+y' by simulating multi-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * @@ -47,18 +47,19 @@ * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) - * always returns the correct integer provided it is + * always returns the correct integer provided it is * representable. * * Constants : - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the - * compiler will convert from decimal to binary accurately enough + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #include +#include #include "math_private.h" static const double @@ -66,6 +67,9 @@ bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ zero = 0.0, +half = 0.5, +qrtr = 0.25, +thrd = 3.3333333333333331e-01, /* 0x3fd55555, 0x55555555 */ one = 1.0, two = 2.0, two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ @@ -108,15 +112,15 @@ __ieee754_pow(double x, double y) ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ - if((iy|ly)==0) return one; + if((iy|ly)==0) return one; /* x==1: 1**y = 1, even if y is NaN */ if (hx==0x3ff00000 && lx == 0) return one; /* y!=zero: result is NaN if either arg is NaN */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || - iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return (x+0.0)+(y+0.0); + iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) + return nan_mix(x, y); /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -124,22 +128,22 @@ __ieee754_pow(double x, double y) * yisint = 2 ... y is an even int */ yisint = 0; - if(hx<0) { + if(hx<0) { if(iy>=0x43400000) yisint = 2; /* even integer y */ else if(iy>=0x3ff00000) { k = (iy>>20)-0x3ff; /* exponent */ if(k>20) { j = ly>>(52-k); - if((j<<(52-k))==ly) yisint = 2-(j&1); + if(((u_int32_t)j<<(52-k))==ly) yisint = 2-(j&1); } else if(ly==0) { j = iy>>(20-k); if((j<<(20-k))==iy) yisint = 2-(j&1); } - } - } + } + } /* special value of y */ - if(ly==0) { + if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) return one; /* (-1)**+-inf is 1 */ @@ -147,7 +151,7 @@ __ieee754_pow(double x, double y) return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero; - } + } if(iy==0x3ff00000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } @@ -167,13 +171,13 @@ __ieee754_pow(double x, double y) if(hx<0) { if(((ix-0x3ff00000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) + } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } - + /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be n = (hx>>31)+1; but ANSI C says a right shift of a signed negative quantity is @@ -195,10 +199,10 @@ __ieee754_pow(double x, double y) /* over/underflow if x is not close to one */ if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny; if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute + /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax-one; /* t has 20 trailing zeros */ - w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); + w = (t*t)*(half-t*(thrd-t*qrtr)); u = ivln2_h*t; /* ivln2_h has 21 sig. bits */ v = t*ivln2_l-w*ivln2; t1 = u+v; @@ -235,9 +239,9 @@ __ieee754_pow(double x, double y) r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6))))); r += s_l*(s_h+ss); s2 = s_h*s_h; - t_h = 3.0+s2+r; + t_h = 3+s2+r; SET_LOW_WORD(t_h,0); - t_l = r-((t_h-3.0)-s2); + t_l = r-((t_h-3)-s2); /* u+v = ss*(1+...) */ u = s_h*t_h; v = s_l*t_h+t_l*ss; @@ -248,7 +252,7 @@ __ieee754_pow(double x, double y) z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */ z_l = cp_l*p_h+p_l*cp+dp_l[k]; /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */ - t = (double)n; + t = n; t1 = (((z_h+z_l)+dp_h[k])+t); SET_LOW_WORD(t1,0); t2 = z_l-(((t1-t)-dp_h[k])-z_h); @@ -288,7 +292,7 @@ __ieee754_pow(double x, double y) n = ((n&0x000fffff)|0x00100000)>>(20-k); if(j<0) n = -n; p_h -= t; - } + } t = p_l+p_h; SET_LOW_WORD(t,0); u = t*lg2_h; diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp index a0699fa4a..9394c8fd8 100644 --- a/modules/fdlibm/src/k_exp.cpp +++ b/modules/fdlibm/src/k_exp.cpp @@ -1,4 +1,6 @@ /*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * * Copyright (c) 2011 David Schultz * All rights reserved. * diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h index 3f7b03291..d9ec44817 100644 --- a/modules/fdlibm/src/math_private.h +++ b/modules/fdlibm/src/math_private.h @@ -45,6 +45,47 @@ #define u_int64_t uint64_t #endif +/* A union which permits us to convert between a long double and + four 32 bit ints. */ + +#if MOZ_BIG_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t mswhi; + u_int32_t mswlo; + u_int32_t lswhi; + u_int32_t lswlo; + } parts32; + struct { + u_int64_t msw; + u_int64_t lsw; + } parts64; +} ieee_quad_shape_type; + +#endif + +#if MOZ_LITTLE_ENDIAN + +typedef union +{ + long double value; + struct { + u_int32_t lswlo; + u_int32_t lswhi; + u_int32_t mswlo; + u_int32_t mswhi; + } parts32; + struct { + u_int64_t lsw; + u_int64_t msw; + } parts64; +} ieee_quad_shape_type; + +#endif + /* * A union which permits us to convert between a double and two 32 bit * ints. @@ -307,8 +348,9 @@ do { \ /* Support switching the mode to FP_PE if necessary. */ #if defined(__i386__) && !defined(NO_FPSETPREC) -#define ENTERI() \ - long double __retval; \ +#define ENTERI() ENTERIT(long double) +#define ENTERIT(returntype) \ + returntype __retval; \ fp_prec_t __oprec; \ \ if ((__oprec = fpgetprec()) != FP_PE) \ @@ -319,9 +361,22 @@ do { \ fpsetprec(__oprec); \ RETURNF(__retval); \ } while (0) +#define ENTERV() \ + fp_prec_t __oprec; \ + \ + if ((__oprec = fpgetprec()) != FP_PE) \ + fpsetprec(FP_PE) +#define RETURNV() do { \ + if (__oprec != FP_PE) \ + fpsetprec(__oprec); \ + return; \ +} while (0) #else -#define ENTERI(x) +#define ENTERI() +#define ENTERIT(x) #define RETURNI(x) RETURNF(x) +#define ENTERV() +#define RETURNV() return #endif /* Default return statement if hack*_t() is not used. */ @@ -436,6 +491,31 @@ do { \ */ void _scan_nan(uint32_t *__words, int __num_words, const char *__s); +/* + * Mix 0, 1 or 2 NaNs. First add 0 to each arg. This normally just turns + * signaling NaNs into quiet NaNs by setting a quiet bit. We do this + * because we want to never return a signaling NaN, and also because we + * don't want the quiet bit to affect the result. Then mix the converted + * args using the specified operation. + * + * When one arg is NaN, the result is typically that arg quieted. When both + * args are NaNs, the result is typically the quietening of the arg whose + * mantissa is largest after quietening. When neither arg is NaN, the + * result may be NaN because it is indeterminate, or finite for subsequent + * construction of a NaN as the indeterminate 0.0L/0.0L. + * + * Technical complications: the result in bits after rounding to the final + * precision might depend on the runtime precision and/or on compiler + * optimizations, especially when different register sets are used for + * different precisions. Try to make the result not depend on at least the + * runtime precision by always doing the main mixing step in long double + * precision. Try to reduce dependencies on optimizations by adding the + * the 0's in different precisions (unless everything is in long double + * precision). + */ +#define nan_mix(x, y) (nan_mix_op((x), (y), +)) +#define nan_mix_op(x, y, op) (((x) + 0.0L) op ((y) + 0)) + #ifdef _COMPLEX_H /* @@ -511,48 +591,6 @@ CMPLXL(long double x, long double y) #endif /* _COMPLEX_H */ -#ifdef __GNUCLIKE_ASM - -/* Asm versions of some functions. */ - -#ifdef __amd64__ -static __inline int -irint(double x) -{ - int n; - - asm("cvtsd2si %1,%0" : "=r" (n) : "x" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#ifdef __i386__ -static __inline int -irint(double x) -{ - int n; - - asm("fistl %0" : "=m" (n) : "t" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINT -#endif - -#if defined(__amd64__) || defined(__i386__) -static __inline int -irintl(long double x) -{ - int n; - - asm("fistl %0" : "=m" (n) : "t" (x)); - return (n); -} -#define HAVE_EFFICIENT_IRINTL -#endif - -#endif /* __GNUCLIKE_ASM */ - #ifdef DEBUG #if defined(__amd64__) || defined(__i386__) #define breakpoint() asm("int $3") diff --git a/modules/fdlibm/src/moz.build b/modules/fdlibm/src/moz.build index cae0e5bc9..be5bf3d9b 100644 --- a/modules/fdlibm/src/moz.build +++ b/modules/fdlibm/src/moz.build @@ -10,26 +10,35 @@ EXPORTS += [ FINAL_LIBRARY = 'js' -if CONFIG['GNU_CXX']: +if CONFIG['CC_TYPE'] in ('clang', 'gcc'): CXXFLAGS += [ '-Wno-parentheses', '-Wno-sign-compare', ] -if CONFIG['CLANG_CXX']: +if CONFIG['CC_TYPE'] == 'clang': CXXFLAGS += [ '-Wno-dangling-else', ] -if CONFIG['_MSC_VER']: +if CONFIG['CC_TYPE'] in ('msvc', 'clang-cl'): CXXFLAGS += [ - '-wd4018', # signed/unsigned mismatch '-wd4146', # unary minus operator applied to unsigned type '-wd4305', # truncation from 'double' to 'const float' '-wd4723', # potential divide by 0 '-wd4756', # overflow in constant arithmetic ] +if CONFIG['CC_TYPE'] == 'msvc': + CXXFLAGS += [ + '-wd4018', # signed/unsigned mismatch + ] + +if CONFIG['CC_TYPE'] == 'clang-cl': + CXXFLAGS += [ + '-Wno-sign-compare', # signed/unsigned mismatch + ] + SOURCES += [ 'e_acos.cpp', 'e_acosh.cpp', diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp index a2de24af7..fe3747e81 100644 --- a/modules/fdlibm/src/s_cbrt.cpp +++ b/modules/fdlibm/src/s_cbrt.cpp @@ -15,6 +15,7 @@ //#include //__FBSDID("$FreeBSD$"); +#include #include "math_private.h" /* cbrt(x) diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp index 4c19485de..90ebc1698 100644 --- a/modules/fdlibm/src/s_expm1.cpp +++ b/modules/fdlibm/src/s_expm1.cpp @@ -187,7 +187,7 @@ expm1(double x) e = hxs*((r1-t)/(6.0 - x*t)); if(k==0) return x - (x*e-hxs); /* c is 0 */ else { - INSERT_WORDS(twopk,0x3ff00000+(k<<20),0); /* 2^k */ + INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20,0); /* 2^k */ e = (x*(e-c)-c); e -= hxs; if(k== -1) return 0.5*(x-e)-0.5; diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp index 3bea0478a..6ca84d71b 100644 --- a/modules/fdlibm/src/s_fabs.cpp +++ b/modules/fdlibm/src/s_fabs.cpp @@ -1,4 +1,4 @@ - /* @(#)s_fabs.c 5.1 93/09/24 */ +/* @(#)s_fabs.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint - //static char rcsid[] = "$FreeBSD$"; -#endif +//#include +//__FBSDID("$FreeBSD$"); /* * fabs(x) returns the absolute value of x. diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp index 532bb5d8d..6c04212d3 100644 --- a/modules/fdlibm/src/s_nearbyint.cpp +++ b/modules/fdlibm/src/s_nearbyint.cpp @@ -1,4 +1,6 @@ /*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD + * * Copyright (c) 2004 David Schultz * All rights reserved. * diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp index 5dbf58c23..dfbcf5c57 100644 --- a/modules/fdlibm/src/s_scalbn.cpp +++ b/modules/fdlibm/src/s_scalbn.cpp @@ -10,9 +10,8 @@ * ==================================================== */ -#ifndef lint -//static char rcsid[] = "$FreeBSD$"; -#endif +//#include +//__FBSDID("$FreeBSD$"); /* * scalbn (double x, int n) @@ -21,7 +20,6 @@ * exponentiation or a multiplication. */ -//#include #include #include "math_private.h" @@ -50,10 +48,12 @@ scalbn (double x, int n) if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ if (k > 0) /* normal result */ {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;} - if (k <= -54) + if (k <= -54) { if (n > 50000) /* in case integer overflow in n+k */ return huge*copysign(huge,x); /*overflow*/ - else return tiny*copysign(tiny,x); /*underflow*/ + else + return tiny*copysign(tiny,x); /*underflow*/ + } k += 54; /* subnormal result */ SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x*twom54; -- cgit v1.2.3