Bug 1583958 - Update our fdlibm import to upstream rev cf4707bb2f78ecf56ba350bdc24e3135b4339622, adjusting local patchwork as necessary for changes upstream's made since last import. r=arai
authorJeff Walden <jwalden@mit.edu>
Thu, 26 Sep 2019 07:07:16 +0000
changeset 495057 e1df5d37a8dd9092f11c36d7d76b04c3591b5f89
parent 495056 36c6fe0b96ab5c2df25daf480c58ce97d2084d3e
child 495058 2769fdb65549c8678e17ab14bc763fb0866d3519
push id36622
push usershindli@mozilla.com
push dateThu, 26 Sep 2019 21:35:42 +0000
treeherdermozilla-central@1d189ae70326 [default view] [failures only]
perfherder[talos] [build metrics] [platform microbench] (compared to previous push)
reviewersarai
bugs1583958
milestone71.0a1
first release with
nightly linux32
nightly linux64
nightly mac
nightly win32
nightly win64
last release without
nightly linux32
nightly linux64
nightly mac
nightly win32
nightly win64
Bug 1583958 - Update our fdlibm import to upstream rev cf4707bb2f78ecf56ba350bdc24e3135b4339622, adjusting local patchwork as necessary for changes upstream's made since last import. r=arai Differential Revision: https://phabricator.services.mozilla.com/D47186
modules/fdlibm/README.mozilla
modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch
modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch
modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch
modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch
modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch
modules/fdlibm/patches/08_remove_weak_reference_macro.patch
modules/fdlibm/patches/09_comment_out_rcsid_variable.patch
modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch
modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch
modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch
modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch
modules/fdlibm/patches/18_use_stdlib_sqrt.patch
modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch
modules/fdlibm/src/e_atan2.cpp
modules/fdlibm/src/e_exp.cpp
modules/fdlibm/src/e_hypot.cpp
modules/fdlibm/src/e_pow.cpp
modules/fdlibm/src/math_private.h
modules/fdlibm/src/s_cbrt.cpp
modules/fdlibm/src/s_expm1.cpp
--- a/modules/fdlibm/README.mozilla
+++ b/modules/fdlibm/README.mozilla
@@ -8,12 +8,12 @@ Each file is downloaded separately, as c
 resources.
 
 The in-tree copy is updated by running
   sh update.sh
 or
   sh update.sh <sha-commit>
 from within the modules/fdlibm directory.
 
-Current version: [commit b21ccf63f28a3a4692d8a31419e0a725a1b1a800 (2018-02-14T07:59:30Z)].
+Current version: [commit cf4707bb2f78ecf56ba350bdc24e3135b4339622 (2019-09-25T18:50:57Z)].
 
 patches 01-18 fixes files to be usable within mozilla-central tree.
 See https://bugzilla.mozilla.org/show_bug.cgi?id=933257
--- a/modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch
+++ b/modules/fdlibm/patches/02_change_include_guard_in_fdlibm_h.patch
@@ -17,17 +17,17 @@ diff --git a/modules/fdlibm/src/fdlibm.h
  
  double	acos(double);
  double	asin(double);
  double	atan(double);
  double	atan2(double, double);
  
  double	cosh(double);
  double	sinh(double);
-@@ -53,9 +53,9 @@ double	scalbn(double, int);
+@@ -52,9 +52,9 @@ double	scalbn(double, int);
  
  float	ceilf(float);
  float	floorf(float);
  
  float	nearbyintf(float);
  float	rintf(float);
  float	truncf(float);
  
--- a/modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch
+++ b/modules/fdlibm/patches/03_put_fdlibm_functions_into_fdlibm_namespace.patch
@@ -15,17 +15,17 @@ diff --git a/modules/fdlibm/src/fdlibm.h
  double	acos(double);
  double	asin(double);
  double	atan(double);
  double	atan2(double, double);
  
  double	cosh(double);
  double	sinh(double);
  double	tanh(double);
