Backed out changeset 2f7c7d7ddacc (bug 933257)
authorTooru Fujisawa <arai_a@mac.com>
Sun, 13 Mar 2016 04:55:00 +0900
changeset 288469 7009000732a2ceffece3786738ecf2f46abbcea7
parent 288468 ff0cb614dbdbd6ac705b728795b4f6539b1ffde0
child 288470 5cbe6e726f6c725fe205f7bce29c4a76e47bfc1a
push id30079
push userryanvm@gmail.com
push dateSat, 12 Mar 2016 20:24:19 +0000
treeherdermozilla-central@d1d47ba19ce9 [default view] [failures only]
perfherder[talos] [build metrics] [platform microbench] (compared to previous push)
bugs933257
milestone48.0a1
backs out2f7c7d7ddacc51bcd6ad6a7d78fc99b80ef75724
first release with
nightly linux32
nightly linux64
nightly mac
nightly win32
nightly win64
last release without
nightly linux32
nightly linux64
nightly mac
nightly win32
nightly win64
Backed out changeset 2f7c7d7ddacc (bug 933257)
modules/fdlibm/README.mozilla
modules/fdlibm/src/e_acos.cpp
modules/fdlibm/src/e_acosh.cpp
modules/fdlibm/src/e_asin.cpp
modules/fdlibm/src/e_atan2.cpp
modules/fdlibm/src/e_atanh.cpp
modules/fdlibm/src/e_cosh.cpp
modules/fdlibm/src/e_exp.cpp
modules/fdlibm/src/e_hypot.cpp
modules/fdlibm/src/e_log.cpp
modules/fdlibm/src/e_log10.cpp
modules/fdlibm/src/e_log2.cpp
modules/fdlibm/src/e_pow.cpp
modules/fdlibm/src/e_rem_pio2.cpp
modules/fdlibm/src/e_sinh.cpp
modules/fdlibm/src/e_sqrt.cpp
modules/fdlibm/src/fdlibm.h
modules/fdlibm/src/k_cos.cpp
modules/fdlibm/src/k_exp.cpp
modules/fdlibm/src/k_log.h
modules/fdlibm/src/k_rem_pio2.cpp
modules/fdlibm/src/k_sin.cpp
modules/fdlibm/src/k_tan.cpp
modules/fdlibm/src/math_private.h
modules/fdlibm/src/s_asinh.cpp
modules/fdlibm/src/s_atan.cpp
modules/fdlibm/src/s_cbrt.cpp
modules/fdlibm/src/s_ceil.cpp
modules/fdlibm/src/s_ceilf.cpp
modules/fdlibm/src/s_copysign.cpp
modules/fdlibm/src/s_cos.cpp
modules/fdlibm/src/s_expm1.cpp
modules/fdlibm/src/s_fabs.cpp
modules/fdlibm/src/s_floor.cpp
modules/fdlibm/src/s_floorf.cpp
modules/fdlibm/src/s_log1p.cpp
modules/fdlibm/src/s_scalbn.cpp
modules/fdlibm/src/s_sin.cpp
modules/fdlibm/src/s_tan.cpp
modules/fdlibm/src/s_tanh.cpp
modules/fdlibm/src/s_trunc.cpp
--- a/modules/fdlibm/README.mozilla
+++ b/modules/fdlibm/README.mozilla
@@ -6,12 +6,12 @@ Upstream code can be viewed at
 
 Each file is downloaded separately, as cloning whole repository takes so much
 resources.
 
 The in-tree copy is updated by running
   sh update.sh
 from within the modules/fdlibm directory.
 
-Current version: [commit bcea9d50b15e4f0027a5dd526e0e2a612238471e].
+Current version: [commit 0000000000000000000000000000000000000000].
 
 patches 01-16 fixes files to be usable within mozilla-central tree.
 See https://bugzilla.mozilla.org/show_bug.cgi?id=933257
