Issue #1291 - Part 3: Update fdlibm to Sept 2019 version

13 files changed, 144 insertions, 89 deletions

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 <cmath>
 
+#include <float.h>
 #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 <das@FreeBSD.ORG>
  * 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 <sys/cdefs.h>
 //__FBSDID("$FreeBSD$");
 
+#include <float.h>
 #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 <sys/cdefs.h>
+//__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 <das@FreeBSD.ORG>
  * 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 <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
 
 /*
  * scalbn (double x, int n)
@@ -21,7 +20,6 @@
  * exponentiation or a multiplication.
  */
 
-//#include <sys/cdefs.h>
 #include <float.h>
 
 #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;