-@@ -53,9 +55,11 @@ double	scalbn(double, int);
+@@ -52,9 +54,11 @@ double	scalbn(double, int);
  
  float	ceilf(float);
  float	floorf(float);
  
  float	nearbyintf(float);
  float	rintf(float);
  float	truncf(float);
  
--- a/modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch
+++ b/modules/fdlibm/patches/04_include_fdlibm_h_from_math_private_h.patch
@@ -227,34 +227,34 @@ diff --git a/modules/fdlibm/src/e_log2.c
  static const double
  two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
  ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
  
 diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp
 --- a/modules/fdlibm/src/e_pow.cpp
 +++ b/modules/fdlibm/src/e_pow.cpp
-@@ -52,17 +52,16 @@
-  *
+@@ -53,17 +53,16 @@
   * 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
   * to produce the hexadecimal values shown.
   */
  
+ #include <float.h>
 -#include "math.h"
  #include "math_private.h"
  
  static const double
  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,
- one	=  1.0,
+ half    =  0.5,
 diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
 --- a/modules/fdlibm/src/e_sinh.cpp
 +++ b/modules/fdlibm/src/e_sinh.cpp
 @@ -29,17 +29,16 @@
   *
   * Special cases:
   *	sinh(x) is |x| if x is +INF, -INF, or NaN.
   *	only sinh(0)=0 is exact for finite x.
@@ -354,25 +354,25 @@ diff --git a/modules/fdlibm/src/s_atan.c
    4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
    7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
    9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
    1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
  };
 diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp
 --- a/modules/fdlibm/src/s_cbrt.cpp
 +++ b/modules/fdlibm/src/s_cbrt.cpp
-@@ -10,17 +10,16 @@
-  * ====================================================
+@@ -11,17 +11,16 @@
   *
   * Optimized by Bruce D. Evans.
   */
  
  #include <sys/cdefs.h>
  __FBSDID("$FreeBSD$");
  
+ #include <float.h>
 -#include "math.h"
  #include "math_private.h"
  
  /* cbrt(x)
   * Return cube root of x
   */
  static const u_int32_t
  	B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
--- a/modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch
+++ b/modules/fdlibm/patches/06_use_mfbt_endian_h_in_math_private_h.patch
@@ -1,12 +1,12 @@
 diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
 --- a/modules/fdlibm/src/math_private.h
 +++ b/modules/fdlibm/src/math_private.h
-@@ -14,20 +14,21 @@
+@@ -14,52 +14,38 @@
   * $FreeBSD$
   */
  
  #ifndef _MATH_PRIVATE_H_
  #define	_MATH_PRIVATE_H_
  
  #include <stdint.h>
  #include <sys/types.h>
@@ -19,46 +19,93 @@ diff --git a/modules/fdlibm/src/math_pri
  /*
   * The original fdlibm code used statements like:
   *	n0 = ((*(int*)&one)>>29)^1;		* index of high word *
   *	ix0 = *(n0+(int*)&x);			* high word of x *
   *	ix1 = *((1-n0)+(int*)&x);		* low word of x *
   * to dig two 32 bit words out of the 64 bit IEEE floating point
   * value.  That is non-ANSI, and, moreover, the gcc instruction
   * scheduler gets it wrong.  We instead use the following macros.
-@@ -36,27 +37,17 @@
+  * Unlike the original code, we determine the endianness at compile
+  * time, not at run time; I don't see much benefit to selecting
   * endianness at run time.
   */
  
- /*
-  * A union which permits us to convert between a double and two 32 bit
-  * ints.
-  */
- 
+-/*
+- * A union which permits us to convert between a double and two 32 bit
+- * ints.
+- */
+-
 -#ifdef __arm__
 -#if defined(__VFP_FP__) || defined(__ARM_EABI__)
 -#define	IEEE_WORD_ORDER	BYTE_ORDER
 -#else
 -#define	IEEE_WORD_ORDER	BIG_ENDIAN
 -#endif
 -#else /* __arm__ */
 -#define	IEEE_WORD_ORDER	BYTE_ORDER
 -#endif
 -