deleted file mode 100644
--- a/modules/fdlibm/src/e_acos.cpp
+++ /dev/null
@@ -1,106 +0,0 @@
-
-/* @(#)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);
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_acosh.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-
-/* @(#)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));
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_asin.cpp
+++ /dev/null
@@ -1,112 +0,0 @@
-
-/* @(#)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;    
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_atan2.cpp
+++ /dev/null
@@ -1,124 +0,0 @@
-
-/* @(#)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(-,-) */
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_atanh.cpp
+++ /dev/null
@@ -1,63 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_cosh.cpp
+++ /dev/null
@@ -1,80 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_exp.cpp
+++ /dev/null
@@ -1,162 +0,0 @@
-
-/* @(#)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;
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_hypot.cpp
+++ /dev/null
@@ -1,126 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_log.cpp
+++ /dev/null
@@ -1,142 +0,0 @@
-
-/* @(#)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);
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_log10.cpp
+++ /dev/null
@@ -1,89 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_log2.cpp
+++ /dev/null
@@ -1,112 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_pow.cpp
+++ /dev/null
@@ -1,305 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_rem_pio2.cpp
+++ /dev/null
@@ -1,187 +0,0 @@
-
-/* @(#)e_rem_pio2.c 1.4 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.
- * ====================================================
- *
- * Optimized by Bruce D. Evans.
- */
-
-//#include <sys/cdefs.h>
-//__FBSDID("$FreeBSD$");
-
-/* __ieee754_rem_pio2(x,y)
- * 
- * return the remainder of x rem pi/2 in y[0]+y[1] 
- * use __kernel_rem_pio2()
- */
-
-#include <float.h>
-
-#include "mozilla/Attributes.h"
-#include "math_private.h"
-
-/*
- * invpio2:  53 bits of 2/pi
- * pio2_1:   first  33 bit of pi/2
- * pio2_1t:  pi/2 - pio2_1
- * pio2_2:   second 33 bit of pi/2
- * pio2_2t:  pi/2 - (pio2_1+pio2_2)
- * pio2_3:   third  33 bit of pi/2
- * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
- */
-
-static const double
-zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
-two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
-invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
-pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
-pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
-pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
-pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
-pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
-pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
-
-#ifdef INLINE_REM_PIO2
-static MOZ_ALWAYS_INLINE
-#endif
-int
-__ieee754_rem_pio2(double x, double *y)
-{
-	double z,w,t,r,fn;
-	double tx[3],ty[2];
-	int32_t e0,i,j,nx,n,ix,hx;
-	u_int32_t low;
-
-	GET_HIGH_WORD(hx,x);		/* high word of x */
-	ix = hx&0x7fffffff;
-#if 0 /* Must be handled in caller. */
-	if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
-	    {y[0] = x; y[1] = 0; return 0;}
-#endif
-	if (ix <= 0x400f6a7a) {		/* |x| ~<= 5pi/4 */
-	    if ((ix & 0xfffff) == 0x921fb)  /* |x| ~= pi/2 or 2pi/2 */
-		goto medium;		/* cancellation -- use medium case */
-	    if (ix <= 0x4002d97c) {	/* |x| ~<= 3pi/4 */
-		if (hx > 0) {
-		    z = x - pio2_1;	/* one round good to 85 bits */
-		    y[0] = z - pio2_1t;
-		    y[1] = (z-y[0])-pio2_1t;
-		    return 1;
-		} else {
-		    z = x + pio2_1;
-		    y[0] = z + pio2_1t;
-		    y[1] = (z-y[0])+pio2_1t;
-		    return -1;
-		}
-	    } else {
-		if (hx > 0) {
-		    z = x - 2*pio2_1;
-		    y[0] = z - 2*pio2_1t;
-		    y[1] = (z-y[0])-2*pio2_1t;
-		    return 2;
-		} else {
-		    z = x + 2*pio2_1;
-		    y[0] = z + 2*pio2_1t;
-		    y[1] = (z-y[0])+2*pio2_1t;
-		    return -2;
-		}
-	    }
-	}
-	if (ix <= 0x401c463b) {		/* |x| ~<= 9pi/4 */
-	    if (ix <= 0x4015fdbc) {	/* |x| ~<= 7pi/4 */
-		if (ix == 0x4012d97c)	/* |x| ~= 3pi/2 */
-		    goto medium;
-		if (hx > 0) {
-		    z = x - 3*pio2_1;
-		    y[0] = z - 3*pio2_1t;
-		    y[1] = (z-y[0])-3*pio2_1t;
-		    return 3;
-		} else {
-		    z = x + 3*pio2_1;
-		    y[0] = z + 3*pio2_1t;
-		    y[1] = (z-y[0])+3*pio2_1t;
-		    return -3;
-		}
-	    } else {
-		if (ix == 0x401921fb)	/* |x| ~= 4pi/2 */
-		    goto medium;
-		if (hx > 0) {
-		    z = x - 4*pio2_1;
-		    y[0] = z - 4*pio2_1t;
-		    y[1] = (z-y[0])-4*pio2_1t;
-		    return 4;
-		} else {
-		    z = x + 4*pio2_1;
-		    y[0] = z + 4*pio2_1t;
-		    y[1] = (z-y[0])+4*pio2_1t;
-		    return -4;
-		}
-	    }
-	}
-	if(ix<0x413921fb) {	/* |x| ~< 2^20*(pi/2), medium size */
-medium:
-	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-	    double const_0x1_8p52 = pow(2, 52) + pow(2, 51);
-	    STRICT_ASSIGN(double,fn,x*invpio2+const_0x1_8p52);
-	    fn = fn-const_0x1_8p52;
-#ifdef HAVE_EFFICIENT_IRINT
-	    n  = irint(fn);
-#else
-	    n  = (int32_t)fn;
-#endif
-	    r  = x-fn*pio2_1;
-	    w  = fn*pio2_1t;	/* 1st round good to 85 bit */
-	    {
-	        u_int32_t high;
-	        j  = ix>>20;
-	        y[0] = r-w; 
-		GET_HIGH_WORD(high,y[0]);
-	        i = j-((high>>20)&0x7ff);
-	        if(i>16) {  /* 2nd iteration needed, good to 118 */
-		    t  = r;
-		    w  = fn*pio2_2;	
-		    r  = t-w;
-		    w  = fn*pio2_2t-((t-r)-w);	
-		    y[0] = r-w;
-		    GET_HIGH_WORD(high,y[0]);
-		    i = j-((high>>20)&0x7ff);
-		    if(i>49)  {	/* 3rd iteration need, 151 bits acc */
-		    	t  = r;	/* will cover all possible cases */
-		    	w  = fn*pio2_3;	
-		    	r  = t-w;
-		    	w  = fn*pio2_3t-((t-r)-w);	
-		    	y[0] = r-w;
-		    }
-		}
-	    }
-	    y[1] = (r-y[0])-w;
-	    return n;
-	}
-    /* 
-     * all other (large) arguments
-     */
-	if(ix>=0x7ff00000) {		/* x is inf or NaN */
-	    y[0]=y[1]=x-x; return 0;
-	}
-    /* set z = scalbn(|x|,ilogb(x)-23) */
-	GET_LOW_WORD(low,x);
-	e0 	= (ix>>20)-1046;	/* e0 = ilogb(z)-23; */
-	INSERT_WORDS(z, ix - ((int32_t)(e0<<20)), low);
-	for(i=0;i<2;i++) {
-		tx[i] = (double)((int32_t)(z));
-		z     = (z-tx[i])*two24;
-	}
-	tx[2] = z;
-	nx = 3;
-	while(tx[nx-1]==zero) nx--;	/* skip zero term */
-	n  =  __kernel_rem_pio2(tx,ty,e0,nx,1);
-	if(hx<0) {y[0] = -ty[0]; y[1] = -ty[1]; return -n;}
-	y[0] = ty[0]; y[1] = ty[1]; return n;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_sinh.cpp
+++ /dev/null
@@ -1,74 +0,0 @@
-
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/e_sqrt.cpp
+++ /dev/null
@@ -1,446 +0,0 @@
-
-/* @(#)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).	
- 
- */
- 
deleted file mode 100644
--- a/modules/fdlibm/src/fdlibm.h
+++ /dev/null
@@ -1,63 +0,0 @@
-/*
- * ====================================================
- * 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	cos(double);
-double	sin(double);
-double	tan(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	ceil(double);
-float	ceilf(float);
-double	fabs(double);
-double	floor(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	copysign(double, double);
-double	scalbn(double, int);
-double	trunc(double);
-
-float	floorf(float);
-
-} /* namespace fdlibm */
-
-#endif /* mozilla_imported_fdlibm_h */
deleted file mode 100644
--- a/modules/fdlibm/src/k_cos.cpp
+++ /dev/null
@@ -1,78 +0,0 @@
-
-/* @(#)k_cos.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$");
-
-/*
- * __kernel_cos( x,  y )
- * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
- * Input x is assumed to be bounded by ~pi/4 in magnitude.
- * Input y is the tail of x. 
- *
- * Algorithm
- *	1. Since cos(-x) = cos(x), we need only to consider positive x.
- *	2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
- *	3. cos(x) is approximated by a polynomial of degree 14 on
- *	   [0,pi/4]
- *		  	                 4            14
- *	   	cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
- *	   where the remez error is
- *	
- * 	|              2     4     6     8     10    12     14 |     -58
- * 	|cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
- * 	|    					               | 
- * 
- * 	               4     6     8     10    12     14 
- *	4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
- *	       cos(x) ~ 1 - x*x/2 + r
- *	   since cos(x+y) ~ cos(x) - sin(x)*y 
- *			  ~ cos(x) - x*y,
- *	   a correction term is necessary in cos(x) and hence
- *		cos(x+y) = 1 - (x*x/2 - (r - x*y))
- *	   For better accuracy, rearrange to
- *		cos(x+y) ~ w + (tmp + (r-x*y))
- *	   where w = 1 - x*x/2 and tmp is a tiny correction term
- *	   (1 - x*x/2 == w + tmp exactly in infinite precision).
- *	   The exactness of w + tmp in infinite precision depends on w
- *	   and tmp having the same precision as x.  If they have extra
- *	   precision due to compiler bugs, then the extra precision is
- *	   only good provided it is retained in all terms of the final
- *	   expression for cos().  Retention happens in all cases tested
- *	   under FreeBSD, so don't pessimize things by forcibly clipping
- *	   any extra precision in w.
- */
-
-#include "math_private.h"
-
-static const double
-one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-C1  =  4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
-C2  = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
-C3  =  2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
-C4  = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
-C5  =  2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
-C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
-
-double
-__kernel_cos(double x, double y)
-{
-	double hz,z,r,w;
-
-	z  = x*x;
-	w  = z*z;
-	r  = z*(C1+z*(C2+z*C3)) + w*w*(C4+z*(C5+z*C6));
-	hz = 0.5*z;
-	w  = one-hz;
-	return w + (((one-w)-hz) + (z*r-x*y));
-}
deleted file mode 100644
--- a/modules/fdlibm/src/k_exp.cpp
+++ /dev/null
@@ -1,81 +0,0 @@
-/*-
- * 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);
-}
deleted file mode 100644
--- a/modules/fdlibm/src/k_log.h
+++ /dev/null
@@ -1,100 +0,0 @@
-
-/* @(#)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);
-}
deleted file mode 100644
--- a/modules/fdlibm/src/k_rem_pio2.cpp
+++ /dev/null
@@ -1,442 +0,0 @@
-
-/* @(#)k_rem_pio2.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$");
-
-/*
- * __kernel_rem_pio2(x,y,e0,nx,prec)
- * double x[],y[]; int e0,nx,prec;
- * 
- * __kernel_rem_pio2 return the last three digits of N with 
- *		y = x - N*pi/2
- * so that |y| < pi/2.
- *
- * The method is to compute the integer (mod 8) and fraction parts of 
- * (2/pi)*x without doing the full multiplication. In general we
- * skip the part of the product that are known to be a huge integer (
- * more accurately, = 0 mod 8 ). Thus the number of operations are
- * independent of the exponent of the input.
- *
- * (2/pi) is represented by an array of 24-bit integers in ipio2[].
- *
- * Input parameters:
- * 	x[]	The input value (must be positive) is broken into nx 
- *		pieces of 24-bit integers in double precision format.
- *		x[i] will be the i-th 24 bit of x. The scaled exponent 
- *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 
- *		match x's up to 24 bits.
- *
- *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
- *			e0 = ilogb(z)-23
- *			z  = scalbn(z,-e0)
- *		for i = 0,1,2
- *			x[i] = floor(z)
- *			z    = (z-x[i])*2**24
- *
- *
- *	y[]	output result in an array of double precision numbers.
- *		The dimension of y[] is:
- *			24-bit  precision	1
- *			53-bit  precision	2
- *			64-bit  precision	2
- *			113-bit precision	3
- *		The actual value is the sum of them. Thus for 113-bit
- *		precison, one may have to do something like:
- *
- *		long double t,w,r_head, r_tail;
- *		t = (long double)y[2] + (long double)y[1];
- *		w = (long double)y[0];
- *		r_head = t+w;
- *		r_tail = w - (r_head - t);
- *
- *	e0	The exponent of x[0]. Must be <= 16360 or you need to
- *              expand the ipio2 table.
- *
- *	nx	dimension of x[]
- *
- *  	prec	an integer indicating the precision:
- *			0	24  bits (single)
- *			1	53  bits (double)
- *			2	64  bits (extended)
- *			3	113 bits (quad)
- *
- * External function:
- *	double scalbn(), floor();
- *
- *
- * Here is the description of some local variables:
- *
- * 	jk	jk+1 is the initial number of terms of ipio2[] needed
- *		in the computation. The minimum and recommended value
- *		for jk is 3,4,4,6 for single, double, extended, and quad.
- *		jk+1 must be 2 larger than you might expect so that our
- *		recomputation test works. (Up to 24 bits in the integer
- *		part (the 24 bits of it that we compute) and 23 bits in
- *		the fraction part may be lost to cancelation before we
- *		recompute.)
- *
- * 	jz	local integer variable indicating the number of 
- *		terms of ipio2[] used. 
- *
- *	jx	nx - 1
- *
- *	jv	index for pointing to the suitable ipio2[] for the
- *		computation. In general, we want
- *			( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
- *		is an integer. Thus
- *			e0-3-24*jv >= 0 or (e0-3)/24 >= jv
- *		Hence jv = max(0,(e0-3)/24).
- *
- *	jp	jp+1 is the number of terms in PIo2[] needed, jp = jk.
- *
- * 	q[]	double array with integral value, representing the
- *		24-bits chunk of the product of x and 2/pi.
- *
- *	q0	the corresponding exponent of q[0]. Note that the
- *		exponent for q[i] would be q0-24*i.
- *
- *	PIo2[]	double precision array, obtained by cutting pi/2
- *		into 24 bits chunks. 
- *
- *	f[]	ipio2[] in floating point 
- *
- *	iq[]	integer array by breaking up q[] in 24-bits chunk.
- *
- *	fq[]	final product of x*(2/pi) in fq[0],..,fq[jk]
- *
- *	ih	integer. If >0 it indicates q[] is >= 0.5, hence
- *		it also indicates the *sign* of the result.
- *
- */
-
-
-/*
- * 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 int init_jk[] = {3,4,4,6}; /* initial value for jk */
-
-/*
- * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
- *
- *		integer array, contains the (24*i)-th to (24*i+23)-th 
- *		bit of 2/pi after binary point. The corresponding 
- *		floating value is
- *
- *			ipio2[i] * 2^(-24(i+1)).
- *
- * NB: This table must have at least (e0-3)/24 + jk terms.
- *     For quad precision (e0 <= 16360, jk = 6), this is 686.
- */
-static const int32_t ipio2[] = {
-0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
-0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
-0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
-0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
-0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
-0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
-0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
-0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
-0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
-0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
-0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
-
-#if LDBL_MAX_EXP > 1024
-#if LDBL_MAX_EXP > 16384
-#error "ipio2 table needs to be expanded"
-#endif
-0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
-0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
-0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
-0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
-0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
-0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
-0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
-0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
-0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
-0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
-0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
-0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
-0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
-0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
-0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
-0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
-0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
-0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
-0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
-0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
-0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
-0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
-0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
-0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
-0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
-0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
-0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
-0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
-0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
-0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
-0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
-0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
-0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
-0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
-0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
-0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
-0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
-0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
-0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
-0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
-0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
-0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
-0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
-0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
-0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
-0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
-0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
-0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
-0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
-0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
-0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
-0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
-0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
-0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
-0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
-0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
-0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
-0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
-0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
-0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
-0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
-0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
-0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
-0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
-0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
-0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
-0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
-0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
-0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
-0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
-0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
-0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
-0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
-0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
-0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
-0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
-0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
-0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
-0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
-0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
-0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
-0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
-0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
-0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
-0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
-0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
-0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
-0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
-0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
-0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
-0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
-0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
-0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
-0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
-0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
-0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
-0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
-0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
-0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
-0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
-0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
-0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
-0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
-0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
-#endif
-
-};
-
-static const double PIo2[] = {
-  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
-  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
-  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
-  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
-  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
-  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
-  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
-  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
-};
-
-static const double			
-zero   = 0.0,
-one    = 1.0,
-two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
-twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
-
-int
-__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
-{
-	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
-	double z,fw,f[20],fq[20],q[20];
-
-    /* initialize jk*/
-	jk = init_jk[prec];
-	jp = jk;
-
-    /* determine jx,jv,q0, note that 3>q0 */
-	jx =  nx-1;
-	jv = (e0-3)/24; if(jv<0) jv=0;
-	q0 =  e0-24*(jv+1);
-
-    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
-	j = jv-jx; m = jx+jk;
-	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
-
-    /* compute q[0],q[1],...q[jk] */
-	for (i=0;i<=jk;i++) {
-	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
-	}
-
-	jz = jk;
-recompute:
-    /* distill q[] into iq[] reversingly */
-	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
-	    fw    =  (double)((int32_t)(twon24* z));
-	    iq[i] =  (int32_t)(z-two24*fw);
-	    z     =  q[j-1]+fw;
-	}
-
-    /* compute n */
-	z  = scalbn(z,q0);		/* actual value of z */
-	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
-	n  = (int32_t) z;
-	z -= (double)n;
-	ih = 0;
-	if(q0>0) {	/* need iq[jz-1] to determine n */
-	    i  = (iq[jz-1]>>(24-q0)); n += i;
-	    iq[jz-1] -= i<<(24-q0);
-	    ih = iq[jz-1]>>(23-q0);
-	} 
-	else if(q0==0) ih = iq[jz-1]>>23;
-	else if(z>=0.5) ih=2;
-
-	if(ih>0) {	/* q > 0.5 */
-	    n += 1; carry = 0;
-	    for(i=0;i<jz ;i++) {	/* compute 1-q */
-		j = iq[i];
-		if(carry==0) {
-		    if(j!=0) {
-			carry = 1; iq[i] = 0x1000000- j;
-		    }
-		} else  iq[i] = 0xffffff - j;
-	    }
-	    if(q0>0) {		/* rare case: chance is 1 in 12 */
-	        switch(q0) {
-	        case 1:
-	    	   iq[jz-1] &= 0x7fffff; break;
-	    	case 2:
-	    	   iq[jz-1] &= 0x3fffff; break;
-	        }
-	    }
-	    if(ih==2) {
-		z = one - z;
-		if(carry!=0) z -= scalbn(one,q0);
-	    }
-	}
-
-    /* check if recomputation is needed */
-	if(z==zero) {
-	    j = 0;
-	    for (i=jz-1;i>=jk;i--) j |= iq[i];
-	    if(j==0) { /* need recomputation */
-		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
-
-		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
-		    f[jx+i] = (double) ipio2[jv+i];
-		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
-		    q[i] = fw;
-		}
-		jz += k;
-		goto recompute;
-	    }
-	}
-
-    /* chop off zero terms */
-	if(z==0.0) {
-	    jz -= 1; q0 -= 24;
-	    while(iq[jz]==0) { jz--; q0-=24;}
-	} else { /* break z into 24-bit if necessary */
-	    z = scalbn(z,-q0);
-	    if(z>=two24) { 
-		fw = (double)((int32_t)(twon24*z));
-		iq[jz] = (int32_t)(z-two24*fw);
-		jz += 1; q0 += 24;
-		iq[jz] = (int32_t) fw;
-	    } else iq[jz] = (int32_t) z ;
-	}
-
-    /* convert integer "bit" chunk to floating-point value */
-	fw = scalbn(one,q0);
-	for(i=jz;i>=0;i--) {
-	    q[i] = fw*(double)iq[i]; fw*=twon24;
-	}
-
-    /* compute PIo2[0,...,jp]*q[jz,...,0] */
-	for(i=jz;i>=0;i--) {
-	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
-	    fq[jz-i] = fw;
-	}
-
-    /* compress fq[] into y[] */
-	switch(prec) {
-	    case 0:
-		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i];
-		y[0] = (ih==0)? fw: -fw; 
-		break;
-	    case 1:
-	    case 2:
-		fw = 0.0;
-		for (i=jz;i>=0;i--) fw += fq[i]; 
-		STRICT_ASSIGN(double,fw,fw);
-		y[0] = (ih==0)? fw: -fw; 
-		fw = fq[0]-fw;
-		for (i=1;i<=jz;i++) fw += fq[i];
-		y[1] = (ih==0)? fw: -fw; 
-		break;
-	    case 3:	/* painful */
-		for (i=jz;i>0;i--) {
-		    fw      = fq[i-1]+fq[i]; 
-		    fq[i]  += fq[i-1]-fw;
-		    fq[i-1] = fw;
-		}
-		for (i=jz;i>1;i--) {
-		    fw      = fq[i-1]+fq[i]; 
-		    fq[i]  += fq[i-1]-fw;
-		    fq[i-1] = fw;
-		}
-		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 
-		if(ih==0) {
-		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
-		} else {
-		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
-		}
-	}
-	return n&7;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/k_sin.cpp
+++ /dev/null
@@ -1,69 +0,0 @@
-
-/* @(#)k_sin.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$");
-
-/* __kernel_sin( x, y, iy)
- * kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
- * Input x is assumed to be bounded by ~pi/4 in magnitude.
- * Input y is the tail of x.
- * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). 
- *
- * Algorithm
- *	1. Since sin(-x) = -sin(x), we need only to consider positive x. 
- *	2. Callers must return sin(-0) = -0 without calling here since our
- *	   odd polynomial is not evaluated in a way that preserves -0.
- *	   Callers may do the optimization sin(x) ~ x for tiny x.
- *	3. sin(x) is approximated by a polynomial of degree 13 on
- *	   [0,pi/4]
- *		  	         3            13
- *	   	sin(x) ~ x + S1*x + ... + S6*x
- *	   where
- *	
- * 	|sin(x)         2     4     6     8     10     12  |     -58
- * 	|----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
- * 	|  x 					           | 
- * 
- *	4. sin(x+y) = sin(x) + sin'(x')*y
- *		    ~ sin(x) + (1-x*x/2)*y
- *	   For better accuracy, let 
- *		     3      2      2      2      2
- *		r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
- *	   then                   3    2
- *		sin(x) = x + (S1*x + (x *(r-y/2)+y))
- */
-
-#include "math_private.h"
-
-static const double
-half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-S1  = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
-S2  =  8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
-S3  = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
-S4  =  2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
-S5  = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
-S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
-
-double
-__kernel_sin(double x, double y, int iy)
-{
-	double z,r,v,w;
-
-	z	=  x*x;
-	w	=  z*z;
-	r	=  S2+z*(S3+z*S4) + z*w*(S5+z*S6);
-	v	=  z*x;
-	if(iy==0) return x+v*(S1+z*r);
-	else      return x-((z*(half*y-v*r)-y)-v*S1);
-}
deleted file mode 100644
--- a/modules/fdlibm/src/k_tan.cpp
+++ /dev/null
@@ -1,131 +0,0 @@
-/* @(#)k_tan.c 1.5 04/04/22 SMI */
-
-/*
- * ====================================================
- * Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
- *
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
-
-/* INDENT OFF */
-//#include <sys/cdefs.h>
-//__FBSDID("$FreeBSD$");
-
-/* __kernel_tan( x, y, k )
- * kernel tan function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854
- * Input x is assumed to be bounded by ~pi/4 in magnitude.
- * Input y is the tail of x.
- * Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.
- *
- * Algorithm
- *	1. Since tan(-x) = -tan(x), we need only to consider positive x.
- *	2. Callers must return tan(-0) = -0 without calling here since our
- *	   odd polynomial is not evaluated in a way that preserves -0.
- *	   Callers may do the optimization tan(x) ~ x for tiny x.
- *	3. tan(x) is approximated by a odd polynomial of degree 27 on
- *	   [0,0.67434]
- *		  	         3             27
- *	   	tan(x) ~ x + T1*x + ... + T13*x
- *	   where
- *
- * 	        |tan(x)         2     4            26   |     -59.2
- * 	        |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
- * 	        |  x 					|
- *
- *	   Note: tan(x+y) = tan(x) + tan'(x)*y
- *		          ~ tan(x) + (1+x*x)*y
- *	   Therefore, for better accuracy in computing tan(x+y), let
- *		     3      2      2       2       2
- *		r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
- *	   then
- *		 		    3    2
- *		tan(x+y) = x + (T1*x + (x *(r+y)+y))
- *
- *      4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
- *		tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
- *		       = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
- */
-
-#include "math_private.h"
-static const double xxx[] = {
-		 3.33333333333334091986e-01,	/* 3FD55555, 55555563 */
-		 1.33333333333201242699e-01,	/* 3FC11111, 1110FE7A */
-		 5.39682539762260521377e-02,	/* 3FABA1BA, 1BB341FE */
-		 2.18694882948595424599e-02,	/* 3F9664F4, 8406D637 */
-		 8.86323982359930005737e-03,	/* 3F8226E3, E96E8493 */
-		 3.59207910759131235356e-03,	/* 3F6D6D22, C9560328 */
-		 1.45620945432529025516e-03,	/* 3F57DBC8, FEE08315 */
-		 5.88041240820264096874e-04,	/* 3F4344D8, F2F26501 */
-		 2.46463134818469906812e-04,	/* 3F3026F7, 1A8D1068 */
-		 7.81794442939557092300e-05,	/* 3F147E88, A03792A6 */
-		 7.14072491382608190305e-05,	/* 3F12B80F, 32F0A7E9 */
-		-1.85586374855275456654e-05,	/* BEF375CB, DB605373 */
-		 2.59073051863633712884e-05,	/* 3EFB2A70, 74BF7AD4 */
-/* one */	 1.00000000000000000000e+00,	/* 3FF00000, 00000000 */
-/* pio4 */	 7.85398163397448278999e-01,	/* 3FE921FB, 54442D18 */
-/* pio4lo */	 3.06161699786838301793e-17	/* 3C81A626, 33145C07 */
-};
-#define	one	xxx[13]
-#define	pio4	xxx[14]
-#define	pio4lo	xxx[15]
-#define	T	xxx
-/* INDENT ON */
-
-double
-__kernel_tan(double x, double y, int iy) {
-	double z, r, v, w, s;
-	int32_t ix, hx;
-
-	GET_HIGH_WORD(hx,x);
-	ix = hx & 0x7fffffff;			/* high word of |x| */
-	if (ix >= 0x3FE59428) {	/* |x| >= 0.6744 */
-		if (hx < 0) {
-			x = -x;
-			y = -y;
-		}
-		z = pio4 - x;
-		w = pio4lo - y;
-		x = z + w;
-		y = 0.0;
-	}
-	z = x * x;
-	w = z * z;
-	/*
-	 * Break x^5*(T[1]+x^2*T[2]+...) into
-	 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
-	 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
-	 */
-	r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] +
-		w * T[11]))));
-	v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] +
-		w * T[12])))));
-	s = z * x;
-	r = y + z * (s * (r + v) + y);
-	r += T[0] * s;
-	w = x + r;
-	if (ix >= 0x3FE59428) {
-		v = (double) iy;
-		return (double) (1 - ((hx >> 30) & 2)) *
-			(v - 2.0 * (x - (w * w / (w + v) - r)));
-	}
-	if (iy == 1)
-		return w;
-	else {
-		/*
-		 * if allow error up to 2 ulp, simply return
-		 * -1.0 / (x+r) here
-		 */
-		/* compute -1.0 / (x+r) accurately */
-		double a, t;
-		z = w;
-		SET_LOW_WORD(z,0);
-		v = r - (z - x);	/* z+v = r+x */
-		t = a = -1.0 / w;	/* a = -1.0/w */
-		SET_LOW_WORD(t,0);
-		s = 1.0 + t * z;
-		return t + a * (s + t * v);
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/math_private.h
+++ /dev/null
@@ -1,820 +0,0 @@
-/*
- * ====================================================
- * 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/Endian.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 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 cos fdlibm::cos
-#define sin fdlibm::sin
-#define tan fdlibm::tan
-#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 floorf fdlibm::floorf
-
-/* 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_ */
deleted file mode 100644
--- a/modules/fdlibm/src/s_asinh.cpp
+++ /dev/null
@@ -1,57 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_atan.cpp
+++ /dev/null
@@ -1,119 +0,0 @@
-/* @(#)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;
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_cbrt.cpp
+++ /dev/null
@@ -1,112 +0,0 @@
-/* @(#)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 divison 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);
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_ceil.cpp
+++ /dev/null
@@ -1,72 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_ceilf.cpp
+++ /dev/null
@@ -1,51 +0,0 @@
-/* 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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_copysign.cpp
+++ /dev/null
@@ -1,32 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_cos.cpp
+++ /dev/null
@@ -1,84 +0,0 @@
-/* @(#)s_cos.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$");
-
-/* cos(x)
- * Return cosine function of x.
- *
- * kernel function:
- *	__kernel_sin		... sine function on [-pi/4,pi/4]
- *	__kernel_cos		... cosine function on [-pi/4,pi/4]
- *	__ieee754_rem_pio2	... argument reduction routine
- *
- * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
- *	in [-pi/4 , +pi/4], and let n = k mod 4.
- *	We have
- *
- *          n        sin(x)      cos(x)        tan(x)
- *     ----------------------------------------------------------
- *	    0	       S	   C		 T
- *	    1	       C	  -S		-1/T
- *	    2	      -S	  -C		 T
- *	    3	      -C	   S		-1/T
- *     ----------------------------------------------------------
- *
- * Special cases:
- *      Let trig be any of sin, cos, or tan.
- *      trig(+-INF)  is NaN, with signals;
- *      trig(NaN)    is that NaN;
- *
- * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded
- */
-
-#include <float.h>
-
-#define INLINE_REM_PIO2
-#include "math_private.h"
-#include "e_rem_pio2.cpp"
-
-double
-cos(double x)
-{
-	double y[2],z=0.0;
-	int32_t n, ix;
-
-    /* High word of x. */
-	GET_HIGH_WORD(ix,x);
-
-    /* |x| ~< pi/4 */
-	ix &= 0x7fffffff;
-	if(ix <= 0x3fe921fb) {
-	    if(ix<0x3e46a09e)			/* if x < 2**-27 * sqrt(2) */
-		if(((int)x)==0) return 1.0;	/* generate inexact */
-	    return __kernel_cos(x,z);
-	}
-
-    /* cos(Inf or NaN) is NaN */
-	else if (ix>=0x7ff00000) return x-x;
-
-    /* argument reduction needed */
-	else {
-	    n = __ieee754_rem_pio2(x,y);
-	    switch(n&3) {
-		case 0: return  __kernel_cos(y[0],y[1]);
-		case 1: return -__kernel_sin(y[0],y[1],1);
-		case 2: return -__kernel_cos(y[0],y[1]);
-		default:
-		        return  __kernel_sin(y[0],y[1],1);
-	    }
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_expm1.cpp
+++ /dev/null
@@ -1,220 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_fabs.cpp
+++ /dev/null
@@ -1,30 +0,0 @@
- /* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_floor.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_floorf.cpp
+++ /dev/null
@@ -1,60 +0,0 @@
-/* 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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_log1p.cpp
+++ /dev/null
@@ -1,175 +0,0 @@
-/* @(#)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);
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_scalbn.cpp
+++ /dev/null
@@ -1,60 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_sin.cpp
+++ /dev/null
@@ -1,84 +0,0 @@
-/* @(#)s_sin.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$");
-
-/* sin(x)
- * Return sine function of x.
- *
- * kernel function:
- *	__kernel_sin		... sine function on [-pi/4,pi/4]
- *	__kernel_cos		... cose function on [-pi/4,pi/4]
- *	__ieee754_rem_pio2	... argument reduction routine
- *
- * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
- *	in [-pi/4 , +pi/4], and let n = k mod 4.
- *	We have
- *
- *          n        sin(x)      cos(x)        tan(x)
- *     ----------------------------------------------------------
- *	    0	       S	   C		 T
- *	    1	       C	  -S		-1/T
- *	    2	      -S	  -C		 T
- *	    3	      -C	   S		-1/T
- *     ----------------------------------------------------------
- *
- * Special cases:
- *      Let trig be any of sin, cos, or tan.
- *      trig(+-INF)  is NaN, with signals;
- *      trig(NaN)    is that NaN;
- *
- * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded
- */
-
-#include <float.h>
-
-#define INLINE_REM_PIO2
-#include "math_private.h"
-#include "e_rem_pio2.cpp"
-
-double
-sin(double x)
-{
-	double y[2],z=0.0;
-	int32_t n, ix;
-
-    /* High word of x. */
-	GET_HIGH_WORD(ix,x);
-
-    /* |x| ~< pi/4 */
-	ix &= 0x7fffffff;
-	if(ix <= 0x3fe921fb) {
-	    if(ix<0x3e500000)			/* |x| < 2**-26 */
-	       {if((int)x==0) return x;}	/* generate inexact */
-	    return __kernel_sin(x,z,0);
-	}
-
-    /* sin(Inf or NaN) is NaN */
-	else if (ix>=0x7ff00000) return x-x;
-
-    /* argument reduction needed */
-	else {
-	    n = __ieee754_rem_pio2(x,y);
-	    switch(n&3) {
-		case 0: return  __kernel_sin(y[0],y[1],1);
-		case 1: return  __kernel_cos(y[0],y[1]);
-		case 2: return -__kernel_sin(y[0],y[1],1);
-		default:
-			return -__kernel_cos(y[0],y[1]);
-	    }
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_tan.cpp
+++ /dev/null
@@ -1,78 +0,0 @@
-/* @(#)s_tan.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$");
-
-/* tan(x)
- * Return tangent function of x.
- *
- * kernel function:
- *	__kernel_tan		... tangent function on [-pi/4,pi/4]
- *	__ieee754_rem_pio2	... argument reduction routine
- *
- * Method.
- *      Let S,C and T denote the sin, cos and tan respectively on
- *	[-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
- *	in [-pi/4 , +pi/4], and let n = k mod 4.
- *	We have
- *
- *          n        sin(x)      cos(x)        tan(x)
- *     ----------------------------------------------------------
- *	    0	       S	   C		 T
- *	    1	       C	  -S		-1/T
- *	    2	      -S	  -C		 T
- *	    3	      -C	   S		-1/T
- *     ----------------------------------------------------------
- *
- * Special cases:
- *      Let trig be any of sin, cos, or tan.
- *      trig(+-INF)  is NaN, with signals;
- *      trig(NaN)    is that NaN;
- *
- * Accuracy:
- *	TRIG(x) returns trig(x) nearly rounded
- */
-
-#include <float.h>
-
-#define INLINE_REM_PIO2
-#include "math_private.h"
-#include "e_rem_pio2.cpp"
-
-double
-tan(double x)
-{
-	double y[2],z=0.0;
-	int32_t n, ix;
-
-    /* High word of x. */
-	GET_HIGH_WORD(ix,x);
-
-    /* |x| ~< pi/4 */
-	ix &= 0x7fffffff;
-	if(ix <= 0x3fe921fb) {
-	    if(ix<0x3e400000)			/* x < 2**-27 */
-		if((int)x==0) return x;		/* generate inexact */
-	    return __kernel_tan(x,z,1);
-	}
-
-    /* tan(Inf or NaN) is NaN */
-	else if (ix>=0x7ff00000) return x-x;		/* NaN */
-
-    /* argument reduction needed */
-	else {
-	    n = __ieee754_rem_pio2(x,y);
-	    return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /*   1 -- n even
-							-1 -- n odd */
-	}
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_tanh.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/* @(#)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;
-}
deleted file mode 100644
--- a/modules/fdlibm/src/s_trunc.cpp
+++ /dev/null
@@ -1,62 +0,0 @@
-/* @(#)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;
-}