Add m-esr52 at 52.6.0

37 files changed, 4590 insertions, 0 deletions

diff --git a/modules/fdlibm/src/e_acos.cpp b/modules/fdlibm/src/e_acos.cpp
new file mode 100644
index 000000000..12be296cb
--- /dev/null
+++ b/modules/fdlibm/src/e_acos.cpp
@@ -0,0 +1,106 @@
+
+/* @(#)e_acos.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_acos(x)
+ * Method :                  
+ *	acos(x)  = pi/2 - asin(x)
+ *	acos(-x) = pi/2 + asin(x)
+ * For |x|<=0.5
+ *	acos(x) = pi/2 - (x + x*x^2*R(x^2))	(see asin.c)
+ * For x>0.5
+ * 	acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
+ *		= 2asin(sqrt((1-x)/2))  
+ *		= 2s + 2s*z*R(z) 	...z=(1-x)/2, s=sqrt(z)
+ *		= 2f + (2c + 2s*z*R(z))
+ *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
+ *     for f so that f+c ~ sqrt(z).
+ * For x<-0.5
+ *	acos(x) = pi - 2asin(sqrt((1-|x|)/2))
+ *		= pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ * Function needed: sqrt
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+one=  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+pi =  3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+pio2_hi =  1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+static volatile double
+pio2_lo =  6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
+static const double
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double
+__ieee754_acos(double x)
+{
+	double z,p,q,r,w,s,c,df;
+	int32_t hx,ix;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x3ff00000) {	/* |x| >= 1 */
+	    u_int32_t lx;
+	    GET_LOW_WORD(lx,x);
+	    if(((ix-0x3ff00000)|lx)==0) {	/* |x|==1 */
+		if(hx>0) return 0.0;		/* acos(1) = 0  */
+		else return pi+2.0*pio2_lo;	/* acos(-1)= pi */
+	    }
+	    return (x-x)/(x-x);		/* acos(|x|>1) is NaN */
+	}
+	if(ix<0x3fe00000) {	/* |x| < 0.5 */
+	    if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
+	    z = x*x;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    return pio2_hi - (x - (pio2_lo-x*r));
+	} else  if (hx<0) {		/* x < -0.5 */
+	    z = (one+x)*0.5;
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    s = sqrt(z);
+	    r = p/q;
+	    w = r*s-pio2_lo;
+	    return pi - 2.0*(s+w);
+	} else {			/* x > 0.5 */
+	    z = (one-x)*0.5;
+	    s = sqrt(z);
+	    df = s;
+	    SET_LOW_WORD(df,0);
+	    c  = (z-df*df)/(s+df);
+	    p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+	    q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+	    r = p/q;
+	    w = r*s+c;
+	    return 2.0*(df+w);
+	}
+}
diff --git a/modules/fdlibm/src/e_acosh.cpp b/modules/fdlibm/src/e_acosh.cpp
new file mode 100644
index 000000000..bdabcec3e
--- /dev/null
+++ b/modules/fdlibm/src/e_acosh.cpp
@@ -0,0 +1,63 @@
+
+/* @(#)e_acosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_acosh(x)
+ * Method :
+ *	Based on 
+ *		acosh(x) = log [ x + sqrt(x*x-1) ]
+ *	we have
+ *		acosh(x) := log(x)+ln2,	if x is large; else
+ *		acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+ *		acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ *	acosh(x) is NaN with signal if x<1.
+ *	acosh(NaN) is NaN without signal.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+one	= 1.0,
+ln2	= 6.93147180559945286227e-01;  /* 0x3FE62E42, 0xFEFA39EF */
+
+double
+__ieee754_acosh(double x)
+{
+	double t;
+	int32_t hx;
+	u_int32_t lx;
+	EXTRACT_WORDS(hx,lx,x);
+	if(hx<0x3ff00000) {		/* x < 1 */
+	    return (x-x)/(x-x);
+	} else if(hx >=0x41b00000) {	/* x > 2**28 */
+	    if(hx >=0x7ff00000) {	/* x is inf of NaN */
+	        return x+x;
+	    } else 
+		return __ieee754_log(x)+ln2;	/* acosh(huge)=log(2x) */
+	} else if(((hx-0x3ff00000)|lx)==0) {
+	    return 0.0;			/* acosh(1) = 0 */
+	} else if (hx > 0x40000000) {	/* 2**28 > x > 2 */
+	    t=x*x;
+	    return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
+	} else {			/* 1<x<2 */
+	    t = x-one;
+	    return log1p(t+sqrt(2.0*t+t*t));
+	}
+}
diff --git a/modules/fdlibm/src/e_asin.cpp b/modules/fdlibm/src/e_asin.cpp
new file mode 100644
index 000000000..396f49449
--- /dev/null
+++ b/modules/fdlibm/src/e_asin.cpp
@@ -0,0 +1,112 @@
+
+/* @(#)e_asin.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_asin(x)
+ * Method :                  
+ *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *	we approximate asin(x) on [0,0.5] by
+ *		asin(x) = x + x*x^2*R(x^2)
+ *	where
+ *		R(x^2) is a rational approximation of (asin(x)-x)/x^3 
+ *	and its remez error is bounded by
+ *		|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
+ *
+ *	For x in [0.5,1]
+ *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
+ *	then for x>0.98
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
+ *	For x<=0.98, let pio4_hi = pio2_hi/2, then
+ *		f = hi part of s;
+ *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
+ *	and
+ *		asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
+ *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
+ *
+ * Special cases:
+ *	if x is NaN, return x itself;
+ *	if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+huge =  1.000e+300,
+pio2_hi =  1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
+pio2_lo =  6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
+pio4_hi =  7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
+	/* coefficient for R(x^2) */
+pS0 =  1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 =  2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 =  7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 =  3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 =  2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double
+__ieee754_asin(double x)
+{
+	double t=0.0,w,p,q,c,r,s;
+	int32_t hx,ix;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>= 0x3ff00000) {		/* |x|>= 1 */
+	    u_int32_t lx;
+	    GET_LOW_WORD(lx,x);
+	    if(((ix-0x3ff00000)|lx)==0)
+		    /* asin(1)=+-pi/2 with inexact */
+		return x*pio2_hi+x*pio2_lo;	
+	    return (x-x)/(x-x);		/* asin(|x|>1) is NaN */   
+	} else if (ix<0x3fe00000) {	/* |x|<0.5 */
+	    if(ix<0x3e500000) {		/* if |x| < 2**-26 */
+		if(huge+x>one) return x;/* return x with inexact if x!=0*/
+	    }
+	    t = x*x;
+	    p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+	    q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+	    w = p/q;
+	    return x+x*w;
+	}
+	/* 1> |x|>= 0.5 */
+	w = one-fabs(x);
+	t = w*0.5;
+	p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+	q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+	s = sqrt(t);
+	if(ix>=0x3FEF3333) { 	/* if |x| > 0.975 */
+	    w = p/q;
+	    t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+	} else {
+	    w  = s;
+	    SET_LOW_WORD(w,0);
+	    c  = (t-w*w)/(s+w);
+	    r  = p/q;
+	    p  = 2.0*s*r-(pio2_lo-2.0*c);
+	    q  = pio4_hi-2.0*w;
+	    t  = pio4_hi-(p-q);
+	}    
+	if(hx>0) return t; else return -t;    
+}
diff --git a/modules/fdlibm/src/e_atan2.cpp b/modules/fdlibm/src/e_atan2.cpp
new file mode 100644
index 000000000..9990072cf
--- /dev/null
+++ b/modules/fdlibm/src/e_atan2.cpp
@@ -0,0 +1,124 @@
+
+/* @(#)e_atan2.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_atan2(y,x)
+ * Method :
+ *	1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+ *	2. Reduce x to positive by (if x and y are unexceptional): 
+ *		ARG (x+iy) = arctan(y/x)   	   ... if x > 0,
+ *		ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
+ *
+ * Special cases:
+ *
+ *	ATAN2((anything), NaN ) is NaN;
+ *	ATAN2(NAN , (anything) ) is NaN;
+ *	ATAN2(+-0, +(anything but NaN)) is +-0  ;
+ *	ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ *	ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ *	ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ *	ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ *	ATAN2(+-INF,+INF ) is +-pi/4 ;
+ *	ATAN2(+-INF,-INF ) is +-3pi/4;
+ *	ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * 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_private.h"
+
+static volatile double
+tiny  = 1.0e-300;
+static const double
+zero  = 0.0,
+pi_o_4  = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
+pi_o_2  = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
+pi      = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
+static volatile double
+pi_lo   = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
+
+double
+__ieee754_atan2(double y, double x)
+{
+	double z;
+	int32_t k,m,hx,hy,ix,iy;
+	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;
+	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 */
+		case 2: return  pi+tiny;/* atan(+0,-anything) = pi */
+		case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
+	    }
+	}
+    /* when x = 0 */
+	if((ix|lx)==0) return (hy<0)?  -pi_o_2-tiny: pi_o_2+tiny;
+	    
+    /* when x is INF */
+	if(ix==0x7ff00000) {
+	    if(iy==0x7ff00000) {
+		switch(m) {
+		    case 0: return  pi_o_4+tiny;/* atan(+INF,+INF) */
+		    case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
+		    case 2: return  3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
+		    case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
+		}
+	    } else {
+		switch(m) {
+		    case 0: return  zero  ;	/* atan(+...,+INF) */
+		    case 1: return -zero  ;	/* atan(-...,+INF) */
+		    case 2: return  pi+tiny  ;	/* atan(+...,-INF) */
+		    case 3: return -pi-tiny  ;	/* atan(-...,-INF) */
+		}
+	    }
+	}
+    /* when y is INF */
+	if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
+
+    /* compute y/x */
+	k = (iy-ix)>>20;
+	if(k > 60) {		 	/* |y/x| >  2**60 */
+	    z=pi_o_2+0.5*pi_lo;
+	    m&=1;
+	}
+	else if(hx<0&&k<-60) z=0.0; 	/* 0 > |y|/x > -2**-60 */
+	else z=atan(fabs(y/x));		/* safe to do y/x */
+	switch (m) {
+	    case 0: return       z  ;	/* atan(+,+) */
+	    case 1: return      -z  ;	/* atan(-,+) */
+	    case 2: return  pi-(z-pi_lo);/* atan(+,-) */
+	    default: /* case 3 */
+	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
+	}
+}
diff --git a/modules/fdlibm/src/e_atanh.cpp b/modules/fdlibm/src/e_atanh.cpp
new file mode 100644
index 000000000..a8f0f8deb
--- /dev/null
+++ b/modules/fdlibm/src/e_atanh.cpp
@@ -0,0 +1,63 @@
+
+/* @(#)e_atanh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ *
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_atanh(x)
+ * Method :
+ *    1.Reduced x to positive by atanh(-x) = -atanh(x)
+ *    2.For x>=0.5
+ *                  1              2x                          x
+ *	atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ *                  2             1 - x                      1 - x
+ *	
+ * 	For x<0.5
+ *	atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ *	atanh(x) is NaN if |x| > 1 with signal;
+ *	atanh(NaN) is that NaN with no signal;
+ *	atanh(+-1) is +-INF with signal.
+ *
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double one = 1.0, huge = 1e300;
+static const double zero = 0.0;
+
+double
+__ieee754_atanh(double x)
+{
+	double t;
+	int32_t hx,ix;
+	u_int32_t lx;
+	EXTRACT_WORDS(hx,lx,x);
+	ix = hx&0x7fffffff;
+	if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
+	    return (x-x)/(x-x);
+	if(ix==0x3ff00000) 
+	    return x/zero;
+	if(ix<0x3e300000&&(huge+x)>zero) return x;	/* x<2**-28 */
+	SET_HIGH_WORD(x,ix);
+	if(ix<0x3fe00000) {		/* x < 0.5 */
+	    t = x+x;
+	    t = 0.5*log1p(t+t*x/(one-x));
+	} else 
+	    t = 0.5*log1p((x+x)/(one-x));
+	if(hx>=0) return t; else return -t;
+}
diff --git a/modules/fdlibm/src/e_cosh.cpp b/modules/fdlibm/src/e_cosh.cpp
new file mode 100644
index 000000000..42cb277d4
--- /dev/null
+++ b/modules/fdlibm/src/e_cosh.cpp
@@ -0,0 +1,80 @@
+
+/* @(#)e_cosh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_cosh(x)
+ * Method : 
+ * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
+ *	1. Replace x by |x| (cosh(x) = cosh(-x)). 
+ *	2. 
+ *		                                        [ exp(x) - 1 ]^2 
+ *	    0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
+ *			       			           2*exp(x)
+ *
+ *		                                  exp(x) +  1/exp(x)
+ *	    ln2/2    <= x <= 22     :  cosh(x) := -------------------
+ *			       			          2
+ *	    22       <= x <= lnovft :  cosh(x) := exp(x)/2 
+ *	    lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
+ *	    ln2ovft  <  x	    :  cosh(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ *	cosh(x) is |x| if x is +INF, -INF, or NaN.
+ *	only cosh(0)=1 is exact for finite x.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double one = 1.0, half=0.5, huge = 1.0e300;
+
+double
+__ieee754_cosh(double x)
+{
+	double t,w;
+	int32_t ix;
+
+    /* High word of |x|. */
+	GET_HIGH_WORD(ix,x);
+	ix &= 0x7fffffff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x*x;	
+
+    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+	if(ix<0x3fd62e43) {
+	    t = expm1(fabs(x));
+	    w = one+t;
+	    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
+	    return one+(t*t)/(w+w);
+	}
+
+    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+	if (ix < 0x40360000) {
+		t = __ieee754_exp(fabs(x));
+		return half*t+half/t;
+	}
+
+    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+	if (ix < 0x40862E42)  return half*__ieee754_exp(fabs(x));
+
+    /* |x| in [log(maxdouble), overflowthresold] */
+	if (ix<=0x408633CE)
+	    return __ldexp_exp(fabs(x), -1);
+
+    /* |x| > overflowthresold, cosh(x) overflow */
+	return huge*huge;
+}
diff --git a/modules/fdlibm/src/e_exp.cpp b/modules/fdlibm/src/e_exp.cpp
new file mode 100644
index 000000000..b31979134
--- /dev/null
+++ b/modules/fdlibm/src/e_exp.cpp
@@ -0,0 +1,162 @@
+
+/* @(#)e_exp.c 1.6 04/04/22 */
+/*
+ * ====================================================
+ * 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 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ *   1. Argument reduction:
+ *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ *	Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.  
+ *
+ *      Here r will be represented as r = hi-lo for better 
+ *	accuracy.
+ *
+ *   2. Approximation of exp(r) by a special rational function on
+ *	the interval [0,0.34658]:
+ *	Write
+ *	    R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ *      We use a special Remes algorithm on [0,0.34658] to generate 
+ * 	a polynomial of degree 5 to approximate R. The maximum error 
+ *	of this polynomial approximation is bounded by 2**-59. In
+ *	other words,
+ *	    R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ *  	(where z=r*r, and the values of P1 to P5 are listed below)
+ *	and
+ *	    |                  5          |     -59
+ *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2 
+ *	    |                             |
+ *	The computation of exp(r) thus becomes
+ *                             2*r
+ *		exp(r) = 1 + -------
+ *		              R - r
+ *                                 r*R1(r)	
+ *		       = 1 + r + ----------- (for better accuracy)
+ *		                  2 - R1(r)
+ *	where
+ *			         2       4             10
+ *		R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+ *	
+ *   3. Scale back to obtain exp(x):
+ *	From step 1, we have
+ *	   exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ *	exp(INF) is INF, exp(NaN) is NaN;
+ *	exp(-INF) is 0, and
+ *	for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *	For IEEE double 
+ *	    if x >  7.09782712893383973096e+02 then exp(x) overflow
+ *	    if x < -7.45133219101941108420e+02 then exp(x) underflow
+ *
+ * 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_private.h"
+
+static const double
+one	= 1.0,
+halF[2]	= {0.5,-0.5,},
+o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
+u_threshold= -7.45133219101941108420e+02,  /* 0xc0874910, 0xD52D3051 */
+ln2HI[2]   ={ 6.93147180369123816490e-01,  /* 0x3fe62e42, 0xfee00000 */
+	     -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
+ln2LO[2]   ={ 1.90821492927058770002e-10,  /* 0x3dea39ef, 0x35793c76 */
+	     -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
+invln2 =  1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
+P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+static volatile double
+huge	= 1.0e+300,
+twom1000= 9.33263618503218878990e-302;     /* 2**-1000=0x01700000,0*/
+
+double
+__ieee754_exp(double x)	/* default IEEE double exp */
+{
+	double y,hi=0.0,lo=0.0,c,t,twopk;
+	int32_t k=0,xsb;
+	u_int32_t hx;
+
+	GET_HIGH_WORD(hx,x);
+	xsb = (hx>>31)&1;		/* sign bit of x */
+	hx &= 0x7fffffff;		/* high word of |x| */
+
+    /* filter out non-finite argument */
+	if(hx >= 0x40862E42) {			/* if |x|>=709.78... */
+            if(hx>=0x7ff00000) {
+	        u_int32_t lx;
+		GET_LOW_WORD(lx,x);
+		if(((hx&0xfffff)|lx)!=0)
+		     return x+x; 		/* NaN */
+		else return (xsb==0)? x:0.0;	/* exp(+-inf)={inf,0} */
+	    }
+	    if(x > o_threshold) return huge*huge; /* overflow */
+	    if(x < u_threshold) return twom1000*twom1000; /* underflow */
+	}
+
+    /* argument reduction */
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */ 
+	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
+		hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
+	    } else {
+		k  = (int)(invln2*x+halF[xsb]);
+		t  = k;
+		hi = x - t*ln2HI[0];	/* t*ln2HI is exact here */
+		lo = t*ln2LO[0];
+	    }
+	    STRICT_ASSIGN(double, x, hi - lo);
+	} 
+	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);
+	else
+	    INSERT_WORDS(twopk,0x3ff00000+((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;
+	    }
+	    return y*twopk;
+	} else {
+	    return y*twopk*twom1000;
+	}
+}
diff --git a/modules/fdlibm/src/e_hypot.cpp b/modules/fdlibm/src/e_hypot.cpp
new file mode 100644
index 000000000..f5c7037bb
--- /dev/null
+++ b/modules/fdlibm/src/e_hypot.cpp
@@ -0,0 +1,126 @@
+
+/* @(#)e_hypot.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_hypot(x,y)
+ *
+ * Method :                  
+ *	If (assume round-to-nearest) z=x*x+y*y 
+ *	has error less than sqrt(2)/2 ulp, than 
+ *	sqrt(z) has error less than 1 ulp (exercise).
+ *
+ *	So, compute sqrt(x*x+y*y) with some care as 
+ *	follows to get the error below 1 ulp:
+ *
+ *	Assume x>y>0;
+ *	(if possible, set rounding to round-to-nearest)
+ *	1. if x > 2y  use
+ *		x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
+ *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
+ *	2. if x <= 2y use
+ *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
+ *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, 
+ *	y1= y with lower 32 bits chopped, y2 = y-y1.
+ *		
+ *	NOTE: scaling may be necessary if some argument is too 
+ *	      large or too tiny
+ *
+ * Special cases:
+ *	hypot(x,y) is INF if x or y is +INF or -INF; else
+ *	hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * 	hypot(x,y) returns sqrt(x^2+y^2) with error less 
+ * 	than 1 ulps (units in the last place) 
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+double
+__ieee754_hypot(double x, double y)
+{
+	double a,b,t1,t2,y1,y2,w;
+	int32_t j,k,ha,hb;
+
+	GET_HIGH_WORD(ha,x);
+	ha &= 0x7fffffff;
+	GET_HIGH_WORD(hb,y);
+	hb &= 0x7fffffff;
+	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
+	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);
+	       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;
+	   SET_HIGH_WORD(a,ha);
+	   SET_HIGH_WORD(b,hb);
+	}
+	if(hb < 0x20b00000) {	/* b < 2**-500 */
+	    if(hb <= 0x000fffff) {	/* subnormal b or 0 */
+	        u_int32_t low;
+		GET_LOW_WORD(low,b);
+		if((hb|low)==0) return a;
+		t1=0;
+		SET_HIGH_WORD(t1,0x7fd00000);	/* t1=2^1022 */
+		b *= t1;
+		a *= t1;
+		k -= 1022;
+	    } else {		/* scale a and b by 2^600 */
+	        ha += 0x25800000; 	/* a *= 2^600 */
+		hb += 0x25800000;	/* b *= 2^600 */
+		k -= 600;
+		SET_HIGH_WORD(a,ha);
+		SET_HIGH_WORD(b,hb);
+	    }
+	}
+    /* medium size a and b */
+	w = a-b;
+	if (w>b) {
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha);
+	    t2 = a-t1;
+	    w  = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
+	} else {
+	    a  = a+a;
+	    y1 = 0;
+	    SET_HIGH_WORD(y1,hb);
+	    y2 = b - y1;
+	    t1 = 0;
+	    SET_HIGH_WORD(t1,ha+0x00100000);
+	    t2 = a - t1;
+	    w  = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+	}
+	if(k!=0) {
+	    u_int32_t high;
+	    t1 = 1.0;
+	    GET_HIGH_WORD(high,t1);
+	    SET_HIGH_WORD(t1,high+(k<<20));
+	    return t1*w;
+	} else return w;
+}
diff --git a/modules/fdlibm/src/e_log.cpp b/modules/fdlibm/src/e_log.cpp
new file mode 100644
index 000000000..fa2da8fcb
--- /dev/null
+++ b/modules/fdlibm/src/e_log.cpp
@@ -0,0 +1,142 @@
+
+/* @(#)e_log.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_log(x)
+ * Return the logrithm of x
+ *
+ * Method :                  
+ *   1. Argument Reduction: find k and f such that 
+ *			x = 2^k * (1+f), 
+ *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *	     	 = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate 
+ * 	a polynomial of degree 14 to approximate R The maximum error 
+ *	of this polynomial approximation is bounded by 2**-58.45. In
+ *	other words,
+ *		        2      4      6      8      10      12      14
+ *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *  	(the values of Lg1 to Lg7 are listed in the program)
+ *	and
+ *	    |      2          14          |     -58.45
+ *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
+ *	    |                             |
+ *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *	In order to guarantee error in log below 1ulp, we compute log
+ *	by
+ *		log(1+f) = f - s*(f - R)	(if f is not too large)
+ *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
+ *	
+ *	3. Finally,  log(x) = k*ln2 + log(1+f).  
+ *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *	   Here ln2 is split into two floating point number: 
+ *			ln2_hi + ln2_lo,
+ *	   where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *	log(x) is NaN with signal if x < 0 (including -INF) ; 
+ *	log(+INF) is +INF; log(0) is -INF with signal;
+ *	log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * 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_private.h"
+
+static const double
+ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
+ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
+two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
+Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+
+static const double zero   =  0.0;
+static volatile double vzero = 0.0;
+
+double
+__ieee754_log(double x)
+{
+	double hfsq,f,s,z,R,w,t1,t2,dk;
+	int32_t k,hx,i,j;
+	u_int32_t lx;
+
+	EXTRACT_WORDS(hx,lx,x);
+
+	k=0;
+	if (hx < 0x00100000) {			/* x < 2**-1022  */
+	    if (((hx&0x7fffffff)|lx)==0) 
+		return -two54/vzero;		/* log(+-0)=-inf */
+	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
+	    k -= 54; x *= two54; /* subnormal number, scale up x */
+	    GET_HIGH_WORD(hx,x);
+	} 
+	if (hx >= 0x7ff00000) return x+x;
+	k += (hx>>20)-1023;
+	hx &= 0x000fffff;
+	i = (hx+0x95f64)&0x100000;
+	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
+	k += (i>>20);
+	f = x-1.0;
+	if((0x000fffff&(2+hx))<3) {	/* -2**-20 <= f < 2**-20 */
+	    if(f==zero) {
+		if(k==0) {
+		    return zero;
+		} else {
+		    dk=(double)k;
+		    return dk*ln2_hi+dk*ln2_lo;
+		}
+	    }
+	    R = f*f*(0.5-0.33333333333333333*f);
+	    if(k==0) return f-R; else {dk=(double)k;
+	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
+	}
+ 	s = f/(2.0+f); 
+	dk = (double)k;
+	z = s*s;
+	i = hx-0x6147a;
+	w = z*z;
+	j = 0x6b851-hx;
+	t1= w*(Lg2+w*(Lg4+w*Lg6)); 
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
+	i |= j;
+	R = t2+t1;
+	if(i>0) {
+	    hfsq=0.5*f*f;
+	    if(k==0) return f-(hfsq-s*(hfsq+R)); else
+		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
+	} else {
+	    if(k==0) return f-s*(f-R); else
+		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
+	}
+}
diff --git a/modules/fdlibm/src/e_log10.cpp b/modules/fdlibm/src/e_log10.cpp
new file mode 100644
index 000000000..ed6879885
--- /dev/null
+++ b/modules/fdlibm/src/e_log10.cpp
@@ -0,0 +1,89 @@
+
+/* @(#)e_log10.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * Return the base 10 logarithm of x.  See e_log.c and k_log.h for most
+ * comments.
+ *
+ *    log10(x) = (f - 0.5*f*f + k_log1p(f)) / ln10 + k * log10(2)
+ * in not-quite-routine extra precision.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+#include "k_log.h"
+ 
+static const double
+two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+ivln10hi   =  4.34294481878168880939e-01, /* 0x3fdbcb7b, 0x15200000 */
+ivln10lo   =  2.50829467116452752298e-11, /* 0x3dbb9438, 0xca9aadd5 */
+log10_2hi  =  3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
+log10_2lo  =  3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
+
+static const double zero   =  0.0;
+static volatile double vzero = 0.0;
+
+double
+__ieee754_log10(double x)
+{
+	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y,y2;
+	int32_t i,k,hx;
+	u_int32_t lx;
+
+	EXTRACT_WORDS(hx,lx,x);
+
+	k=0;
+	if (hx < 0x00100000) {			/* x < 2**-1022  */
+	    if (((hx&0x7fffffff)|lx)==0)
+		return -two54/vzero;		/* log(+-0)=-inf */
+	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
+	    k -= 54; x *= two54; /* subnormal number, scale up x */
+	    GET_HIGH_WORD(hx,x);
+	}
+	if (hx >= 0x7ff00000) return x+x;
+	if (hx == 0x3ff00000 && lx == 0)
+	    return zero;			/* log(1) = +0 */
+	k += (hx>>20)-1023;
+	hx &= 0x000fffff;
+	i = (hx+0x95f64)&0x100000;
+	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
+	k += (i>>20);
+	y = (double)k;
+	f = x - 1.0;
+	hfsq = 0.5*f*f;
+	r = k_log1p(f);
+
+	/* See e_log2.c for most details. */
+	hi = f - hfsq;
+	SET_LOW_WORD(hi,0);
+	lo = (f - hi) - hfsq + r;
+	val_hi = hi*ivln10hi;
+	y2 = y*log10_2hi;
+	val_lo = y*log10_2lo + (lo+hi)*ivln10lo + lo*ivln10hi;
+
+	/*
+	 * Extra precision in for adding y*log10_2hi is not strictly needed
+	 * since there is no very large cancellation near x = sqrt(2) or
+	 * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
+	 * with some parallelism and it reduces the error for many args.
+	 */
+	w = y2 + val_hi;
+	val_lo += (y2 - w) + val_hi;
+	val_hi = w;
+
+	return val_lo + val_hi;
+}
diff --git a/modules/fdlibm/src/e_log2.cpp b/modules/fdlibm/src/e_log2.cpp
new file mode 100644
index 000000000..5649fec44
--- /dev/null
+++ b/modules/fdlibm/src/e_log2.cpp
@@ -0,0 +1,112 @@
+
+/* @(#)e_log10.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * Return the base 2 logarithm of x.  See e_log.c and k_log.h for most
+ * comments.
+ *
+ * This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
+ * then does the combining and scaling steps
+ *    log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
+ * in not-quite-routine extra precision.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+#include "k_log.h"
+
+static const double
+two54      =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+ivln2hi    =  1.44269504072144627571e+00, /* 0x3ff71547, 0x65200000 */
+ivln2lo    =  1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
+
+static const double zero   =  0.0;
+static volatile double vzero = 0.0;
+
+double
+__ieee754_log2(double x)
+{
+	double f,hfsq,hi,lo,r,val_hi,val_lo,w,y;
+	int32_t i,k,hx;
+	u_int32_t lx;
+
+	EXTRACT_WORDS(hx,lx,x);
+
+	k=0;
+	if (hx < 0x00100000) {			/* x < 2**-1022  */
+	    if (((hx&0x7fffffff)|lx)==0)
+		return -two54/vzero;		/* log(+-0)=-inf */
+	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
+	    k -= 54; x *= two54; /* subnormal number, scale up x */
+	    GET_HIGH_WORD(hx,x);
+	}
+	if (hx >= 0x7ff00000) return x+x;
+	if (hx == 0x3ff00000 && lx == 0)
+	    return zero;			/* log(1) = +0 */
+	k += (hx>>20)-1023;
+	hx &= 0x000fffff;
+	i = (hx+0x95f64)&0x100000;
+	SET_HIGH_WORD(x,hx|(i^0x3ff00000));	/* normalize x or x/2 */
+	k += (i>>20);
+	y = (double)k;
+	f = x - 1.0;
+	hfsq = 0.5*f*f;
+	r = k_log1p(f);
+
+	/*
+	 * f-hfsq must (for args near 1) be evaluated in extra precision
+	 * to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
+	 * This is fairly efficient since f-hfsq only depends on f, so can
+	 * be evaluated in parallel with R.  Not combining hfsq with R also
+	 * keeps R small (though not as small as a true `lo' term would be),
+	 * so that extra precision is not needed for terms involving R.
+	 *
+	 * Compiler bugs involving extra precision used to break Dekker's
+	 * theorem for spitting f-hfsq as hi+lo, unless double_t was used
+	 * or the multi-precision calculations were avoided when double_t
+	 * has extra precision.  These problems are now automatically
+	 * avoided as a side effect of the optimization of combining the
+	 * Dekker splitting step with the clear-low-bits step.
+	 *
+	 * y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
+	 * precision to avoid a very large cancellation when x is very near
+	 * these values.  Unlike the above cancellations, this problem is
+	 * specific to base 2.  It is strange that adding +-1 is so much
+	 * harder than adding +-ln2 or +-log10_2.
+	 *
+	 * This uses Dekker's theorem to normalize y+val_hi, so the
+	 * compiler bugs are back in some configurations, sigh.  And I
+	 * don't want to used double_t to avoid them, since that gives a
+	 * pessimization and the support for avoiding the pessimization
+	 * is not yet available.
+	 *
+	 * The multi-precision calculations for the multiplications are
+	 * routine.
+	 */
+	hi = f - hfsq;
+	SET_LOW_WORD(hi,0);
+	lo = (f - hi) - hfsq + r;
+	val_hi = hi*ivln2hi;
+	val_lo = (lo+hi)*ivln2lo + lo*ivln2hi;
+
+	/* spadd(val_hi, val_lo, y), except for not using double_t: */
+	w = y + val_hi;
+	val_lo += (y - w) + val_hi;
+	val_hi = w;
+
+	return val_lo + val_hi;
+}
diff --git a/modules/fdlibm/src/e_pow.cpp b/modules/fdlibm/src/e_pow.cpp
new file mode 100644
index 000000000..366e3933b
--- /dev/null
+++ b/modules/fdlibm/src/e_pow.cpp
@@ -0,0 +1,305 @@
+/* @(#)e_pow.c 1.5 04/04/22 SMI */
+/*
+ * ====================================================
+ * 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 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_pow(x,y) return x**y
+ *
+ *		      n
+ * Method:  Let x =  2   * (1+f)
+ *	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 
+ *	   arithmetic, where |y'|<=0.5.
+ *	3. Return x**y = 2**n*exp(y'*log2)
+ *
+ * Special cases:
+ *	1.  (anything) ** 0  is 1
+ *	2.  (anything) ** 1  is itself
+ *	3.  (anything) ** NAN is NAN except 1 ** NAN = 1
+ *	4.  NAN ** (anything except 0) is NAN
+ *	5.  +-(|x| > 1) **  +INF is +INF
+ *	6.  +-(|x| > 1) **  -INF is +0
+ *	7.  +-(|x| < 1) **  +INF is +0
+ *	8.  +-(|x| < 1) **  -INF is +INF
+ *	9.  +-1         ** +-INF is 1
+ *	10. +0 ** (+anything except 0, NAN)               is +0
+ *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
+ *	12. +0 ** (-anything except 0, NAN)               is +INF
+ *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
+ *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
+ *	15. +INF ** (+anything except 0,NAN) is +INF
+ *	16. +INF ** (-anything except 0,NAN) is +0
+ *	17. -INF ** (anything)  = -0 ** (-anything)
+ *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+ *	19. (-anything except 0 and inf) ** (non-integer) is NAN
+ *
+ * Accuracy:
+ *	pow(x,y) returns x**y nearly rounded. In particular
+ *			pow(integer,integer)
+ *	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 
+ * to produce the hexadecimal values shown.
+ */
+
+#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,
+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 */
+L3  =  3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
+L4  =  2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
+L5  =  2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
+L6  =  2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
+P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5   =  4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
+lg2  =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
+lg2_h  =  6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
+lg2_l  = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
+ovt =  8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
+cp    =  9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
+cp_h  =  9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
+cp_l  = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
+ivln2    =  1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
+ivln2_h  =  1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
+ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
+
+double
+__ieee754_pow(double x, double y)
+{
+	double z,ax,z_h,z_l,p_h,p_l;
+	double y1,t1,t2,r,s,t,u,v,w;
+	int32_t i,j,k,yisint,n;
+	int32_t hx,hy,ix,iy;
+	u_int32_t lx,ly;
+
+	EXTRACT_WORDS(hx,lx,x);
+	EXTRACT_WORDS(hy,ly,y);
+	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
+
+    /* y==zero: x**0 = 1 */
+	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);
+
+    /* 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);
+		} 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 (iy==0x7ff00000) {	/* y is +-inf */
+	        if(((ix-0x3ff00000)|lx)==0)
+		    return  one;	/* (-1)**+-inf is 1 */
+	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
+		    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;
+	    }
+	    if(hy==0x40000000) return x*x; /* y is  2 */
+	    if(hy==0x3fe00000) {	/* y is  0.5 */
+		if(hx>=0)	/* x >= +0 */
+		return sqrt(x);	
+	    }
+	}
+
+	ax   = fabs(x);
+    /* special value of x */
+	if(lx==0) {
+	    if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
+		z = ax;			/*x is +-0,+-inf,+-1*/
+		if(hy<0) z = one/z;	/* z = (1/|x|) */
+		if(hx<0) {
+		    if(((ix-0x3ff00000)|yisint)==0) {
+			z = (z-z)/(z-z); /* (-1)**non-int is NaN */
+		    } 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
+       implementation defined.  */
+	n = ((u_int32_t)hx>>31)-1;
+
+    /* (x<0)**(non-int) is NaN */
+	if((n|yisint)==0) return (x-x)/(x-x);
+
+	s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+	if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
+
+    /* |y| is huge */
+	if(iy>0x41e00000) { /* if |y| > 2**31 */
+	    if(iy>0x43f00000){	/* if |y| > 2**64, must o/uflow */
+		if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
+		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));
+	    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;
+	/* take care subnormal number */
+	    if(ix<0x00100000)
+		{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
+	    n  += ((ix)>>20)-0x3ff;
+	    j  = ix&0x000fffff;
+	/* determine interval */
+	    ix = j|0x3ff00000;		/* normalize ix */
+	    if(j<=0x3988E) k=0;		/* |x|<sqrt(3/2) */
+	    else if(j<0xBB67A) k=1;	/* |x|<sqrt(3)   */
+	    else {k=0;n+=1;ix -= 0x00100000;}
+	    SET_HIGH_WORD(ax,ix);
+
+	/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+	    u = ax-bp[k];		/* bp[0]=1.0, bp[1]=1.5 */
+	    v = one/(ax+bp[k]);
+	    ss = u*v;
+	    s_h = ss;
+	    SET_LOW_WORD(s_h,0);
+	/* t_h=ax+bp[k] High */
+	    t_h = zero;
+	    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;
+	    SET_LOW_WORD(t_h,0);
+	    t_l = r-((t_h-3.0)-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;
+	    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);
+	p_l = (y-y1)*t1+y*t2;
+	p_h = y1*t1;
+	z = p_l+p_h;
+	EXTRACT_WORDS(j,i,z);
+	if (j>=0x40900000) {				/* z >= 1024 */
+	    if(((j-0x40900000)|i)!=0)			/* if z > 1024 */
+		return s*huge*huge;			/* overflow */
+	    else {
+		if(p_l+ovt>z-p_h) return s*huge*huge;	/* overflow */
+	    }
+	} else if((j&0x7fffffff)>=0x4090cc00 ) {	/* z <= -1075 */
+	    if(((j-0xc090cc00)|i)!=0) 		/* z < -1075 */
+		return s*tiny*tiny;		/* underflow */
+	    else {
+		if(p_l<=z-p_h) return s*tiny*tiny;	/* underflow */
+	    }
+	}
+    /*
+     * compute 2**(p_h+p_l)
+     */
+	i = j&0x7fffffff;
+	k = (i>>20)-0x3ff;
+	n = 0;
+	if(i>0x3fe00000) {		/* if |z| > 0.5, set n = [z+0.5] */
+	    n = j+(0x00100000>>(k+1));
+	    k = ((n&0x7fffffff)>>20)-0x3ff;	/* new k for n */
+	    t = zero;
+	    SET_HIGH_WORD(t,n&~(0x000fffff>>k));
+	    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;
+	v = (p_l-(t-p_h))*lg2+t*lg2_l;
+	z = u+v;
+	w = v-(z-u);
+	t  = z*z;
+	t1  = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	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;
+}
diff --git a/modules/fdlibm/src/e_sinh.cpp b/modules/fdlibm/src/e_sinh.cpp
new file mode 100644
index 000000000..c3418e687
--- /dev/null
+++ b/modules/fdlibm/src/e_sinh.cpp
@@ -0,0 +1,74 @@
+
+/* @(#)e_sinh.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_sinh(x)
+ * Method : 
+ * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
+ *	1. Replace x by |x| (sinh(-x) = -sinh(x)). 
+ *	2. 
+ *		                                    E + E/(E+1)
+ *	    0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
+ *			       			        2
+ *
+ *	    22       <= x <= lnovft :  sinh(x) := exp(x)/2 
+ *	    lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
+ *	    ln2ovft  <  x	    :  sinh(x) := x*shuge (overflow)
+ *
+ * Special cases:
+ *	sinh(x) is |x| if x is +INF, -INF, or NaN.
+ *	only sinh(0)=0 is exact for finite x.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double one = 1.0, shuge = 1.0e307;
+
+double
+__ieee754_sinh(double x)
+{
+	double t,h;
+	int32_t ix,jx;
+
+    /* High word of |x|. */
+	GET_HIGH_WORD(jx,x);
+	ix = jx&0x7fffffff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) return x+x;	
+
+	h = 0.5;
+	if (jx<0) h = -h;
+    /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
+	if (ix < 0x40360000) {		/* |x|<22 */
+	    if (ix<0x3e300000) 		/* |x|<2**-28 */
+		if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
+	    t = expm1(fabs(x));
+	    if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
+	    return h*(t+t/(t+one));
+	}
+
+    /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
+	if (ix < 0x40862E42)  return h*__ieee754_exp(fabs(x));
+
+    /* |x| in [log(maxdouble), overflowthresold] */
+	if (ix<=0x408633CE)
+	    return h*2.0*__ldexp_exp(fabs(x), -1);
+
+    /* |x| > overflowthresold, sinh(x) overflow */
+	return x*shuge;
+}
diff --git a/modules/fdlibm/src/e_sqrt.cpp b/modules/fdlibm/src/e_sqrt.cpp
new file mode 100644
index 000000000..681505390
--- /dev/null
+++ b/modules/fdlibm/src/e_sqrt.cpp
@@ -0,0 +1,446 @@
+
+/* @(#)e_sqrt.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* __ieee754_sqrt(x)
+ * Return correctly rounded sqrt.
+ *           ------------------------------------------
+ *	     |  Use the hardware sqrt if you have one |
+ *           ------------------------------------------
+ * Method: 
+ *   Bit by bit method using integer arithmetic. (Slow, but portable) 
+ *   1. Normalization
+ *	Scale x to y in [1,4) with even powers of 2: 
+ *	find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+ *		sqrt(x) = 2^k * sqrt(y)
+ *   2. Bit by bit computation
+ *	Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+ *	     i							 0
+ *                                     i+1         2
+ *	    s  = 2*q , and	y  =  2   * ( y - q  ).		(1)
+ *	     i      i            i                 i
+ *                                                        
+ *	To compute q    from q , one checks whether 
+ *		    i+1       i                       
+ *
+ *			      -(i+1) 2
+ *			(q + 2      ) <= y.			(2)
+ *     			  i
+ *							      -(i+1)
+ *	If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+ *		 	       i+1   i             i+1   i
+ *
+ *	With some algebric manipulation, it is not difficult to see
+ *	that (2) is equivalent to 
+ *                             -(i+1)
+ *			s  +  2       <= y			(3)
+ *			 i                i
+ *
+ *	The advantage of (3) is that s  and y  can be computed by 
+ *				      i      i
+ *	the following recurrence formula:
+ *	    if (3) is false
+ *
+ *	    s     =  s  ,	y    = y   ;			(4)
+ *	     i+1      i		 i+1    i
+ *
+ *	    otherwise,
+ *                         -i                     -(i+1)
+ *	    s	  =  s  + 2  ,  y    = y  -  s  - 2  		(5)
+ *           i+1      i          i+1    i     i
+ *				
+ *	One may easily use induction to prove (4) and (5). 
+ *	Note. Since the left hand side of (3) contain only i+2 bits,
+ *	      it does not necessary to do a full (53-bit) comparison 
+ *	      in (3).
+ *   3. Final rounding
+ *	After generating the 53 bits result, we compute one more bit.
+ *	Together with the remainder, we can decide whether the
+ *	result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ *	(it will never equal to 1/2ulp).
+ *	The rounding mode can be detected by checking whether
+ *	huge + tiny is equal to huge, and whether huge - tiny is
+ *	equal to huge for some floating point number "huge" and "tiny".
+ *		
+ * Special cases:
+ *	sqrt(+-0) = +-0 	... exact
+ *	sqrt(inf) = inf
+ *	sqrt(-ve) = NaN		... with invalid signal
+ *	sqrt(NaN) = NaN		... with invalid signal for signaling NaN
+ *
+ * Other methods : see the appended file at the end of the program below.
+ *---------------
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static	const double	one	= 1.0, tiny=1.0e-300;
+
+double
+__ieee754_sqrt(double x)
+{
+	double z;
+	int32_t sign = (int)0x80000000;
+	int32_t ix0,s0,q,m,t,i;
+	u_int32_t r,t1,s1,ix1,q1;
+
+	EXTRACT_WORDS(ix0,ix1,x);
+
+    /* take care of Inf and NaN */
+	if((ix0&0x7ff00000)==0x7ff00000) {			
+	    return x*x+x;		/* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+					   sqrt(-inf)=sNaN */
+	} 
+    /* take care of zero */
+	if(ix0<=0) {
+	    if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
+	    else if(ix0<0)
+		return (x-x)/(x-x);		/* sqrt(-ve) = sNaN */
+	}
+    /* normalize x */
+	m = (ix0>>20);
+	if(m==0) {				/* subnormal x */
+	    while(ix0==0) {
+		m -= 21;
+		ix0 |= (ix1>>11); ix1 <<= 21;
+	    }
+	    for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
+	    m -= i-1;
+	    ix0 |= (ix1>>(32-i));
+	    ix1 <<= i;
+	}
+	m -= 1023;	/* unbias exponent */
+	ix0 = (ix0&0x000fffff)|0x00100000;
+	if(m&1){	/* odd m, double x to make it even */
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	}
+	m >>= 1;	/* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+	ix0 += ix0 + ((ix1&sign)>>31);
+	ix1 += ix1;
+	q = q1 = s0 = s1 = 0;	/* [q,q1] = sqrt(x) */
+	r = 0x00200000;		/* r = moving bit from right to left */
+
+	while(r!=0) {
+	    t = s0+r; 
+	    if(t<=ix0) { 
+		s0   = t+r; 
+		ix0 -= t; 
+		q   += r; 
+	    } 
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	    r>>=1;
+	}
+
+	r = sign;
+	while(r!=0) {
+	    t1 = s1+r; 
+	    t  = s0;
+	    if((t<ix0)||((t==ix0)&&(t1<=ix1))) { 
+		s1  = t1+r;
+		if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
+		ix0 -= t;
+		if (ix1 < t1) ix0 -= 1;
+		ix1 -= t1;
+		q1  += r;
+	    }
+	    ix0 += ix0 + ((ix1&sign)>>31);
+	    ix1 += ix1;
+	    r>>=1;
+	}
+
+    /* use floating add to find out rounding direction */
+	if((ix0|ix1)!=0) {
+	    z = one-tiny; /* trigger inexact flag */
+	    if (z>=one) {
+	        z = one+tiny;
+	        if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
+		else if (z>one) {
+		    if (q1==(u_int32_t)0xfffffffe) q+=1;
+		    q1+=2; 
+		} else
+	            q1 += (q1&1);
+	    }
+	}
+	ix0 = (q>>1)+0x3fe00000;
+	ix1 =  q1>>1;
+	if ((q&1)==1) ix1 |= sign;
+	ix0 += (m <<20);
+	INSERT_WORDS(z,ix0,ix1);
+	return z;
+}
+
+/*
+Other methods  (use floating-point arithmetic)
+-------------
+(This is a copy of a drafted paper by Prof W. Kahan 
+and K.C. Ng, written in May, 1986)
+
+	Two algorithms are given here to implement sqrt(x) 
+	(IEEE double precision arithmetic) in software.
+	Both supply sqrt(x) correctly rounded. The first algorithm (in
+	Section A) uses newton iterations and involves four divisions.
+	The second one uses reciproot iterations to avoid division, but
+	requires more multiplications. Both algorithms need the ability
+	to chop results of arithmetic operations instead of round them, 
+	and the INEXACT flag to indicate when an arithmetic operation
+	is executed exactly with no roundoff error, all part of the 
+	standard (IEEE 754-1985). The ability to perform shift, add,
+	subtract and logical AND operations upon 32-bit words is needed
+	too, though not part of the standard.
+
+A.  sqrt(x) by Newton Iteration
+
+   (1)	Initial approximation
+
+	Let x0 and x1 be the leading and the trailing 32-bit words of
+	a floating point number x (in IEEE double format) respectively 
+
+	    1    11		     52				  ...widths
+	   ------------------------------------------------------
+	x: |s|	  e     |	      f				|
+	   ------------------------------------------------------
+	      msb    lsb  msb				      lsb ...order
+
+ 
+	     ------------------------  	     ------------------------
+	x0:  |s|   e    |    f1     |	 x1: |          f2           |
+	     ------------------------  	     ------------------------
+
+	By performing shifts and subtracts on x0 and x1 (both regarded
+	as integers), we obtain an 8-bit approximation of sqrt(x) as
+	follows.
+
+		k  := (x0>>1) + 0x1ff80000;
+		y0 := k - T1[31&(k>>15)].	... y ~ sqrt(x) to 8 bits
+	Here k is a 32-bit integer and T1[] is an integer array containing
+	correction terms. Now magically the floating value of y (y's
+	leading 32-bit word is y0, the value of its trailing word is 0)
+	approximates sqrt(x) to almost 8-bit.
+
+	Value of T1:
+	static int T1[32]= {
+	0,	1024,	3062,	5746,	9193,	13348,	18162,	23592,
+	29598,	36145,	43202,	50740,	58733,	67158,	75992,	85215,
+	83599,	71378,	60428,	50647,	41945,	34246,	27478,	21581,
+	16499,	12183,	8588,	5674,	3403,	1742,	661,	130,};
+
+    (2)	Iterative refinement
+
+	Apply Heron's rule three times to y, we have y approximates 
+	sqrt(x) to within 1 ulp (Unit in the Last Place):
+
+		y := (y+x/y)/2		... almost 17 sig. bits
+		y := (y+x/y)/2		... almost 35 sig. bits
+		y := y-(y-x/y)/2	... within 1 ulp
+
+
+	Remark 1.
+	    Another way to improve y to within 1 ulp is:
+
+		y := (y+x/y)		... almost 17 sig. bits to 2*sqrt(x)
+		y := y - 0x00100006	... almost 18 sig. bits to sqrt(x)
+
+				2
+			    (x-y )*y
+		y := y + 2* ----------	...within 1 ulp
+			       2
+			     3y  + x
+
+
+	This formula has one division fewer than the one above; however,
+	it requires more multiplications and additions. Also x must be
+	scaled in advance to avoid spurious overflow in evaluating the
+	expression 3y*y+x. Hence it is not recommended uless division
+	is slow. If division is very slow, then one should use the 
+	reciproot algorithm given in section B.
+
+    (3) Final adjustment
+
+	By twiddling y's last bit it is possible to force y to be 
+	correctly rounded according to the prevailing rounding mode
+	as follows. Let r and i be copies of the rounding mode and
+	inexact flag before entering the square root program. Also we
+	use the expression y+-ulp for the next representable floating
+	numbers (up and down) of y. Note that y+-ulp = either fixed
+	point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+	mode.
+
+		I := FALSE;	... reset INEXACT flag I
+		R := RZ;	... set rounding mode to round-toward-zero
+		z := x/y;	... chopped quotient, possibly inexact
+		If(not I) then {	... if the quotient is exact
+		    if(z=y) {
+		        I := i;	 ... restore inexact flag
+		        R := r;  ... restore rounded mode
+		        return sqrt(x):=y.
+		    } else {
+			z := z - ulp;	... special rounding
+		    }
+		}
+		i := TRUE;		... sqrt(x) is inexact
+		If (r=RN) then z=z+ulp	... rounded-to-nearest
+		If (r=RP) then {	... round-toward-+inf
+		    y = y+ulp; z=z+ulp;
+		}
+		y := y+z;		... chopped sum
+		y0:=y0-0x00100000;	... y := y/2 is correctly rounded.
+	        I := i;	 		... restore inexact flag
+	        R := r;  		... restore rounded mode
+	        return sqrt(x):=y.
+		    
+    (4)	Special cases
+
+	Square root of +inf, +-0, or NaN is itself;
+	Square root of a negative number is NaN with invalid signal.
+
+
+B.  sqrt(x) by Reciproot Iteration
+
+   (1)	Initial approximation
+
+	Let x0 and x1 be the leading and the trailing 32-bit words of
+	a floating point number x (in IEEE double format) respectively
+	(see section A). By performing shifs and subtracts on x0 and y0,
+	we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
+
+	    k := 0x5fe80000 - (x0>>1);
+	    y0:= k - T2[63&(k>>14)].	... y ~ 1/sqrt(x) to 7.8 bits
+
+	Here k is a 32-bit integer and T2[] is an integer array 
+	containing correction terms. Now magically the floating
+	value of y (y's leading 32-bit word is y0, the value of
+	its trailing word y1 is set to zero) approximates 1/sqrt(x)
+	to almost 7.8-bit.
+
+	Value of T2:
+	static int T2[64]= {
+	0x1500,	0x2ef8,	0x4d67,	0x6b02,	0x87be,	0xa395,	0xbe7a,	0xd866,
+	0xf14a,	0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
+	0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
+	0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
+	0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
+	0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
+	0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
+	0x1527f,0x1334a,0x11051,0xe951,	0xbe01,	0x8e0d,	0x5924,	0x1edd,};
+
+    (2)	Iterative refinement
+
+	Apply Reciproot iteration three times to y and multiply the
+	result by x to get an approximation z that matches sqrt(x)
+	to about 1 ulp. To be exact, we will have 
+		-1ulp < sqrt(x)-z<1.0625ulp.
+	
+	... set rounding mode to Round-to-nearest
+	   y := y*(1.5-0.5*x*y*y)	... almost 15 sig. bits to 1/sqrt(x)
+	   y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
+	... special arrangement for better accuracy
+	   z := x*y			... 29 bits to sqrt(x), with z*y<1
+	   z := z + 0.5*z*(1-z*y)	... about 1 ulp to sqrt(x)
+
+	Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
+	(a) the term z*y in the final iteration is always less than 1; 
+	(b) the error in the final result is biased upward so that
+		-1 ulp < sqrt(x) - z < 1.0625 ulp
+	    instead of |sqrt(x)-z|<1.03125ulp.
+
+    (3)	Final adjustment
+
+	By twiddling y's last bit it is possible to force y to be 
+	correctly rounded according to the prevailing rounding mode
+	as follows. Let r and i be copies of the rounding mode and
+	inexact flag before entering the square root program. Also we
+	use the expression y+-ulp for the next representable floating
+	numbers (up and down) of y. Note that y+-ulp = either fixed
+	point y+-1, or multiply y by nextafter(1,+-inf) in chopped
+	mode.
+
+	R := RZ;		... set rounding mode to round-toward-zero
+	switch(r) {
+	    case RN:		... round-to-nearest
+	       if(x<= z*(z-ulp)...chopped) z = z - ulp; else
+	       if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
+	       break;
+	    case RZ:case RM:	... round-to-zero or round-to--inf
+	       R:=RP;		... reset rounding mod to round-to-+inf
+	       if(x<z*z ... rounded up) z = z - ulp; else
+	       if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
+	       break;
+	    case RP:		... round-to-+inf
+	       if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
+	       if(x>z*z ...chopped) z = z+ulp;
+	       break;
+	}
+
+	Remark 3. The above comparisons can be done in fixed point. For
+	example, to compare x and w=z*z chopped, it suffices to compare
+	x1 and w1 (the trailing parts of x and w), regarding them as
+	two's complement integers.
+
+	...Is z an exact square root?
+	To determine whether z is an exact square root of x, let z1 be the
+	trailing part of z, and also let x0 and x1 be the leading and
+	trailing parts of x.
+
+	If ((z1&0x03ffffff)!=0)	... not exact if trailing 26 bits of z!=0
+	    I := 1;		... Raise Inexact flag: z is not exact
+	else {
+	    j := 1 - [(x0>>20)&1]	... j = logb(x) mod 2
+	    k := z1 >> 26;		... get z's 25-th and 26-th 
+					    fraction bits
+	    I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
+	}
+	R:= r		... restore rounded mode
+	return sqrt(x):=z.
+
+	If multiplication is cheaper then the foregoing red tape, the 
+	Inexact flag can be evaluated by
+
+	    I := i;
+	    I := (z*z!=x) or I.
+
+	Note that z*z can overwrite I; this value must be sensed if it is 
+	True.
+
+	Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
+	zero.
+
+		    --------------------
+		z1: |        f2        | 
+		    --------------------
+		bit 31		   bit 0
+
+	Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
+	or even of logb(x) have the following relations:
+
+	-------------------------------------------------
+	bit 27,26 of z1		bit 1,0 of x1	logb(x)
+	-------------------------------------------------
+	00			00		odd and even
+	01			01		even
+	10			10		odd
+	10			00		even
+	11			01		even
+	-------------------------------------------------
+
+    (4)	Special cases (see (4) of Section A).	
+ 
+ */
+ 
diff --git a/modules/fdlibm/src/fdlibm.h b/modules/fdlibm/src/fdlibm.h
new file mode 100644
index 000000000..0ad215911
--- /dev/null
+++ b/modules/fdlibm/src/fdlibm.h
@@ -0,0 +1,65 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ * $FreeBSD$
+ */
+
+#ifndef mozilla_imported_fdlibm_h
+#define mozilla_imported_fdlibm_h
+
+namespace fdlibm {
+
+double	acos(double);
+double	asin(double);
+double	atan(double);
+double	atan2(double, double);
+
+double	cosh(double);
+double	sinh(double);
+double	tanh(double);
+
+double	exp(double);
+double	log(double);
+double	log10(double);
+
+double	pow(double, double);
+double	sqrt(double);
+double	fabs(double);
+
+double	floor(double);
+double	trunc(double);
+double	ceil(double);
+
+double	acosh(double);
+double	asinh(double);
+double	atanh(double);
+double	cbrt(double);
+double	expm1(double);
+double	hypot(double, double);
+double	log1p(double);
+double	log2(double);
+double	rint(double);
+double	copysign(double, double);
+double	nearbyint(double);
+double	scalbn(double, int);
+
+float	ceilf(float);
+float	floorf(float);
+
+float	nearbyintf(float);
+float	rintf(float);
+float	truncf(float);
+
+} /* namespace fdlibm */
+
+#endif /* mozilla_imported_fdlibm_h */
diff --git a/modules/fdlibm/src/k_exp.cpp b/modules/fdlibm/src/k_exp.cpp
new file mode 100644
index 000000000..a0699fa4a
--- /dev/null
+++ b/modules/fdlibm/src/k_exp.cpp
@@ -0,0 +1,81 @@
+/*-
+ * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+ #include "math_private.h"
+
+static const uint32_t k = 1799;		/* constant for reduction */
+static const double kln2 =  1246.97177782734161156;	/* k * ln2 */
+
+/*
+ * Compute exp(x), scaled to avoid spurious overflow.  An exponent is
+ * returned separately in 'expt'.
+ *
+ * Input:  ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91
+ * Output: 2**1023 <= y < 2**1024
+ */
+static double
+__frexp_exp(double x, int *expt)
+{
+	double exp_x;
+	uint32_t hx;
+
+	/*
+	 * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to
+	 * minimize |exp(kln2) - 2**k|.  We also scale the exponent of
+	 * exp_x to MAX_EXP so that the result can be multiplied by
+	 * a tiny number without losing accuracy due to denormalization.
+	 */
+	exp_x = exp(x - kln2);
+	GET_HIGH_WORD(hx, exp_x);
+	*expt = (hx >> 20) - (0x3ff + 1023) + k;
+	SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20));
+	return (exp_x);
+}
+
+/*
+ * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt.
+ * They are intended for large arguments (real part >= ln(DBL_MAX))
+ * where care is needed to avoid overflow.
+ *
+ * The present implementation is narrowly tailored for our hyperbolic and
+ * exponential functions.  We assume expt is small (0 or -1), and the caller
+ * has filtered out very large x, for which overflow would be inevitable.
+ */
+
+double
+__ldexp_exp(double x, int expt)
+{
+	double exp_x, scale;
+	int ex_expt;
+
+	exp_x = __frexp_exp(x, &ex_expt);
+	expt += ex_expt;
+	INSERT_WORDS(scale, (0x3ff + expt) << 20, 0);
+	return (exp_x * scale);
+}
diff --git a/modules/fdlibm/src/k_log.h b/modules/fdlibm/src/k_log.h
new file mode 100644
index 000000000..0efa020f6
--- /dev/null
+++ b/modules/fdlibm/src/k_log.h
@@ -0,0 +1,100 @@
+
+/* @(#)e_log.c 1.3 95/01/18 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * k_log1p(f):
+ * Return log(1+f) - f for 1+f in ~[sqrt(2)/2, sqrt(2)].
+ *
+ * The following describes the overall strategy for computing
+ * logarithms in base e.  The argument reduction and adding the final
+ * term of the polynomial are done by the caller for increased accuracy
+ * when different bases are used.
+ *
+ * Method :                  
+ *   1. Argument Reduction: find k and f such that 
+ *			x = 2^k * (1+f), 
+ *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *	     	 = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate 
+ * 	a polynomial of degree 14 to approximate R The maximum error 
+ *	of this polynomial approximation is bounded by 2**-58.45. In
+ *	other words,
+ *		        2      4      6      8      10      12      14
+ *	    R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *  	(the values of Lg1 to Lg7 are listed in the program)
+ *	and
+ *	    |      2          14          |     -58.45
+ *	    | Lg1*s +...+Lg7*s    -  R(z) | <= 2 
+ *	    |                             |
+ *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *	In order to guarantee error in log below 1ulp, we compute log
+ *	by
+ *		log(1+f) = f - s*(f - R)	(if f is not too large)
+ *		log(1+f) = f - (hfsq - s*(hfsq+R)).	(better accuracy)
+ *	
+ *	3. Finally,  log(x) = k*ln2 + log(1+f).  
+ *			    = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *	   Here ln2 is split into two floating point number: 
+ *			ln2_hi + ln2_lo,
+ *	   where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *	log(x) is NaN with signal if x < 0 (including -INF) ; 
+ *	log(+INF) is +INF; log(0) is -INF with signal;
+ *	log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * 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.
+ */
+
+static const double
+Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
+Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
+Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
+Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+
+/*
+ * We always inline k_log1p(), since doing so produces a
+ * substantial performance improvement (~40% on amd64).
+ */
+static inline double
+k_log1p(double f)
+{
+	double hfsq,s,z,R,w,t1,t2;
+
+ 	s = f/(2.0+f);
+	z = s*s;
+	w = z*z;
+	t1= w*(Lg2+w*(Lg4+w*Lg6));
+	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
+	R = t2+t1;
+	hfsq=0.5*f*f;
+	return s*(hfsq+R);
+}
diff --git a/modules/fdlibm/src/math_private.h b/modules/fdlibm/src/math_private.h
new file mode 100644
index 000000000..6947cecc0
--- /dev/null
+++ b/modules/fdlibm/src/math_private.h
@@ -0,0 +1,822 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * from: @(#)fdlibm.h 5.1 93/09/24
+ * $FreeBSD$
+ */
+
+#ifndef _MATH_PRIVATE_H_
+#define	_MATH_PRIVATE_H_
+
+#include <cfloat>
+#include <stdint.h>
+#include <sys/types.h>
+
+#include "fdlibm.h"
+
+#include "mozilla/EndianUtils.h"
+
+/*
+ * 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.
+ * 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.
+ */
+
+#ifdef WIN32
+#define u_int32_t uint32_t
+#define u_int64_t uint64_t
+#endif
+
+/*
+ * 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;
+  } parts;
+  struct
+  {
+    u_int64_t w;
+  } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+#if MOZ_LITTLE_ENDIAN
+
+typedef union
+{
+  double value;
+  struct
+  {
+    u_int32_t lsw;
+    u_int32_t msw;
+  } parts;
+  struct
+  {
+    u_int64_t w;
+  } xparts;
+} ieee_double_shape_type;
+
+#endif
+
+/* Get two 32 bit ints from a double.  */
+
+#define EXTRACT_WORDS(ix0,ix1,d)				\
+do {								\
+  ieee_double_shape_type ew_u;					\
+  ew_u.value = (d);						\
+  (ix0) = ew_u.parts.msw;					\
+  (ix1) = ew_u.parts.lsw;					\
+} while (0)
+
+/* Get a 64-bit int from a double. */
+#define EXTRACT_WORD64(ix,d)					\
+do {								\
+  ieee_double_shape_type ew_u;					\
+  ew_u.value = (d);						\
+  (ix) = ew_u.xparts.w;						\
+} while (0)
+
+/* Get the more significant 32 bit int from a double.  */
+
+#define GET_HIGH_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gh_u;					\
+  gh_u.value = (d);						\
+  (i) = gh_u.parts.msw;						\
+} while (0)
+
+/* Get the less significant 32 bit int from a double.  */
+
+#define GET_LOW_WORD(i,d)					\
+do {								\
+  ieee_double_shape_type gl_u;					\
+  gl_u.value = (d);						\
+  (i) = gl_u.parts.lsw;						\
+} while (0)
+
+/* Set a double from two 32 bit ints.  */
+
+#define INSERT_WORDS(d,ix0,ix1)					\
+do {								\
+  ieee_double_shape_type iw_u;					\
+  iw_u.parts.msw = (ix0);					\
+  iw_u.parts.lsw = (ix1);					\
+  (d) = iw_u.value;						\
+} while (0)
+
+/* Set a double from a 64-bit int. */
+#define INSERT_WORD64(d,ix)					\
+do {								\
+  ieee_double_shape_type iw_u;					\
+  iw_u.xparts.w = (ix);						\
+  (d) = iw_u.value;						\
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int.  */
+
+#define SET_HIGH_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sh_u;					\
+  sh_u.value = (d);						\
+  sh_u.parts.msw = (v);						\
+  (d) = sh_u.value;						\
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int.  */
+
+#define SET_LOW_WORD(d,v)					\
+do {								\
+  ieee_double_shape_type sl_u;					\
+  sl_u.value = (d);						\
+  sl_u.parts.lsw = (v);						\
+  (d) = sl_u.value;						\
+} while (0)
+
+/*
+ * A union which permits us to convert between a float and a 32 bit
+ * int.
+ */
+
+typedef union
+{
+  float value;
+  /* FIXME: Assumes 32 bit int.  */
+  unsigned int word;
+} ieee_float_shape_type;
+
+/* Get a 32 bit int from a float.  */
+
+#define GET_FLOAT_WORD(i,d)					\
+do {								\
+  ieee_float_shape_type gf_u;					\
+  gf_u.value = (d);						\
+  (i) = gf_u.word;						\
+} while (0)
+
+/* Set a float from a 32 bit int.  */
+
+#define SET_FLOAT_WORD(d,i)					\
+do {								\
+  ieee_float_shape_type sf_u;					\
+  sf_u.word = (i);						\
+  (d) = sf_u.value;						\
+} while (0)
+
+/*
+ * Get expsign and mantissa as 16 bit and 64 bit ints from an 80 bit long
+ * double.
+ */
+
+#define	EXTRACT_LDBL80_WORDS(ix0,ix1,d)				\
+do {								\
+  union IEEEl2bits ew_u;					\
+  ew_u.e = (d);							\
+  (ix0) = ew_u.xbits.expsign;					\
+  (ix1) = ew_u.xbits.man;					\
+} while (0)
+
+/*
+ * Get expsign and mantissa as one 16 bit and two 64 bit ints from a 128 bit
+ * long double.
+ */
+
+#define	EXTRACT_LDBL128_WORDS(ix0,ix1,ix2,d)			\
+do {								\
+  union IEEEl2bits ew_u;					\
+  ew_u.e = (d);							\
+  (ix0) = ew_u.xbits.expsign;					\
+  (ix1) = ew_u.xbits.manh;					\
+  (ix2) = ew_u.xbits.manl;					\
+} while (0)
+
+/* Get expsign as a 16 bit int from a long double.  */
+
+#define	GET_LDBL_EXPSIGN(i,d)					\
+do {								\
+  union IEEEl2bits ge_u;					\
+  ge_u.e = (d);							\
+  (i) = ge_u.xbits.expsign;					\
+} while (0)
+
+/*
+ * Set an 80 bit long double from a 16 bit int expsign and a 64 bit int
+ * mantissa.
+ */
+
+#define	INSERT_LDBL80_WORDS(d,ix0,ix1)				\
+do {								\
+  union IEEEl2bits iw_u;					\
+  iw_u.xbits.expsign = (ix0);					\
+  iw_u.xbits.man = (ix1);					\
+  (d) = iw_u.e;							\
+} while (0)
+
+/*
+ * Set a 128 bit long double from a 16 bit int expsign and two 64 bit ints
+ * comprising the mantissa.
+ */
+
+#define	INSERT_LDBL128_WORDS(d,ix0,ix1,ix2)			\
+do {								\
+  union IEEEl2bits iw_u;					\
+  iw_u.xbits.expsign = (ix0);					\
+  iw_u.xbits.manh = (ix1);					\
+  iw_u.xbits.manl = (ix2);					\
+  (d) = iw_u.e;							\
+} while (0)
+
+/* Set expsign of a long double from a 16 bit int.  */
+
+#define	SET_LDBL_EXPSIGN(d,v)					\
+do {								\
+  union IEEEl2bits se_u;					\
+  se_u.e = (d);							\
+  se_u.xbits.expsign = (v);					\
+  (d) = se_u.e;							\
+} while (0)
+
+#ifdef __i386__
+/* Long double constants are broken on i386. */
+#define	LD80C(m, ex, v) {						\
+	.xbits.man = __CONCAT(m, ULL),					\
+	.xbits.expsign = (0x3fff + (ex)) | ((v) < 0 ? 0x8000 : 0),	\
+}
+#else
+/* 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.
+ */
+#if !defined(_MSC_VER) && (FLT_EVAL_METHOD == 0 || __GNUC__ == 0)
+#define	STRICT_ASSIGN(type, lval, rval)	((lval) = (rval))
+#else
+#define	STRICT_ASSIGN(type, lval, rval) do {	\
+	volatile type __lval;			\
+						\
+	if (sizeof(type) >= sizeof(long double))	\
+		(lval) = (rval);		\
+	else {					\
+		__lval = (rval);		\
+		(lval) = __lval;		\
+	}					\
+} while (0)
+#endif
+#else
+#define	STRICT_ASSIGN(type, lval, rval) do {	\
+	volatile type __lval;			\
+						\
+	if (sizeof(type) >= sizeof(long double))	\
+		(lval) = (rval);		\
+	else {					\
+		__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;			\
+	fp_prec_t __oprec;			\
+						\
+	if ((__oprec = fpgetprec()) != FP_PE)	\
+		fpsetprec(FP_PE)
+#define	RETURNI(x) do {				\
+	__retval = (x);				\
+	if (__oprec != FP_PE)			\
+		fpsetprec(__oprec);		\
+	RETURNF(__retval);			\
+} while (0)
+#else
+#define	ENTERI(x)
+#define	RETURNI(x)	RETURNF(x)
+#endif
+
+/* Default return statement if hack*_t() is not used. */
+#define      RETURNF(v)      return (v)
+
+/*
+ * 2sum gives the same result as 2sumF without requiring |a| >= |b| or
+ * a == 0, but is slower.
+ */
+#define	_2sum(a, b) do {	\
+	__typeof(a) __s, __w;	\
+				\
+	__w = (a) + (b);	\
+	__s = __w - (a);	\
+	(b) = ((a) - (__w - __s)) + ((b) - __s); \
+	(a) = __w;		\
+} while (0)
+
+/*
+ * 2sumF algorithm.
+ *
+ * "Normalize" the terms in the infinite-precision expression a + b for
+ * the sum of 2 floating point values so that b is as small as possible
+ * relative to 'a'.  (The resulting 'a' is the value of the expression in
+ * the same precision as 'a' and the resulting b is the rounding error.)
+ * |a| must be >= |b| or 0, b's type must be no larger than 'a's type, and
+ * exponent overflow or underflow must not occur.  This uses a Theorem of
+ * Dekker (1971).  See Knuth (1981) 4.2.2 Theorem C.  The name "TwoSum"
+ * is apparently due to Skewchuk (1997).
+ *
+ * For this to always work, assignment of a + b to 'a' must not retain any
+ * extra precision in a + b.  This is required by C standards but broken
+ * in many compilers.  The brokenness cannot be worked around using
+ * STRICT_ASSIGN() like we do elsewhere, since the efficiency of this
+ * algorithm would be destroyed by non-null strict assignments.  (The
+ * compilers are correct to be broken -- the efficiency of all floating
+ * point code calculations would be destroyed similarly if they forced the
+ * conversions.)
+ *
+ * Fortunately, a case that works well can usually be arranged by building
+ * any extra precision into the type of 'a' -- 'a' should have type float_t,
+ * double_t or long double.  b's type should be no larger than 'a's type.
+ * Callers should use these types with scopes as large as possible, to
+ * reduce their own extra-precision and efficiciency problems.  In
+ * particular, they shouldn't convert back and forth just to call here.
+ */
+#ifdef DEBUG
+#define	_2sumF(a, b) do {				\
+	__typeof(a) __w;				\
+	volatile __typeof(a) __ia, __ib, __r, __vw;	\
+							\
+	__ia = (a);					\
+	__ib = (b);					\
+	assert(__ia == 0 || fabsl(__ia) >= fabsl(__ib));	\
+							\
+	__w = (a) + (b);				\
+	(b) = ((a) - __w) + (b);			\
+	(a) = __w;					\
+							\
+	/* The next 2 assertions are weak if (a) is already long double. */ \
+	assert((long double)__ia + __ib == (long double)(a) + (b));	\
+	__vw = __ia + __ib;				\
+	__r = __ia - __vw;				\
+	__r += __ib;					\
+	assert(__vw == (a) && __r == (b));		\
+} while (0)
+#else /* !DEBUG */
+#define	_2sumF(a, b) do {	\
+	__typeof(a) __w;	\
+				\
+	__w = (a) + (b);	\
+	(b) = ((a) - __w) + (b); \
+	(a) = __w;		\
+} while (0)
+#endif /* DEBUG */
+
+/*
+ * Set x += c, where x is represented in extra precision as a + b.
+ * x must be sufficiently normalized and sufficiently larger than c,
+ * and the result is then sufficiently normalized.
+ *
+ * The details of ordering are that |a| must be >= |c| (so that (a, c)
+ * can be normalized without extra work to swap 'a' with c).  The details of
+ * the normalization are that b must be small relative to the normalized 'a'.
+ * Normalization of (a, c) makes the normalized c tiny relative to the
+ * normalized a, so b remains small relative to 'a' in the result.  However,
+ * b need not ever be tiny relative to 'a'.  For example, b might be about
+ * 2**20 times smaller than 'a' to give about 20 extra bits of precision.
+ * That is usually enough, and adding c (which by normalization is about
+ * 2**53 times smaller than a) cannot change b significantly.  However,
+ * cancellation of 'a' with c in normalization of (a, c) may reduce 'a'
+ * significantly relative to b.  The caller must ensure that significant
+ * cancellation doesn't occur, either by having c of the same sign as 'a',
+ * or by having |c| a few percent smaller than |a|.  Pre-normalization of
+ * (a, b) may help.
+ *
+ * This is is a variant of an algorithm of Kahan (see Knuth (1981) 4.2.2
+ * exercise 19).  We gain considerable efficiency by requiring the terms to
+ * be sufficiently normalized and sufficiently increasing.
+ */
+#define	_3sumF(a, b, c) do {	\
+	__typeof(a) __tmp;	\
+				\
+	__tmp = (c);		\
+	_2sumF(__tmp, (a));	\
+	(b) += (a);		\
+	(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);
+
+#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 {
+	float complex f;
+	float a[2];
+} float_complex;
+typedef union {
+	double complex f;
+	double a[2];
+} double_complex;
+typedef union {
+	long double complex f;
+	long double a[2];
+} long_double_complex;
+#define	REALPART(z)	((z).a[0])
+#define	IMAGPART(z)	((z).a[1])
+
+/*
+ * Inline functions that can be used to construct complex values.
+ *
+ * The C99 standard intends x+I*y to be used for this, but x+I*y is
+ * currently unusable in general since gcc introduces many overflow,
+ * underflow, sign and efficiency bugs by rewriting I*y as
+ * (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
+ * In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
+ * to -0.0+I*0.0.
+ *
+ * The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
+ * to construct complex values.  Compilers that conform to the C99
+ * standard require the following functions to avoid the above issues.
+ */
+
+#ifndef CMPLXF
+static __inline float complex
+CMPLXF(float x, float y)
+{
+	float_complex z;
+
+	REALPART(z) = x;
+	IMAGPART(z) = y;
+	return (z.f);
+}
+#endif
+
+#ifndef CMPLX
+static __inline double complex
+CMPLX(double x, double y)
+{
+	double_complex z;
+
+	REALPART(z) = x;
+	IMAGPART(z) = y;
+	return (z.f);
+}
+#endif
+
+#ifndef CMPLXL
+static __inline long double complex
+CMPLXL(long double x, long double y)
+{
+	long_double_complex z;
+
+	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
+#endif
+
+/* Write a pari script to test things externally. */
+#ifdef DOPRINT
+#include <stdio.h>
+
+#ifndef DOPRINT_SWIZZLE
+#define	DOPRINT_SWIZZLE		0
+#endif
+
+#ifdef DOPRINT_LD80
+
+#define	DOPRINT_START(xp) do {						\
+	uint64_t __lx;							\
+	uint16_t __hx;							\
+									\
+	/* Hack to give more-problematic args. */			\
+	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);				\
+	__lx ^= DOPRINT_SWIZZLE;					\
+	INSERT_LDBL80_WORDS(*xp, __hx, __lx);				\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_D64)
+
+#define	DOPRINT_START(xp) do {						\
+	uint32_t __hx, __lx;						\
+									\
+	EXTRACT_WORDS(__hx, __lx, *xp);					\
+	__lx ^= DOPRINT_SWIZZLE;					\
+	INSERT_WORDS(*xp, __hx, __lx);					\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#elif defined(DOPRINT_F32)
+
+#define	DOPRINT_START(xp) do {						\
+	uint32_t __hx;							\
+									\
+	GET_FLOAT_WORD(__hx, *xp);					\
+	__hx ^= DOPRINT_SWIZZLE;					\
+	SET_FLOAT_WORD(*xp, __hx);					\
+	printf("x = %.21Lg; ", (long double)*xp);			\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
+
+#ifndef DOPRINT_SWIZZLE_HIGH
+#define	DOPRINT_SWIZZLE_HIGH	0
+#endif
+
+#define	DOPRINT_START(xp) do {						\
+	uint64_t __lx, __llx;						\
+	uint16_t __hx;							\
+									\
+	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);			\
+	__llx ^= DOPRINT_SWIZZLE;					\
+	__lx ^= DOPRINT_SWIZZLE_HIGH;					\
+	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);			\
+	printf("x = %.36Lg; ", (long double)*xp);					\
+} while (0)
+#define	DOPRINT_END1(v)							\
+	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
+#define	DOPRINT_END2(hi, lo)						\
+	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",		\
+	    (long double)(hi), (long double)(lo))
+
+#endif /* DOPRINT_LD80 */
+
+#else /* !DOPRINT */
+#define	DOPRINT_START(xp)
+#define	DOPRINT_END1(v)
+#define	DOPRINT_END2(hi, lo)
+#endif /* DOPRINT */
+
+#define	RETURNP(x) do {			\
+	DOPRINT_END1(x);		\
+	RETURNF(x);			\
+} while (0)
+#define	RETURNPI(x) do {		\
+	DOPRINT_END1(x);		\
+	RETURNI(x);			\
+} while (0)
+#define	RETURN2P(x, y) do {		\
+	DOPRINT_END2((x), (y));		\
+	RETURNF((x) + (y));		\
+} while (0)
+#define	RETURN2PI(x, y) do {		\
+	DOPRINT_END2((x), (y));		\
+	RETURNI((x) + (y));		\
+} while (0)
+#ifdef STRUCT_RETURN
+#define	RETURNSP(rp) do {		\
+	if (!(rp)->lo_set)		\
+		RETURNP((rp)->hi);	\
+	RETURN2P((rp)->hi, (rp)->lo);	\
+} while (0)
+#define	RETURNSPI(rp) do {		\
+	if (!(rp)->lo_set)		\
+		RETURNPI((rp)->hi);	\
+	RETURN2PI((rp)->hi, (rp)->lo);	\
+} while (0)
+#endif
+#define	SUM2P(x, y) ({			\
+	const __typeof (x) __x = (x);	\
+	const __typeof (y) __y = (y);	\
+					\
+	DOPRINT_END2(__x, __y);		\
+	__x + __y;			\
+})
+
+/*
+ * ieee style elementary functions
+ *
+ * We rename functions here to improve other sources' diffability
+ * against fdlibm.
+ */
+#define	__ieee754_sqrt	sqrt
+#define	__ieee754_acos	acos
+#define	__ieee754_acosh	acosh
+#define	__ieee754_log	log
+#define	__ieee754_log2	log2
+#define	__ieee754_atanh	atanh
+#define	__ieee754_asin	asin
+#define	__ieee754_atan2	atan2
+#define	__ieee754_exp	exp
+#define	__ieee754_cosh	cosh
+#define	__ieee754_fmod	fmod
+#define	__ieee754_pow	pow
+#define	__ieee754_lgamma lgamma
+#define	__ieee754_gamma	gamma
+#define	__ieee754_lgamma_r lgamma_r
+#define	__ieee754_gamma_r gamma_r
+#define	__ieee754_log10	log10
+#define	__ieee754_sinh	sinh
+#define	__ieee754_hypot	hypot
+#define	__ieee754_j0	j0
+#define	__ieee754_j1	j1
+#define	__ieee754_y0	y0
+#define	__ieee754_y1	y1
+#define	__ieee754_jn	jn
+#define	__ieee754_yn	yn
+#define	__ieee754_remainder remainder
+#define	__ieee754_scalb	scalb
+#define	__ieee754_sqrtf	sqrtf
+#define	__ieee754_acosf	acosf
+#define	__ieee754_acoshf acoshf
+#define	__ieee754_logf	logf
+#define	__ieee754_atanhf atanhf
+#define	__ieee754_asinf	asinf
+#define	__ieee754_atan2f atan2f
+#define	__ieee754_expf	expf
+#define	__ieee754_coshf	coshf
+#define	__ieee754_fmodf	fmodf
+#define	__ieee754_powf	powf
+#define	__ieee754_lgammaf lgammaf
+#define	__ieee754_gammaf gammaf
+#define	__ieee754_lgammaf_r lgammaf_r
+#define	__ieee754_gammaf_r gammaf_r
+#define	__ieee754_log10f log10f
+#define	__ieee754_log2f log2f
+#define	__ieee754_sinhf	sinhf
+#define	__ieee754_hypotf hypotf
+#define	__ieee754_j0f	j0f
+#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
+
+#define acos fdlibm::acos
+#define asin fdlibm::asin
+#define atan fdlibm::atan
+#define atan2 fdlibm::atan2
+#define cosh fdlibm::cosh
+#define sinh fdlibm::sinh
+#define tanh fdlibm::tanh
+#define exp fdlibm::exp
+#define log fdlibm::log
+#define log10 fdlibm::log10
+#define pow fdlibm::pow
+#define sqrt fdlibm::sqrt
+#define ceil fdlibm::ceil
+#define ceilf fdlibm::ceilf
+#define fabs fdlibm::fabs
+#define floor fdlibm::floor
+#define acosh fdlibm::acosh
+#define asinh fdlibm::asinh
+#define atanh fdlibm::atanh
+#define cbrt fdlibm::cbrt
+#define expm1 fdlibm::expm1
+#define hypot fdlibm::hypot
+#define log1p fdlibm::log1p
+#define log2 fdlibm::log2
+#define scalb fdlibm::scalb
+#define copysign fdlibm::copysign
+#define scalbn fdlibm::scalbn
+#define trunc fdlibm::trunc
+#define truncf fdlibm::truncf
+#define floorf fdlibm::floorf
+#define nearbyint fdlibm::nearbyint
+#define nearbyintf fdlibm::nearbyintf
+#define rint fdlibm::rint
+#define rintf fdlibm::rintf
+
+/* fdlibm kernel function */
+int	__kernel_rem_pio2(double*,double*,int,int,int);
+
+/* double precision kernel functions */
+#ifndef INLINE_REM_PIO2
+int	__ieee754_rem_pio2(double,double*);
+#endif
+double	__kernel_sin(double,double,int);
+double	__kernel_cos(double,double);
+double	__kernel_tan(double,double,int);
+double	__ldexp_exp(double,int);
+#ifdef _COMPLEX_H
+double complex __ldexp_cexp(double complex,int);
+#endif
+
+/* float precision kernel functions */
+#ifndef INLINE_REM_PIO2F
+int	__ieee754_rem_pio2f(float,double*);
+#endif
+#ifndef INLINE_KERNEL_SINDF
+float	__kernel_sindf(double);
+#endif
+#ifndef INLINE_KERNEL_COSDF
+float	__kernel_cosdf(double);
+#endif
+#ifndef INLINE_KERNEL_TANDF
+float	__kernel_tandf(double,int);
+#endif
+float	__ldexp_expf(float,int);
+#ifdef _COMPLEX_H
+float complex __ldexp_cexpf(float complex,int);
+#endif
+
+/* long double precision kernel functions */
+long double __kernel_sinl(long double, long double, int);
+long double __kernel_cosl(long double, long double);
+long double __kernel_tanl(long double, long double, int);
+
+#endif /* !_MATH_PRIVATE_H_ */
diff --git a/modules/fdlibm/src/moz.build b/modules/fdlibm/src/moz.build
new file mode 100644
index 000000000..b197881ac
--- /dev/null
+++ b/modules/fdlibm/src/moz.build
@@ -0,0 +1,67 @@
+# -*- Mode: python; indent-tabs-mode: nil; tab-width: 40 -*-
+# vim: set filetype=python:
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+EXPORTS += [
+    'fdlibm.h',
+]
+
+FINAL_LIBRARY = 'js'
+
+if CONFIG['GNU_CXX']:
+    CXXFLAGS += [
+        '-Wno-parentheses',
+        '-Wno-sign-compare',
+    ]
+
+if CONFIG['CLANG_CXX']:
+    CXXFLAGS += [
+        '-Wno-dangling-else',
+    ]
+
+if CONFIG['_MSC_VER']:
+    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
+    ]
+
+SOURCES += [
+    'e_acos.cpp',
+    'e_acosh.cpp',
+    'e_asin.cpp',
+    'e_atan2.cpp',
+    'e_atanh.cpp',
+    'e_cosh.cpp',
+    'e_exp.cpp',
+    'e_hypot.cpp',
+    'e_log.cpp',
+    'e_log10.cpp',
+    'e_log2.cpp',
+    'e_pow.cpp',
+    'e_sinh.cpp',
+    'e_sqrt.cpp',
+    'k_exp.cpp',
+    's_asinh.cpp',
+    's_atan.cpp',
+    's_cbrt.cpp',
+    's_ceil.cpp',
+    's_ceilf.cpp',
+    's_copysign.cpp',
+    's_expm1.cpp',
+    's_fabs.cpp',
+    's_floor.cpp',
+    's_floorf.cpp',
+    's_log1p.cpp',
+    's_nearbyint.cpp',
+    's_rint.cpp',
+    's_rintf.cpp',
+    's_scalbn.cpp',
+    's_tanh.cpp',
+    's_trunc.cpp',
+    's_truncf.cpp',
+]
diff --git a/modules/fdlibm/src/s_asinh.cpp b/modules/fdlibm/src/s_asinh.cpp
new file mode 100644
index 000000000..400fd89c1
--- /dev/null
+++ b/modules/fdlibm/src/s_asinh.cpp
@@ -0,0 +1,57 @@
+/* @(#)s_asinh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* asinh(x)
+ * Method :
+ *	Based on
+ *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
+ *	we have
+ *	asinh(x) := x  if  1+x*x=1,
+ *		 := sign(x)*(log(x)+ln2)) for large |x|, else
+ *		 := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
+ *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
+huge=  1.00000000000000000000e+300;
+
+double
+asinh(double x)
+{
+	double t,w;
+	int32_t hx,ix;
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x7ff00000) return x+x;	/* x is inf or NaN */
+	if(ix< 0x3e300000) {	/* |x|<2**-28 */
+	    if(huge+x>one) return x;	/* return x inexact except 0 */
+	}
+	if(ix>0x41b00000) {	/* |x| > 2**28 */
+	    w = __ieee754_log(fabs(x))+ln2;
+	} else if (ix>0x40000000) {	/* 2**28 > |x| > 2.0 */
+	    t = fabs(x);
+	    w = __ieee754_log(2.0*t+one/(__ieee754_sqrt(x*x+one)+t));
+	} else {		/* 2.0 > |x| > 2**-28 */
+	    t = x*x;
+	    w =log1p(fabs(x)+t/(one+__ieee754_sqrt(one+t)));
+	}
+	if(hx>0) return w; else return -w;
+}
diff --git a/modules/fdlibm/src/s_atan.cpp b/modules/fdlibm/src/s_atan.cpp
new file mode 100644
index 000000000..21bc0d820
--- /dev/null
+++ b/modules/fdlibm/src/s_atan.cpp
@@ -0,0 +1,119 @@
+/* @(#)s_atan.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* atan(x)
+ * Method
+ *   1. Reduce x to positive by atan(x) = -atan(-x).
+ *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ *      is further reduced to one of the following intervals and the
+ *      arctangent of t is evaluated by the corresponding formula:
+ *
+ *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
+ *
+ * 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_private.h"
+
+static const double atanhi[] = {
+  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 */
+};
+
+static const double atanlo[] = {
+  2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+  3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+  1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+  6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+};
+
+static const double aT[] = {
+  3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
+ -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+  1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
+ -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+  9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
+ -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+  6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
+ -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+  4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
+ -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+  1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
+};
+
+	static const double
+one   = 1.0,
+huge   = 1.0e300;
+
+double
+atan(double x)
+{
+	double w,s1,s2,z;
+	int32_t ix,hx,id;
+
+	GET_HIGH_WORD(hx,x);
+	ix = hx&0x7fffffff;
+	if(ix>=0x44100000) {	/* if |x| >= 2^66 */
+	    u_int32_t low;
+	    GET_LOW_WORD(low,x);
+	    if(ix>0x7ff00000||
+		(ix==0x7ff00000&&(low!=0)))
+		return x+x;		/* NaN */
+	    if(hx>0) return  atanhi[3]+*(volatile double *)&atanlo[3];
+	    else     return -atanhi[3]-*(volatile double *)&atanlo[3];
+	} if (ix < 0x3fdc0000) {	/* |x| < 0.4375 */
+	    if (ix < 0x3e400000) {	/* |x| < 2^-27 */
+		if(huge+x>one) return x;	/* raise inexact */
+	    }
+	    id = -1;
+	} else {
+	x = fabs(x);
+	if (ix < 0x3ff30000) {		/* |x| < 1.1875 */
+	    if (ix < 0x3fe60000) {	/* 7/16 <=|x|<11/16 */
+		id = 0; x = (2.0*x-one)/(2.0+x);
+	    } else {			/* 11/16<=|x|< 19/16 */
+		id = 1; x  = (x-one)/(x+one);
+	    }
+	} else {
+	    if (ix < 0x40038000) {	/* |x| < 2.4375 */
+		id = 2; x  = (x-1.5)/(one+1.5*x);
+	    } else {			/* 2.4375 <= |x| < 2^66 */
+		id = 3; x  = -1.0/x;
+	    }
+	}}
+    /* end of argument reduction */
+	z = x*x;
+	w = z*z;
+    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+	s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
+	s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
+	if (id<0) return x - x*(s1+s2);
+	else {
+	    z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+	    return (hx<0)? -z:z;
+	}
+}
diff --git a/modules/fdlibm/src/s_cbrt.cpp b/modules/fdlibm/src/s_cbrt.cpp
new file mode 100644
index 000000000..a2de24af7
--- /dev/null
+++ b/modules/fdlibm/src/s_cbrt.cpp
@@ -0,0 +1,112 @@
+/* @(#)s_cbrt.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+#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 */
+
+/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
+static const double
+P0 =  1.87595182427177009643,		/* 0x3ffe03e6, 0x0f61e692 */
+P1 = -1.88497979543377169875,		/* 0xbffe28e0, 0x92f02420 */
+P2 =  1.621429720105354466140,		/* 0x3ff9f160, 0x4a49d6c2 */
+P3 = -0.758397934778766047437,		/* 0xbfe844cb, 0xbee751d9 */
+P4 =  0.145996192886612446982;		/* 0x3fc2b000, 0xd4e4edd7 */
+
+double
+cbrt(double x)
+{
+	int32_t	hx;
+	union {
+	    double value;
+	    uint64_t bits;
+	} u;
+	double r,s,t=0.0,w;
+	u_int32_t sign;
+	u_int32_t high,low;
+
+	EXTRACT_WORDS(hx,low,x);
+	sign=hx&0x80000000; 		/* sign= sign(x) */
+	hx  ^=sign;
+	if(hx>=0x7ff00000) return(x+x); /* cbrt(NaN,INF) is itself */
+
+    /*
+     * Rough cbrt to 5 bits:
+     *    cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
+     * where e is integral and >= 0, m is real and in [0, 1), and "/" and
+     * "%" are integer division and modulus with rounding towards minus
+     * infinity.  The RHS is always >= the LHS and has a maximum relative
+     * error of about 1 in 16.  Adding a bias of -0.03306235651 to the
+     * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
+     * floating point representation, for finite positive normal values,
+     * ordinary integer division of the value in bits magically gives
+     * almost exactly the RHS of the above provided we first subtract the
+     * exponent bias (1023 for doubles) and later add it back.  We do the
+     * subtraction virtually to keep e >= 0 so that ordinary integer
+     * division rounds towards minus infinity; this is also efficient.
+     */
+	if(hx<0x00100000) { 		/* zero or subnormal? */
+	    if((hx|low)==0)
+		return(x);		/* cbrt(0) is itself */
+	    SET_HIGH_WORD(t,0x43500000); /* set t= 2**54 */
+	    t*=x;
+	    GET_HIGH_WORD(high,t);
+	    INSERT_WORDS(t,sign|((high&0x7fffffff)/3+B2),0);
+	} else
+	    INSERT_WORDS(t,sign|(hx/3+B1),0);
+
+    /*
+     * New cbrt to 23 bits:
+     *    cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
+     * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
+     * to within 2**-23.5 when |r - 1| < 1/10.  The rough approximation
+     * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
+     * gives us bounds for r = t**3/x.
+     *
+     * Try to optimize for parallel evaluation as in k_tanf.c.
+     */
+	r=(t*t)*(t/x);
+	t=t*((P0+r*(P1+r*P2))+((r*r)*r)*(P3+r*P4));
+
+    /*
+     * Round t away from zero to 23 bits (sloppily except for ensuring that
+     * the result is larger in magnitude than cbrt(x) but not much more than
+     * 2 23-bit ulps larger).  With rounding towards zero, the error bound
+     * would be ~5/6 instead of ~4/6.  With a maximum error of 2 23-bit ulps
+     * in the rounded t, the infinite-precision error in the Newton
+     * approximation barely affects third digit in the final error
+     * 0.667; the error in the rounded t can be up to about 3 23-bit ulps
+     * before the final error is larger than 0.667 ulps.
+     */
+	u.value=t;
+	u.bits=(u.bits+0x80000000)&0xffffffffc0000000ULL;
+	t=u.value;
+
+    /* one step Newton iteration to 53 bits with error < 0.667 ulps */
+	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);
+}
diff --git a/modules/fdlibm/src/s_ceil.cpp b/modules/fdlibm/src/s_ceil.cpp
new file mode 100644
index 000000000..67e9c1679
--- /dev/null
+++ b/modules/fdlibm/src/s_ceil.cpp
@@ -0,0 +1,72 @@
+/* @(#)s_ceil.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * ceil(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to ceil(x).
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+ceil(double x)
+{
+	int32_t i0,i1,j0;
+	u_int32_t i,j;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0<0) {i0=0x80000000;i1=0;}
+		    else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;}
+		}
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0>0) i0 += (0x00100000)>>j0;
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0) { 		/* raise inexact flag */
+		if(i0>0) {
+		    if(j0==20) i0+=1;
+		    else {
+			j = i1 + (1<<(52-j0));
+			if(j<i1) i0+=1;	/* got a carry */
+			i1 = j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	INSERT_WORDS(x,i0,i1);
+	return x;
+}
diff --git a/modules/fdlibm/src/s_ceilf.cpp b/modules/fdlibm/src/s_ceilf.cpp
new file mode 100644
index 000000000..7b52deeed
--- /dev/null
+++ b/modules/fdlibm/src/s_ceilf.cpp
@@ -0,0 +1,51 @@
+/* s_ceilf.c -- float version of s_ceil.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+#include "math_private.h"
+
+static const float huge = 1.0e30;
+
+float
+ceilf(float x)
+{
+	int32_t i0,j0;
+	u_int32_t i;
+
+	GET_FLOAT_WORD(i0,x);
+	j0 = ((i0>>23)&0xff)-0x7f;
+	if(j0<23) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0<0) {i0=0x80000000;}
+		    else if(i0!=0) { i0=0x3f800000;}
+		}
+	    } else {
+		i = (0x007fffff)>>j0;
+		if((i0&i)==0) return x; /* x is integral */
+		if(huge+x>(float)0.0) {	/* raise inexact flag */
+		    if(i0>0) i0 += (0x00800000)>>j0;
+		    i0 &= (~i);
+		}
+	    }
+	} else {
+	    if(j0==0x80) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	}
+	SET_FLOAT_WORD(x,i0);
+	return x;
+}
diff --git a/modules/fdlibm/src/s_copysign.cpp b/modules/fdlibm/src/s_copysign.cpp
new file mode 100644
index 000000000..b150106fb
--- /dev/null
+++ b/modules/fdlibm/src/s_copysign.cpp
@@ -0,0 +1,32 @@
+/* @(#)s_copysign.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * copysign(double x, double y)
+ * copysign(x,y) returns a value with the magnitude of x and
+ * with the sign bit of y.
+ */
+
+#include "math_private.h"
+
+double
+copysign(double x, double y)
+{
+	u_int32_t hx,hy;
+	GET_HIGH_WORD(hx,x);
+	GET_HIGH_WORD(hy,y);
+	SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
+        return x;
+}
diff --git a/modules/fdlibm/src/s_expm1.cpp b/modules/fdlibm/src/s_expm1.cpp
new file mode 100644
index 000000000..4c19485de
--- /dev/null
+++ b/modules/fdlibm/src/s_expm1.cpp
@@ -0,0 +1,220 @@
+/* @(#)s_expm1.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* expm1(x)
+ * Returns exp(x)-1, the exponential of x minus 1.
+ *
+ * Method
+ *   1. Argument reduction:
+ *	Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2 ~ 0.34658
+ *
+ *      Here a correction term c will be computed to compensate
+ *	the error in r when rounded to a floating-point number.
+ *
+ *   2. Approximating expm1(r) by a special rational function on
+ *	the interval [0,0.34658]:
+ *	Since
+ *	    r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+ *	we define R1(r*r) by
+ *	    r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+ *	That is,
+ *	    R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+ *		     = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+ *		     = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+ *      We use a special Reme algorithm on [0,0.347] to generate
+ * 	a polynomial of degree 5 in r*r to approximate R1. The
+ *	maximum error of this polynomial approximation is bounded
+ *	by 2**-61. In other words,
+ *	    R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+ *	where 	Q1  =  -1.6666666666666567384E-2,
+ * 		Q2  =   3.9682539681370365873E-4,
+ * 		Q3  =  -9.9206344733435987357E-6,
+ * 		Q4  =   2.5051361420808517002E-7,
+ * 		Q5  =  -6.2843505682382617102E-9;
+ *		z   =  r*r,
+ *	with error bounded by
+ *	    |                  5           |     -61
+ *	    | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
+ *	    |                              |
+ *
+ *	expm1(r) = exp(r)-1 is then computed by the following
+ * 	specific way which minimize the accumulation rounding error:
+ *			       2     3
+ *			      r     r    [ 3 - (R1 + R1*r/2)  ]
+ *	      expm1(r) = r + --- + --- * [--------------------]
+ *		              2     2    [ 6 - r*(3 - R1*r/2) ]
+ *
+ *	To compensate the error in the argument reduction, we use
+ *		expm1(r+c) = expm1(r) + c + expm1(r)*c
+ *			   ~ expm1(r) + c + r*c
+ *	Thus c+r*c will be added in as the correction terms for
+ *	expm1(r+c). Now rearrange the term to avoid optimization
+ * 	screw up:
+ *		        (      2                                    2 )
+ *		        ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
+ *	 expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+ *	                ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
+ *                      (                                             )
+ *
+ *		   = r - E
+ *   3. Scale back to obtain expm1(x):
+ *	From step 1, we have
+ *	   expm1(x) = either 2^k*[expm1(r)+1] - 1
+ *		    = or     2^k*[expm1(r) + (1-2^-k)]
+ *   4. Implementation notes:
+ *	(A). To save one multiplication, we scale the coefficient Qi
+ *	     to Qi*2^i, and replace z by (x^2)/2.
+ *	(B). To achieve maximum accuracy, we compute expm1(x) by
+ *	  (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+ *	  (ii)  if k=0, return r-E
+ *	  (iii) if k=-1, return 0.5*(r-E)-0.5
+ *        (iv)	if k=1 if r < -0.25, return 2*((r+0.5)- E)
+ *	       	       else	     return  1.0+2.0*(r-E);
+ *	  (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+ *	  (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+ *	  (vii) return 2^k(1-((E+2^-k)-r))
+ *
+ * Special cases:
+ *	expm1(INF) is INF, expm1(NaN) is NaN;
+ *	expm1(-INF) is -1, and
+ *	for finite argument, only expm1(0)=0 is exact.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *	For IEEE double
+ *	    if x >  7.09782712893383973096e+02 then expm1(x) overflow
+ *
+ * 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_private.h"
+
+static const double
+one		= 1.0,
+tiny		= 1.0e-300,
+o_threshold	= 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
+ln2_hi		= 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
+ln2_lo		= 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
+invln2		= 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
+/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
+Q1  =  -3.33333333333331316428e-02, /* BFA11111 111110F4 */
+Q2  =   1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
+Q3  =  -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
+Q4  =   4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */
+Q5  =  -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
+
+static volatile double huge = 1.0e+300;
+
+double
+expm1(double x)
+{
+	double y,hi,lo,c,t,e,hxs,hfx,r1,twopk;
+	int32_t k,xsb;
+	u_int32_t hx;
+
+	GET_HIGH_WORD(hx,x);
+	xsb = hx&0x80000000;		/* sign bit of x */
+	hx &= 0x7fffffff;		/* high word of |x| */
+
+    /* filter out huge and non-finite argument */
+	if(hx >= 0x4043687A) {			/* if |x|>=56*ln2 */
+	    if(hx >= 0x40862E42) {		/* if |x|>=709.78... */
+                if(hx>=0x7ff00000) {
+		    u_int32_t low;
+		    GET_LOW_WORD(low,x);
+		    if(((hx&0xfffff)|low)!=0)
+		         return x+x; 	 /* NaN */
+		    else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
+	        }
+	        if(x > o_threshold) return huge*huge; /* overflow */
+	    }
+	    if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */
+		if(x+tiny<0.0)		/* raise inexact */
+		return tiny-one;	/* return -1 */
+	    }
+	}
+
+    /* argument reduction */
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */
+	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
+		if(xsb==0)
+		    {hi = x - ln2_hi; lo =  ln2_lo;  k =  1;}
+		else
+		    {hi = x + ln2_hi; lo = -ln2_lo;  k = -1;}
+	    } else {
+		k  = invln2*x+((xsb==0)?0.5:-0.5);
+		t  = k;
+		hi = x - t*ln2_hi;	/* t*ln2_hi is exact here */
+		lo = t*ln2_lo;
+	    }
+	    STRICT_ASSIGN(double, x, hi - lo);
+	    c  = (hi-x)-lo;
+	}
+	else if(hx < 0x3c900000) {  	/* when |x|<2**-54, return x */
+	    t = huge+x;	/* return x with inexact flags when x!=0 */
+	    return x - (t-(huge+x));
+	}
+	else k = 0;
+
+    /* 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 */
+	    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 */
+	        y = one-(e-x);
+		if (k == 1024) {
+		    double const_0x1p1023 = pow(2, 1023);
+		    y = y*2.0*const_0x1p1023;
+		}
+		else y = y*twopk;
+	        return y-one;
+	    }
+	    t = one;
+	    if(k<20) {
+	        SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k));  /* t=1-2^-k */
+	       	y = t-(e-x);
+		y = y*twopk;
+	   } else {
+		SET_HIGH_WORD(t,((0x3ff-k)<<20));	/* 2^-k */
+	       	y = x-(e+t);
+	       	y += one;
+		y = y*twopk;
+	    }
+	}
+	return y;
+}
diff --git a/modules/fdlibm/src/s_fabs.cpp b/modules/fdlibm/src/s_fabs.cpp
new file mode 100644
index 000000000..3bea0478a
--- /dev/null
+++ b/modules/fdlibm/src/s_fabs.cpp
@@ -0,0 +1,30 @@
+ /* @(#)s_fabs.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+ //static char rcsid[] = "$FreeBSD$";
+#endif
+
+/*
+ * fabs(x) returns the absolute value of x.
+ */
+
+#include "math_private.h"
+
+double
+fabs(double x)
+{
+	u_int32_t high;
+	GET_HIGH_WORD(high,x);
+	SET_HIGH_WORD(x,high&0x7fffffff);
+        return x;
+}
diff --git a/modules/fdlibm/src/s_floor.cpp b/modules/fdlibm/src/s_floor.cpp
new file mode 100644
index 000000000..da57fc828
--- /dev/null
+++ b/modules/fdlibm/src/s_floor.cpp
@@ -0,0 +1,73 @@
+/* @(#)s_floor.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * floor(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to floor(x).
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+floor(double x)
+{
+	int32_t i0,i1,j0;
+	u_int32_t i,j;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0>=0) {i0=i1=0;}
+		    else if(((i0&0x7fffffff)|i1)!=0)
+			{ i0=0xbff00000;i1=0;}
+		}
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    if(i0<0) i0 += (0x00100000)>>j0;
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0) { 		/* raise inexact flag */
+		if(i0<0) {
+		    if(j0==20) i0+=1;
+		    else {
+			j = i1+(1<<(52-j0));
+			if(j<i1) i0 +=1 ; 	/* got a carry */
+			i1=j;
+		    }
+		}
+		i1 &= (~i);
+	    }
+	}
+	INSERT_WORDS(x,i0,i1);
+	return x;
+}
diff --git a/modules/fdlibm/src/s_floorf.cpp b/modules/fdlibm/src/s_floorf.cpp
new file mode 100644
index 000000000..88511f209
--- /dev/null
+++ b/modules/fdlibm/src/s_floorf.cpp
@@ -0,0 +1,60 @@
+/* s_floorf.c -- float version of s_floor.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * floorf(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to floorf(x).
+ */
+
+#include "math_private.h"
+
+static const float huge = 1.0e30;
+
+float
+floorf(float x)
+{
+	int32_t i0,j0;
+	u_int32_t i;
+	GET_FLOAT_WORD(i0,x);
+	j0 = ((i0>>23)&0xff)-0x7f;
+	if(j0<23) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>(float)0.0) {/* return 0*sign(x) if |x|<1 */
+		    if(i0>=0) {i0=0;}
+		    else if((i0&0x7fffffff)!=0)
+			{ i0=0xbf800000;}
+		}
+	    } else {
+		i = (0x007fffff)>>j0;
+		if((i0&i)==0) return x; /* x is integral */
+		if(huge+x>(float)0.0) {	/* raise inexact flag */
+		    if(i0<0) i0 += (0x00800000)>>j0;
+		    i0 &= (~i);
+		}
+	    }
+	} else {
+	    if(j0==0x80) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	}
+	SET_FLOAT_WORD(x,i0);
+	return x;
+}
diff --git a/modules/fdlibm/src/s_log1p.cpp b/modules/fdlibm/src/s_log1p.cpp
new file mode 100644
index 000000000..afc6919c6
--- /dev/null
+++ b/modules/fdlibm/src/s_log1p.cpp
@@ -0,0 +1,175 @@
+/* @(#)s_log1p.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* double log1p(double x)
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *			1+x = 2^k * (1+f),
+ *	   where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *      Note. If k=0, then f=x is exact. However, if k!=0, then f
+ *	may not be representable exactly. In that case, a correction
+ *	term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+ *	log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+ *	and add back the correction term c/u.
+ *	(Note: when x > 2**53, one can simply return log(x))
+ *
+ *   2. Approximation of log1p(f).
+ *	Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *		 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *	     	 = 2s + s*R
+ *      We use a special Reme algorithm on [0,0.1716] to generate
+ * 	a polynomial of degree 14 to approximate R The maximum error
+ *	of this polynomial approximation is bounded by 2**-58.45. In
+ *	other words,
+ *		        2      4      6      8      10      12      14
+ *	    R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
+ *  	(the values of Lp1 to Lp7 are listed in the program)
+ *	and
+ *	    |      2          14          |     -58.45
+ *	    | Lp1*s +...+Lp7*s    -  R(z) | <= 2
+ *	    |                             |
+ *	Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *	In order to guarantee error in log below 1ulp, we compute log
+ *	by
+ *		log1p(f) = f - (hfsq - s*(hfsq+R)).
+ *
+ *	3. Finally, log1p(x) = k*ln2 + log1p(f).
+ *		 	     = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *	   Here ln2 is split into two floating point number:
+ *			ln2_hi + ln2_lo,
+ *	   where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *	log1p(x) is NaN with signal if x < -1 (including -INF) ;
+ *	log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+ *	log1p(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * 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.
+ *
+ * Note: Assuming log() return accurate answer, the following
+ * 	 algorithm can be used to compute log1p(x) to within a few ULP:
+ *
+ *		u = 1+x;
+ *		if(u==1.0) return x ; else
+ *			   return log(u)*(x/(u-1.0));
+ *
+ *	 See HP-15C Advanced Functions Handbook, p.193.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
+ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
+two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
+Lp1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
+Lp2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
+Lp3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
+Lp4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
+Lp5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
+Lp6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
+Lp7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */
+
+static const double zero = 0.0;
+static volatile double vzero = 0.0;
+
+double
+log1p(double x)
+{
+	double hfsq,f,c,s,z,R,u;
+	int32_t k,hx,hu,ax;
+
+	GET_HIGH_WORD(hx,x);
+	ax = hx&0x7fffffff;
+
+	k = 1;
+	if (hx < 0x3FDA827A) {			/* 1+x < sqrt(2)+ */
+	    if(ax>=0x3ff00000) {		/* x <= -1.0 */
+		if(x==-1.0) return -two54/vzero; /* log1p(-1)=+inf */
+		else return (x-x)/(x-x);	/* log1p(x<-1)=NaN */
+	    }
+	    if(ax<0x3e200000) {			/* |x| < 2**-29 */
+		if(two54+x>zero			/* raise inexact */
+	            &&ax<0x3c900000) 		/* |x| < 2**-54 */
+		    return x;
+		else
+		    return x - x*x*0.5;
+	    }
+	    if(hx>0||hx<=((int32_t)0xbfd2bec4)) {
+		k=0;f=x;hu=1;}		/* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+	}
+	if (hx >= 0x7ff00000) return x+x;
+	if(k!=0) {
+	    if(hx<0x43400000) {
+		STRICT_ASSIGN(double,u,1.0+x);
+		GET_HIGH_WORD(hu,u);
+	        k  = (hu>>20)-1023;
+	        c  = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
+		c /= u;
+	    } else {
+		u  = x;
+		GET_HIGH_WORD(hu,u);
+	        k  = (hu>>20)-1023;
+		c  = 0;
+	    }
+	    hu &= 0x000fffff;
+	    /*
+	     * The approximation to sqrt(2) used in thresholds is not
+	     * critical.  However, the ones used above must give less
+	     * strict bounds than the one here so that the k==0 case is
+	     * never reached from here, since here we have committed to
+	     * using the correction term but don't use it if k==0.
+	     */
+	    if(hu<0x6a09e) {			/* u ~< sqrt(2) */
+	        SET_HIGH_WORD(u,hu|0x3ff00000);	/* normalize u */
+	    } else {
+	        k += 1;
+		SET_HIGH_WORD(u,hu|0x3fe00000);	/* normalize u/2 */
+	        hu = (0x00100000-hu)>>2;
+	    }
+	    f = u-1.0;
+	}
+	hfsq=0.5*f*f;
+	if(hu==0) {	/* |f| < 2**-20 */
+	    if(f==zero) {
+		if(k==0) {
+		    return zero;
+		} else {
+		    c += k*ln2_lo;
+		    return k*ln2_hi+c;
+		}
+	    }
+	    R = hfsq*(1.0-0.66666666666666666*f);
+	    if(k==0) return f-R; else
+	    	     return k*ln2_hi-((R-(k*ln2_lo+c))-f);
+	}
+ 	s = f/(2.0+f);
+	z = s*s;
+	R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
+	if(k==0) return f-(hfsq-s*(hfsq+R)); else
+		 return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
+}
diff --git a/modules/fdlibm/src/s_nearbyint.cpp b/modules/fdlibm/src/s_nearbyint.cpp
new file mode 100644
index 000000000..532bb5d8d
--- /dev/null
+++ b/modules/fdlibm/src/s_nearbyint.cpp
@@ -0,0 +1,58 @@
+/*-
+ * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+#include <fenv.h>
+#include "math_private.h"
+
+/*
+ * We save and restore the floating-point environment to avoid raising
+ * an inexact exception.  We can get away with using fesetenv()
+ * instead of feclearexcept()/feupdateenv() to restore the environment
+ * because the only exception defined for rint() is overflow, and
+ * rounding can't overflow as long as emax >= p.
+ *
+ * The volatile keyword is needed below because clang incorrectly assumes
+ * that rint won't raise any floating-point exceptions. Declaring ret volatile
+ * is sufficient to trick the compiler into doing the right thing.
+ */
+#define	DECL(type, fn, rint)	\
+type				\
+fn(type x)			\
+{				\
+	volatile type ret;	\
+	fenv_t env;		\
+				\
+	fegetenv(&env);		\
+	ret = rint(x);		\
+	fesetenv(&env);		\
+	return (ret);		\
+}
+
+DECL(double, nearbyint, rint)
+DECL(float, nearbyintf, rintf)
diff --git a/modules/fdlibm/src/s_rint.cpp b/modules/fdlibm/src/s_rint.cpp
new file mode 100644
index 000000000..19171f87f
--- /dev/null
+++ b/modules/fdlibm/src/s_rint.cpp
@@ -0,0 +1,87 @@
+/* @(#)s_rint.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * rint(x)
+ * Return x rounded to integral value according to the prevailing
+ * rounding mode.
+ * Method:
+ *	Using floating addition.
+ * Exception:
+ *	Inexact flag raised if x not equal to rint(x).
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+TWO52[2]={
+  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
+ -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
+};
+
+double
+rint(double x)
+{
+	int32_t i0,j0,sx;
+	u_int32_t i,i1;
+	double w,t;
+	EXTRACT_WORDS(i0,i1,x);
+	sx = (i0>>31)&1;
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) {
+		if(((i0&0x7fffffff)|i1)==0) return x;
+		i1 |= (i0&0x0fffff);
+		i0 &= 0xfffe0000;
+		i0 |= ((i1|-i1)>>12)&0x80000;
+		SET_HIGH_WORD(x,i0);
+	        STRICT_ASSIGN(double,w,TWO52[sx]+x);
+	        t =  w-TWO52[sx];
+		GET_HIGH_WORD(i0,t);
+		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
+	        return t;
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		i>>=1;
+		if(((i0&i)|i1)!=0) {
+		    /*
+		     * Some bit is set after the 0.5 bit.  To avoid the
+		     * possibility of errors from double rounding in
+		     * w = TWO52[sx]+x, adjust the 0.25 bit to a lower
+		     * guard bit.  We do this for all j0<=51.  The
+		     * adjustment is trickiest for j0==18 and j0==19
+		     * since then it spans the word boundary.
+		     */
+		    if(j0==19) i1 = 0x40000000; else
+		    if(j0==18) i1 = 0x80000000; else
+		    i0 = (i0&(~i))|((0x20000)>>j0);
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    i>>=1;
+	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
+	}
+	INSERT_WORDS(x,i0,i1);
+	STRICT_ASSIGN(double,w,TWO52[sx]+x);
+	return w-TWO52[sx];
+}
diff --git a/modules/fdlibm/src/s_rintf.cpp b/modules/fdlibm/src/s_rintf.cpp
new file mode 100644
index 000000000..3a729b005
--- /dev/null
+++ b/modules/fdlibm/src/s_rintf.cpp
@@ -0,0 +1,52 @@
+/* s_rintf.c -- float version of s_rint.c.
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+#include <float.h>
+#include <stdint.h>
+
+#include "math_private.h"
+
+static const float
+TWO23[2]={
+  8.3886080000e+06, /* 0x4b000000 */
+ -8.3886080000e+06, /* 0xcb000000 */
+};
+
+float
+rintf(float x)
+{
+	int32_t i0,j0,sx;
+	float w,t;
+	GET_FLOAT_WORD(i0,x);
+	sx = (i0>>31)&1;
+	j0 = ((i0>>23)&0xff)-0x7f;
+	if(j0<23) {
+	    if(j0<0) {
+		if((i0&0x7fffffff)==0) return x;
+		STRICT_ASSIGN(float,w,TWO23[sx]+x);
+	        t =  w-TWO23[sx];
+		GET_FLOAT_WORD(i0,t);
+		SET_FLOAT_WORD(t,(i0&0x7fffffff)|(sx<<31));
+	        return t;
+	    }
+	    STRICT_ASSIGN(float,w,TWO23[sx]+x);
+	    return w-TWO23[sx];
+	}
+	if(j0==0x80) return x+x;	/* inf or NaN */
+	else return x;			/* x is integral */
+}
diff --git a/modules/fdlibm/src/s_scalbn.cpp b/modules/fdlibm/src/s_scalbn.cpp
new file mode 100644
index 000000000..5dbf58c23
--- /dev/null
+++ b/modules/fdlibm/src/s_scalbn.cpp
@@ -0,0 +1,60 @@
+/* @(#)s_scalbn.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#ifndef lint
+//static char rcsid[] = "$FreeBSD$";
+#endif
+
+/*
+ * scalbn (double x, int n)
+ * scalbn(x,n) returns x* 2**n  computed by  exponent
+ * manipulation rather than by actually performing an
+ * exponentiation or a multiplication.
+ */
+
+//#include <sys/cdefs.h>
+#include <float.h>
+
+#include "math_private.h"
+
+static const double
+two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
+twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
+huge   = 1.0e+300,
+tiny   = 1.0e-300;
+
+double
+scalbn (double x, int n)
+{
+	int32_t k,hx,lx;
+	EXTRACT_WORDS(hx,lx,x);
+        k = (hx&0x7ff00000)>>20;		/* extract exponent */
+        if (k==0) {				/* 0 or subnormal x */
+            if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
+	    x *= two54;
+	    GET_HIGH_WORD(hx,x);
+	    k = ((hx&0x7ff00000)>>20) - 54;
+            if (n< -50000) return tiny*x; 	/*underflow*/
+	    }
+        if (k==0x7ff) return x+x;		/* NaN or Inf */
+        k = k+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 (n > 50000) 	/* in case integer overflow in n+k */
+		return huge*copysign(huge,x);	/*overflow*/
+	    else return tiny*copysign(tiny,x); 	/*underflow*/
+        k += 54;				/* subnormal result */
+	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
+        return x*twom54;
+}
diff --git a/modules/fdlibm/src/s_tanh.cpp b/modules/fdlibm/src/s_tanh.cpp
new file mode 100644
index 000000000..238973fce
--- /dev/null
+++ b/modules/fdlibm/src/s_tanh.cpp
@@ -0,0 +1,79 @@
+/* @(#)s_tanh.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/* Tanh(x)
+ * Return the Hyperbolic Tangent of x
+ *
+ * Method :
+ *				       x    -x
+ *				      e  - e
+ *	0. tanh(x) is defined to be -----------
+ *				       x    -x
+ *				      e  + e
+ *	1. reduce x to non-negative by tanh(-x) = -tanh(x).
+ *	2.  0      <= x <  2**-28 : tanh(x) := x with inexact if x != 0
+ *					        -t
+ *	    2**-28 <= x <  1      : tanh(x) := -----; t = expm1(-2x)
+ *					       t + 2
+ *						     2
+ *	    1      <= x <  22     : tanh(x) := 1 - -----; t = expm1(2x)
+ *						   t + 2
+ *	    22     <= x <= INF    : tanh(x) := 1.
+ *
+ * Special cases:
+ *	tanh(NaN) is NaN;
+ *	only tanh(0)=0 is exact for finite argument.
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const volatile double tiny = 1.0e-300;
+static const double one = 1.0, two = 2.0, huge = 1.0e300;
+
+double
+tanh(double x)
+{
+	double t,z;
+	int32_t jx,ix;
+
+	GET_HIGH_WORD(jx,x);
+	ix = jx&0x7fffffff;
+
+    /* x is INF or NaN */
+	if(ix>=0x7ff00000) {
+	    if (jx>=0) return one/x+one;    /* tanh(+-inf)=+-1 */
+	    else       return one/x-one;    /* tanh(NaN) = NaN */
+	}
+
+    /* |x| < 22 */
+	if (ix < 0x40360000) {		/* |x|<22 */
+	    if (ix<0x3e300000) {	/* |x|<2**-28 */
+		if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */
+	    }
+	    if (ix>=0x3ff00000) {	/* |x|>=1  */
+		t = expm1(two*fabs(x));
+		z = one - two/(t+two);
+	    } else {
+	        t = expm1(-two*fabs(x));
+	        z= -t/(t+two);
+	    }
+    /* |x| >= 22, return +-1 */
+	} else {
+	    z = one - tiny;		/* raise inexact flag */
+	}
+	return (jx>=0)? z: -z;
+}
diff --git a/modules/fdlibm/src/s_trunc.cpp b/modules/fdlibm/src/s_trunc.cpp
new file mode 100644
index 000000000..d2294a272
--- /dev/null
+++ b/modules/fdlibm/src/s_trunc.cpp
@@ -0,0 +1,62 @@
+/* @(#)s_floor.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * trunc(x)
+ * Return x rounded toward 0 to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to trunc(x).
+ */
+
+#include <float.h>
+
+#include "math_private.h"
+
+static const double huge = 1.0e300;
+
+double
+trunc(double x)
+{
+	int32_t i0,i1,j0;
+	u_int32_t i;
+	EXTRACT_WORDS(i0,i1,x);
+	j0 = ((i0>>20)&0x7ff)-0x3ff;
+	if(j0<20) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0) {/* |x|<1, so return 0*sign(x) */
+		    i0 &= 0x80000000U;
+		    i1 = 0;
+		}
+	    } else {
+		i = (0x000fffff)>>j0;
+		if(((i0&i)|i1)==0) return x; /* x is integral */
+		if(huge+x>0.0) {	/* raise inexact flag */
+		    i0 &= (~i); i1=0;
+		}
+	    }
+	} else if (j0>51) {
+	    if(j0==0x400) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	} else {
+	    i = ((u_int32_t)(0xffffffff))>>(j0-20);
+	    if((i1&i)==0) return x;	/* x is integral */
+	    if(huge+x>0.0)		/* raise inexact flag */
+		i1 &= (~i);
+	}
+	INSERT_WORDS(x,i0,i1);
+	return x;
+}
diff --git a/modules/fdlibm/src/s_truncf.cpp b/modules/fdlibm/src/s_truncf.cpp
new file mode 100644
index 000000000..4853a4450
--- /dev/null
+++ b/modules/fdlibm/src/s_truncf.cpp
@@ -0,0 +1,52 @@
+/* @(#)s_floor.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+//#include <sys/cdefs.h>
+//__FBSDID("$FreeBSD$");
+
+/*
+ * truncf(x)
+ * Return x rounded toward 0 to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to truncf(x).
+ */
+
+#include "math_private.h"
+
+static const float huge = 1.0e30F;
+
+float
+truncf(float x)
+{
+	int32_t i0,j0;
+	u_int32_t i;
+	GET_FLOAT_WORD(i0,x);
+	j0 = ((i0>>23)&0xff)-0x7f;
+	if(j0<23) {
+	    if(j0<0) { 	/* raise inexact if x != 0 */
+		if(huge+x>0.0F)		/* |x|<1, so return 0*sign(x) */
+		    i0 &= 0x80000000;
+	    } else {
+		i = (0x007fffff)>>j0;
+		if((i0&i)==0) return x; /* x is integral */
+		if(huge+x>0.0F)		/* raise inexact flag */
+		    i0 &= (~i);
+	    }
+	} else {
+	    if(j0==0x80) return x+x;	/* inf or NaN */
+	    else return x;		/* x is integral */
+	}
+	SET_FLOAT_WORD(x,i0);
+	return x;
+}