+ /* A union which permits us to convert between a long double and
+    four 32 bit ints.  */
+ 
 -#if IEEE_WORD_ORDER == BIG_ENDIAN
 +#if MOZ_BIG_ENDIAN
  
  typedef union
  {
+   long double value;
+   struct {
+     u_int32_t mswhi;
+     u_int32_t mswlo;
+     u_int32_t lswhi;
+@@ -68,17 +54,17 @@ typedef union
+   struct {
+     u_int64_t msw;
+     u_int64_t lsw;
+   } parts64;
+ } ieee_quad_shape_type;
+ 
+ #endif
+ 
+-#if IEEE_WORD_ORDER == LITTLE_ENDIAN
++#if MOZ_LITTLE_ENDIAN
+ 
+ typedef union
+ {
+   long double value;
+   struct {
+     u_int32_t lswlo;
+     u_int32_t lswhi;
+     u_int32_t mswlo;
+@@ -87,17 +73,22 @@ typedef union
+   struct {
+     u_int64_t lsw;
+     u_int64_t msw;
+   } parts64;
+ } ieee_quad_shape_type;
+ 
+ #endif
+ 
+-#if IEEE_WORD_ORDER == BIG_ENDIAN
++/*
++ * A union which permits us to convert between a double and two 32 bit
++ * ints.
++ */
++
++#if MOZ_BIG_ENDIAN
+ 
+ typedef union
+ {
    double value;
    struct
    {
      u_int32_t msw;
      u_int32_t lsw;
-@@ -64,17 +55,17 @@ typedef union
+@@ -105,17 +96,17 @@ typedef union
    struct
    {
      u_int64_t w;
    } xparts;
  } ieee_double_shape_type;
  
  #endif
  
--- a/modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch
+++ b/modules/fdlibm/patches/07_add_fdlibm_namespace_to_functions_defined_and_used_in_fdlibm.patch
@@ -1,12 +1,12 @@
 diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
 --- a/modules/fdlibm/src/math_private.h
 +++ b/modules/fdlibm/src/math_private.h
-@@ -742,16 +742,50 @@ irintl(long double x)
+@@ -872,16 +872,50 @@ irintl(long double x)
  #define	__ieee754_j1f	j1f
  #define	__ieee754_y0f	y0f
  #define	__ieee754_y1f	y1f
  #define	__ieee754_jnf	jnf
  #define	__ieee754_ynf	ynf
  #define	__ieee754_remainderf remainderf
  #define	__ieee754_scalbf scalbf
  
--- a/modules/fdlibm/patches/08_remove_weak_reference_macro.patch
+++ b/modules/fdlibm/patches/08_remove_weak_reference_macro.patch
@@ -169,16 +169,32 @@ diff --git a/modules/fdlibm/src/e_log2.c
  	val_hi = w;
  
  	return val_lo + val_hi;
  }
 -
 -#if (LDBL_MANT_DIG == 53)
 -__weak_reference(log2, log2l);
 -#endif
+diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp
+--- a/modules/fdlibm/src/e_pow.cpp
++++ b/modules/fdlibm/src/e_pow.cpp
+@@ -302,12 +302,8 @@ double
+ 	r  = (z*t1)/(t1-two)-(w+z*w);
+ 	z  = one-(r-z);
+ 	GET_HIGH_WORD(j,z);
+ 	j += (n<<20);
+ 	if((j>>20)<=0) z = scalbn(z,n);	/* subnormal output */
+ 	else SET_HIGH_WORD(z,j);
+ 	return s*z;
+ }
+-
+-#if (LDBL_MANT_DIG == 53)
+-__weak_reference(pow, powl);
+-#endif
 diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
 --- a/modules/fdlibm/src/e_sinh.cpp
 +++ b/modules/fdlibm/src/e_sinh.cpp
 @@ -67,12 +67,8 @@ double
  
      /* |x| in [log(maxdouble), overflowthresold] */
  	if (ix<=0x408633CE)
  	    return h*2.0*__ldexp_exp(fabs(x), -1);
@@ -220,17 +236,17 @@ diff --git a/modules/fdlibm/src/s_atan.c
  }
 -
 -#if LDBL_MANT_DIG == 53
 -__weak_reference(atan, atanl);
 -#endif
 diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp
 --- a/modules/fdlibm/src/s_cbrt.cpp
 +++ b/modules/fdlibm/src/s_cbrt.cpp
