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