-@@ -105,12 +105,8 @@ cbrt(double x)
+@@ -106,12 +106,8 @@ cbrt(double x)
  	s=t*t;				/* t*t is exact */
  	r=x/s;				/* error <= 0.5 ulps; |r| < |t| */
  	w=t+t;				/* t+t is exact */
  	r=(r-t)/(w+r);			/* r-t is exact; w+r ~= 3*t */
  	t=t+t*r;			/* error <= 0.5 + 0.5/3 + epsilon */
  
  	return(t);
  }
--- a/modules/fdlibm/patches/09_comment_out_rcsid_variable.patch
+++ b/modules/fdlibm/patches/09_comment_out_rcsid_variable.patch
@@ -438,23 +438,23 @@ diff --git a/modules/fdlibm/src/s_cbrt.c
   * Optimized by Bruce D. Evans.
   */
  
 -#include <sys/cdefs.h>
 -__FBSDID("$FreeBSD$");
 +//#include <sys/cdefs.h>
 +//__FBSDID("$FreeBSD$");
  
+ #include <float.h>
  #include "math_private.h"
  
  /* cbrt(x)
   * Return cube root of x
   */
  static const u_int32_t
- 	B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
 diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp
 --- a/modules/fdlibm/src/s_ceil.cpp
 +++ b/modules/fdlibm/src/s_ceil.cpp
 @@ -5,18 +5,18 @@
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice
--- a/modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch
+++ b/modules/fdlibm/patches/12_define_u_int32_t_and_u_int64_t_on_windows.patch
@@ -12,16 +12,16 @@ diff --git a/modules/fdlibm/src/math_pri
  
 +#ifndef u_int32_t
 +#define u_int32_t uint32_t
 +#endif
 +#ifndef u_int64_t
 +#define u_int64_t uint64_t
 +#endif
 +
- /*
-  * A union which permits us to convert between a double and two 32 bit
-  * ints.
-  */
+ /* 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;
--- a/modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch
+++ b/modules/fdlibm/patches/13_define_strict_assign_even_if_flt_eval_method_is_not_defined.patch
@@ -1,12 +1,12 @@
 diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
 --- a/modules/fdlibm/src/math_private.h
 +++ b/modules/fdlibm/src/math_private.h
-@@ -287,16 +287,27 @@ do {								\
+@@ -328,16 +328,27 @@ do {								\
  	if (sizeof(type) >= sizeof(long double))	\
  		(lval) = (rval);		\
  	else {					\
  		__lval = (rval);		\
  		(lval) = __lval;		\
  	}					\
  } while (0)
  #endif
@@ -20,12 +20,12 @@ diff --git a/modules/fdlibm/src/math_pri
 +		__lval = (rval);		\
 +		(lval) = __lval;		\
 +	}					\
 +} while (0)
  #endif /* FLT_EVAL_METHOD */
  
  /* 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;			\
- 						\
--- a/modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch
+++ b/modules/fdlibm/patches/14_do_not_use_hexadecimal_floating_point_number.patch
@@ -1,16 +1,16 @@
 diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp
 --- a/modules/fdlibm/src/e_exp.cpp
 +++ b/modules/fdlibm/src/e_exp.cpp
 @@ -146,14 +146,17 @@ double
  	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);
  	if(k >= -1021) {
 -	    if (k==1024) return y*2.0*0x1p1023;
 +	    if (k==1024) {
 +	        double const_0x1p1023 = pow(2, 1023);
 +	        return y*2.0*const_0x1p1023;
--- a/modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch
+++ b/modules/fdlibm/patches/16_use_safer_strict_assign_on_visual_studio.patch
@@ -1,12 +1,12 @@
 diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
 --- a/modules/fdlibm/src/math_private.h
 +++ b/modules/fdlibm/src/math_private.h
-@@ -273,17 +273,17 @@ do {								\
+@@ -314,17 +314,17 @@ do {								\
  /* The above works on non-i386 too, but we use this to check v. */
  #define	LD80C(m, ex, v)	{ .e = (v), }
  #endif
  
  #ifdef FLT_EVAL_METHOD
  /*
   * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
   */
--- a/modules/fdlibm/patches/18_use_stdlib_sqrt.patch
+++ b/modules/fdlibm/patches/18_use_stdlib_sqrt.patch
@@ -182,25 +182,25 @@ diff --git a/modules/fdlibm/src/e_pow.cp
   * 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
  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,
- one	=  1.0,
-@@ -147,17 +149,17 @@ double
+@@ -151,17 +153,17 @@ double
  		    return (hy<0)?-y: zero;
  	    }
  	    if(iy==0x3ff00000) {	/* y is  +-1 */
  		if(hy<0) return one/x; else return x;
  	    }
  	    if(hy==0x40000000) return x*x; /* y is  2 */
  	    if(hy==0x3fe00000) {	/* y is  0.5 */
  		if(hx>=0)	/* x >= +0 */
new file mode 100644
--- /dev/null
+++ b/modules/fdlibm/patches/19_remove_unneeded_round_to_integer_helpers.patch
@@ -0,0 +1,130 @@
+diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
+--- a/modules/fdlibm/src/math_private.h
++++ b/modules/fdlibm/src/math_private.h
+@@ -586,126 +586,16 @@ CMPLXL(long double x, long double y)
+ 	REALPART(z) = x;
+ 	IMAGPART(z) = y;
+ 	return (z.f);
+ }
+ #endif
+ 
+ #endif /* _COMPLEX_H */
+  
+-/*
+- * The rnint() family rounds to the nearest integer for a restricted range
+- * range of args (up to about 2**MANT_DIG).  We assume that the current
+- * rounding mode is FE_TONEAREST so that this can be done efficiently.
+- * Extra precision causes more problems in practice, and we only centralize
+- * this here to reduce those problems, and have not solved the efficiency
+- * problems.  The exp2() family uses a more delicate version of this that
+- * requires extracting bits from the intermediate value, so it is not
+- * centralized here and should copy any solution of the efficiency problems.
+- */
+-
+-static inline double
+-rnint(__double_t x)
+-{
+-	/*
+-	 * This casts to double to kill any extra precision.  This depends
+-	 * on the cast being applied to a double_t to avoid compiler bugs
+-	 * (this is a cleaner version of STRICT_ASSIGN()).  This is
+-	 * inefficient if there actually is extra precision, but is hard
+-	 * to improve on.  We use double_t in the API to minimise conversions
+-	 * for just calling here.  Note that we cannot easily change the
+-	 * magic number to the one that works directly with double_t, since
+-	 * the rounding precision is variable at runtime on x86 so the
+-	 * magic number would need to be variable.  Assuming that the
+-	 * rounding precision is always the default is too fragile.  This
+-	 * and many other complications will move when the default is
+-	 * changed to FP_PE.
+-	 */
+-	return ((double)(x + 0x1.8p52) - 0x1.8p52);
+-}
+-
+-static inline float
+-rnintf(__float_t x)
+-{
+-	/*
+-	 * As for rnint(), except we could just call that to handle the
+-	 * extra precision case, usually without losing efficiency.
+-	 */
+-	return ((float)(x + 0x1.8p23F) - 0x1.8p23F);
+-}
+-
+-#ifdef LDBL_MANT_DIG
+-/*
+- * The complications for extra precision are smaller for rnintl() since it
+- * can safely assume that the rounding precision has been increased from
+- * its default to FP_PE on x86.  We don't exploit that here to get small
+- * optimizations from limiting the rangle to double.  We just need it for
+- * the magic number to work with long doubles.  ld128 callers should use
+- * rnint() instead of this if possible.  ld80 callers should prefer
+- * rnintl() since for amd64 this avoids swapping the register set, while
+- * for i386 it makes no difference (assuming FP_PE), and for other arches
+- * it makes little difference.
+- */
+-static inline long double
+-rnintl(long double x)
+-{
+-	return (x + __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2 -
+-	    __CONCAT(0x1.8p, LDBL_MANT_DIG) / 2);
+-}
+-#endif /* LDBL_MANT_DIG */
+-
+-/*
+- * irint() and i64rint() give the same result as casting to their integer
+- * return type provided their arg is a floating point integer.  They can
+- * sometimes be more efficient because no rounding is required.
+- */
+-#if (defined(amd64) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
+-#define	irint(x)						\
+-    (sizeof(x) == sizeof(float) &&				\
+-    sizeof(__float_t) == sizeof(long double) ? irintf(x) :	\
+-    sizeof(x) == sizeof(double) &&				\
+-    sizeof(__double_t) == sizeof(long double) ? irintd(x) :	\
+-    sizeof(x) == sizeof(long double) ? irintl(x) : (int)(x))
+-#else
+-#define	irint(x)	((int)(x))
+-#endif
+-
+-#define	i64rint(x)	((int64_t)(x))	/* only needed for ld128 so not opt. */
+-
+-#if defined(__i386__) && defined(__GNUCLIKE_ASM)
+-static __inline int
+-irintf(float x)
+-{
+-	int n;
+-
+-	__asm("fistl %0" : "=m" (n) : "t" (x));
+-	return (n);
+-}
+-
+-static __inline int
+-irintd(double x)
+-{
+-	int n;
+-
+-	__asm("fistl %0" : "=m" (n) : "t" (x));
+-	return (n);
+-}
+-#endif
+-
+-#if (defined(__amd64__) || defined(__i386__)) && defined(__GNUCLIKE_ASM)
+-static __inline int
+-irintl(long double x)
+-{
+-	int n;
+-
+-	__asm("fistl %0" : "=m" (n) : "t" (x));
+-	return (n);
+-}
+-#endif
+-
+ #ifdef DEBUG
+ #if defined(__amd64__) || defined(__i386__)
+ #define	breakpoint()	asm("int $3")
+ #else
+ #include <signal.h>
+ 
+ #define	breakpoint()	raise(SIGTRAP)
+ #endif
--- a/modules/fdlibm/src/e_atan2.cpp
+++ b/modules/fdlibm/src/e_atan2.cpp
@@ -64,17 +64,17 @@ double
 	u_int32_t lx,ly;
 
 	EXTRACT_WORDS(hx,lx,x);
 	ix = hx&0x7fffffff;
 	EXTRACT_WORDS(hy,ly,y);
 	iy = hy&0x7fffffff;
 	if(((ix|((lx|-lx)>>31))>0x7ff00000)||
 	   ((iy|((ly|-ly)>>31))>0x7ff00000))	/* x or y is NaN */
-	   return x+y;
+	    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 */
 	if((iy|ly)==0) {
 	    switch(m) {
 		case 0: 
 		case 1: return y; 	/* atan(+-0,+anything)=+-0 */
--- a/modules/fdlibm/src/e_exp.cpp
+++ b/modules/fdlibm/src/e_exp.cpp
@@ -142,19 +142,19 @@ double
 	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
 	    if(huge+x>one) return one+x;/* trigger inexact */
 	}
 	else k = 0;
 
     /* 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);
 	if(k >= -1021) {
 	    if (k==1024) {
 	        double const_0x1p1023 = pow(2, 1023);
 	        return y*2.0*const_0x1p1023;
 	    }
--- a/modules/fdlibm/src/e_hypot.cpp
+++ b/modules/fdlibm/src/e_hypot.cpp
@@ -65,17 +65,17 @@ double
 	a = fabs(a);
 	b = fabs(b);
 	if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
 	k=0;
 	if(ha > 0x5f300000) {	/* a>2**500 */
 	   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);
 	       if(((hb^0x7ff00000)|low)==0) w = b;
 	       return w;
 	   }
 	   /* scale a and b by 2**-600 */
 	   ha -= 0x25800000; hb -= 0x25800000;	k += 600;
--- a/modules/fdlibm/src/e_pow.cpp
+++ b/modules/fdlibm/src/e_pow.cpp
@@ -54,23 +54,27 @@
  * 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
 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 */
 huge	=  1.0e300,
 tiny    =  1.0e-300,
 	/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
 L1  =  5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
 L2  =  4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
@@ -111,31 +115,31 @@ double
 	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);
+	    return nan_mix(x, y);
 
     /* determine if y is an odd int when x < 0
      * yisint = 0	... y is not an integer
      * yisint = 1	... y is an odd int
      * yisint = 2	... y is an even int
      */
 	yisint  = 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 */
@@ -193,17 +197,17 @@ double
 		if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
 	    }
 	/* 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
 	   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;
 	    SET_LOW_WORD(t1,0);
 	    t2 = v-(t1-u);
 	} else {
 	    double ss,s2,s_h,s_l,t_h,t_l;
 	    n = 0;
@@ -230,30 +234,30 @@ double
 	    SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
 	    t_l = ax - (t_h-bp[k]);
 	    s_l = v*((u-s_h*t_h)-s_h*t_l);
 	/* compute log(ax) */
 	    s2 = ss*ss;
 	    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;
 	/* 2/(3log2)*(ss+...) */
 	    p_h = u+v;
 	    SET_LOW_WORD(p_h,0);
 	    p_l = v-(p_h-u);
 	    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);
 	}
 
     /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
 	y1  = y;
 	SET_LOW_WORD(y1,0);
--- a/modules/fdlibm/src/math_private.h
+++ b/modules/fdlibm/src/math_private.h
@@ -40,16 +40,57 @@
 
 #ifndef u_int32_t
 #define u_int32_t uint32_t
 #endif
 #ifndef u_int64_t
 #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.
  */
 
 #if MOZ_BIG_ENDIAN
 
 typedef union
@@ -302,18 +343,19 @@ do {								\
 		__lval = (rval);		\
 		(lval) = __lval;		\
 	}					\
 } while (0)
 #endif /* FLT_EVAL_METHOD */
 
 /* 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)	\
 		fpsetprec(FP_PE)
 #define	RETURNI(x) do {				\
 	__retval = (x);				\
 	if (__oprec != FP_PE)			\
 		fpsetprec(__oprec);		\
@@ -326,16 +368,17 @@ do {								\
 		fpsetprec(FP_PE)
 #define	RETURNV() do {				\
 	if (__oprec != FP_PE)			\
 		fpsetprec(__oprec);		\
 	return;			\
 } while (0)
 #else
 #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. */
 #define      RETURNF(v)      return (v)
 
@@ -443,16 +486,41 @@ do {								\
 	(a) = __tmp;		\
 } while (0)
 
 /*
  * Common routine to process the arguments to nan(), nanf(), and nanl().
  */
 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
 
 /*
  * C99 specifies that complex numbers have the same representation as
  * an array of two elements, where the first element is the real part
  * and the second element is the imaginary part.
  */
 typedef union {
@@ -518,58 +586,16 @@ CMPLXL(long double x, long double y)
 	REALPART(z) = x;
 	IMAGPART(z) = y;
 	return (z.f);
 }
 #endif
 
 #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")
 #else
 #include <signal.h>
 
 #define	breakpoint()	raise(SIGTRAP)
 #endif
--- a/modules/fdlibm/src/s_cbrt.cpp
+++ b/modules/fdlibm/src/s_cbrt.cpp
@@ -10,16 +10,17 @@
  * ====================================================
  *
  * Optimized by Bruce D. Evans.
  */
 
 //#include <sys/cdefs.h>
 //__FBSDID("$FreeBSD$");
 
+#include <float.h>
 #include "math_private.h"
 
 /* cbrt(x)
  * Return cube root of x
  */
 static const u_int32_t
 	B1 = 715094163, /* B1 = (1023-1023/3-0.03306235651)*2**20 */
 	B2 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
--- a/modules/fdlibm/src/s_expm1.cpp
+++ b/modules/fdlibm/src/s_expm1.cpp
@@ -182,17 +182,17 @@ expm1(double x)
     /* x is now in primary range */
 	hfx = 0.5*x;
 	hxs = x*hfx;
 	r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
 	t  = 3.0-r1*hfx;
 	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;
 	    if(k==1) {
 	       	if(x < -0.25) return -2.0*(e-(x+0.5));
 	       	else 	      return  one+2.0*(x-e);
 	    }
 	    if (k <= -2 || k>56) {   /* suffice to return exp(x)-1 */