Merge "Implements some math functions for faster performance"

3 files changed, 1017 insertions, 129 deletions

diff --git a/NOTICE b/NOTICE
index 951e506..5136b4b 100644
--- a/NOTICE
+++ b/NOTICE
@@ -104,3 +104,14 @@ WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 
 See W3C License http://www.w3.org/Consortium/Legal/ for more details.
+
+   =========================================================================
+   ==  NOTICE file for the fdlibm License.                                   ==
+   =========================================================================
+
+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.
diff --git a/luni/src/main/java/java/lang/StrictMath.java b/luni/src/main/java/java/lang/StrictMath.java
index f409c06..2e848f2 100644
--- a/luni/src/main/java/java/lang/StrictMath.java
+++ b/luni/src/main/java/java/lang/StrictMath.java
@@ -15,6 +15,18 @@
  *  limitations under the License.
  */
 
+/*
+ * acos, asin, atan, cosh, sinh, tanh, exp, expm1, log, log10, log1p, and cbrt
+ * have been implemented with the following license.
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 package java.lang;
 
 /**
@@ -102,6 +114,21 @@ public final class StrictMath {
         return Math.abs(l);
     }
 
+    private static final double PIO2_HI = 1.57079632679489655800e+00;
+    private static final double PIO2_LO = 6.12323399573676603587e-17;
+    private static final double PS0 = 1.66666666666666657415e-01;
+    private static final double PS1 = -3.25565818622400915405e-01;
+    private static final double PS2 = 2.01212532134862925881e-01;
+    private static final double PS3 = -4.00555345006794114027e-02;
+    private static final double PS4 = 7.91534994289814532176e-04;
+    private static final double PS5 = 3.47933107596021167570e-05;
+    private static final double QS1 = -2.40339491173441421878e+00;
+    private static final double QS2 = 2.02094576023350569471e+00;
+    private static final double QS3 = -6.88283971605453293030e-01;
+    private static final double QS4 = 7.70381505559019352791e-02;
+    private static final double HUGE = 1.000e+300;
+    private static final double PIO4_HI = 7.85398163397448278999e-01;
+
     /**
      * Returns the closest double approximation of the arc cosine of the
      * argument within the range {@code [0..pi]}.
@@ -113,11 +140,62 @@ public final class StrictMath {
      * <li>{@code acos(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value to compute arc cosine of.
      * @return the arc cosine of the argument.
      */
-    public static native double acos(double d);
+    public static double acos(double x) {
+        double z, p, q, r, w, s, c, df;
+        int hx, ix;
+        final long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >>> 32);
+        ix = hx & 0x7fffffff;
+        if (ix >= 0x3ff00000) { /* |x| >= 1 */
+            if ((((ix - 0x3ff00000) | ((int) bits))) == 0) { /* |x|==1 */
+                if (hx > 0) {
+                    return 0.0; /* ieee_acos(1) = 0 */
+                } else {
+                    return 3.14159265358979311600e+00 + 2.0 * PIO2_LO; /* ieee_acos(-1)= pi */
+                }
+            }
+            return (x - x) / (x - x); /* ieee_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 = 1.00000000000000000000e+00 + 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 = (1.00000000000000000000e+00 + x) * 0.5;
+            p = z * (PS0 + z
+                    * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5)))));
+            q = 1.00000000000000000000e+00 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
+            s = StrictMath.sqrt(z);
+            r = p / q;
+            w = r * s - PIO2_LO;
+            return 3.14159265358979311600e+00 - 2.0 * (s + w);
+        } else { /* x > 0.5 */
+            z = (1.00000000000000000000e+00 - x) * 0.5;
+            s = StrictMath.sqrt(z);
+            df = s;
+            df = Double.longBitsToDouble(
+                    Double.doubleToRawLongBits(df) & 0xffffffffL << 32);
+            c = (z - df * df) / (s + df);
+            p = z * (PS0 + z
+                    * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5)))));
+            q = 1.00000000000000000000e+00 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
+            r = p / q;
+            w = r * s + c;
+            return 2.0 * (df + w);
+        }
+    }
 
     /**
      * Returns the closest double approximation of the arc sine of the argument
@@ -130,11 +208,75 @@ public final class StrictMath {
      * <li>{@code asin(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose arc sine has to be computed.
      * @return the arc sine of the argument.
      */
-    public static native double asin(double d);
+    public static double asin(double x) {
+        double t, w, p, q, c, r, s;
+        int hx, ix;
+        final long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >>> 32);
+        ix = hx & 0x7fffffff;
+        if (ix >= 0x3ff00000) { /* |x|>= 1 */
+            if ((((ix - 0x3ff00000) | ((int) bits))) == 0) {
+                /* ieee_asin(1)=+-pi/2 with inexact */
+                return x * PIO2_HI + x * PIO2_LO;
+            }
+            return (x - x) / (x - x); /* ieee_asin(|x|>1) is NaN */
+        } else if (ix < 0x3fe00000) { /* |x|<0.5 */
+            if (ix < 0x3e400000) { /* if |x| < 2**-27 */
+                if (HUGE + x > 1.00000000000000000000e+00) {
+                    return x;/* return x with inexact if x!=0 */
+                }
+            } else {
+                t = x * x;
+                p = t * (PS0 + t
+                        * (PS1 + t * (PS2 + t * (PS3 + t * (PS4 + t * PS5)))));
+                q = 1.00000000000000000000e+00 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4)));
+                w = p / q;
+                return x + x * w;
+            }
+        }
+        /* 1> |x|>= 0.5 */
+        w = 1.00000000000000000000e+00 - Math.abs(x);
+        t = w * 0.5;
+        p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t * (PS4 + t * PS5)))));
+        q = 1.00000000000000000000e+00 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4)));
+        s = StrictMath.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;
+            w = Double.longBitsToDouble(
+                    Double.doubleToRawLongBits(w) & 0xffffffffL << 32);
+            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);
+        }
+        return (hx > 0) ? t : -t;
+    }
+
+    private static final double[] ATANHI = { 4.63647609000806093515e-01,
+            7.85398163397448278999e-01, 9.82793723247329054082e-01,
+            1.57079632679489655800e+00 };
+    private static final double[] ATANLO = { 2.26987774529616870924e-17,
+            3.06161699786838301793e-17, 1.39033110312309984516e-17,
+            6.12323399573676603587e-17 };
+    private static final double AT0 = 3.33333333333329318027e-01;
+    private static final double AT1 = -1.99999999998764832476e-01;
+    private static final double AT2 = 1.42857142725034663711e-01;
+    private static final double AT3 = -1.11111104054623557880e-01;
+    private static final double AT4 = 9.09088713343650656196e-02;
+    private static final double AT5 = -7.69187620504482999495e-02;
+    private static final double AT6 = 6.66107313738753120669e-02;
+    private static final double AT7= -5.83357013379057348645e-02;
+    private static final double AT8 = 4.97687799461593236017e-02;
+    private static final double AT9 = -3.65315727442169155270e-02;
+    private static final double AT10 = 1.62858201153657823623e-02;
 
     /**
      * Returns the closest double approximation of the arc tangent of the
@@ -149,11 +291,73 @@ public final class StrictMath {
      * <li>{@code atan(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose arc tangent has to be computed.
      * @return the arc tangent of the argument.
      */
-    public static native double atan(double d);
+    public static double atan(double x) {
+       double w, s1, s2, z;
+        int ix, hx, id;
+
+        final long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >>> 32);
+        ix = hx & 0x7fffffff;
+        if (ix >= 0x44100000) { /* if |x| >= 2^66 */
+            if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (((int) bits) != 0))) {
+                return x + x; /* NaN */
+            }
+            if (hx > 0) {
+                return ATANHI[3] + ATANLO[3];
+            } else {
+                return -ATANHI[3] - ATANLO[3];
+            }
+        }
+        if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
+            if (ix < 0x3e200000) { /* |x| < 2^-29 */
+                if (HUGE + x > 1.00000000000000000000e+00) {
+                    return x; /* raise inexact */
+                }
+            }
+            id = -1;
+        } else {
+            x = Math.abs(x);
+            if (ix < 0x3ff30000) { /* |x| < 1.1875 */
+                if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
+                    id = 0;
+                    x = (2.0 * x - 1.00000000000000000000e+00) / (2.0 + x);
+                } else { /* 11/16<=|x|< 19/16 */
+                    id = 1;
+                    x = (x - 1.00000000000000000000e+00) / (x + 1.00000000000000000000e+00);
+                }
+            } else {
+                if (ix < 0x40038000) { /* |x| < 2.4375 */
+                    id = 2;
+                    x = (x - 1.5) / (1.00000000000000000000e+00 + 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 * (AT0 + w * (AT2 + w
+                * (AT4 + w * (AT6 + w * (AT8 + w * AT10)))));
+        s2 = w * (AT1 + w * (AT3 + w * (AT5 + w * (AT7 + w * AT9))));
+        if (id < 0) {
+            return x - x * (s1 + s2);
+        } else {
+            z = ATANHI[id] - ((x * (s1 + s2) - ATANLO[id]) - x);
+            return (hx < 0) ? -z : z;
+        }
+    }
+
+    private static final double PI_O_4 = 7.8539816339744827900E-01;
+    private static final double PI_O_2 = 1.5707963267948965580E+00;
+    private static final double PI_LO = 1.2246467991473531772E-16;
 
     /**
      * Returns the closest double approximation of the arc tangent of
@@ -192,7 +396,108 @@ public final class StrictMath {
      *            the denominator of the value whose atan has to be computed.
      * @return the arc tangent of {@code y/x}.
      */
-    public static native double atan2(double y, double x);
+    public static double atan2(double y, double x) {
+        double z;
+        int k, m, hx, hy, ix, iy;
+        int lx, ly; // watch out, should be unsigned
+
+        final long yBits = Double.doubleToRawLongBits(y);
+        final long xBits = Double.doubleToRawLongBits(x);
+
+        hx = (int) (xBits >>> 32); // __HI(x);
+        ix = hx & 0x7fffffff;
+        lx = (int) xBits; // __LO(x);
+        hy = (int) (yBits >>> 32); // __HI(y);
+        iy = hy & 0x7fffffff;
+        ly = (int) yBits; // __LO(y);
+        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 StrictMath.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; /* ieee_atan(+-0,+anything)=+-0 */
+                case 2:
+                    return 3.14159265358979311600e+00 + TINY;/* ieee_atan(+0,-anything) = pi */
+                case 3:
+                    return -3.14159265358979311600e+00 - TINY;/* ieee_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;/* ieee_atan(+INF,+INF) */
+                    case 1:
+                        return -PI_O_4 - TINY;/* ieee_atan(-INF,+INF) */
+                    case 2:
+                        return 3.0 * PI_O_4 + TINY;/* ieee_atan(+INF,-INF) */
+                    case 3:
+                        return -3.0 * PI_O_4 - TINY;/* ieee_atan(-INF,-INF) */
+                }
+            } else {
+                switch (m) {
+                    case 0:
+                        return 0.0; /* ieee_atan(+...,+INF) */
+                    case 1:
+                        return -0.0; /* ieee_atan(-...,+INF) */
+                    case 2:
+                        return 3.14159265358979311600e+00 + TINY; /* ieee_atan(+...,-INF) */
+                    case 3:
+                        return -3.14159265358979311600e+00 - TINY; /* ieee_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) {
+            z = PI_O_2 + 0.5 * PI_LO; /* |y/x| > 2**60 */
+        } else if (hx < 0 && k < -60) {
+            z = 0.0; /* |y|/x < -2**60 */
+        } else {
+            z = StrictMath.atan(Math.abs(y / x)); /* safe to do y/x */
+        }
+
+        switch (m) {
+            case 0:
+                return z; /* ieee_atan(+,+) */
+            case 1:
+                // __HI(z) ^= 0x80000000;
+                z = Double.longBitsToDouble(
+                        Double.doubleToRawLongBits(z) ^ (0x80000000L << 32));
+                return z; /* ieee_atan(-,+) */
+            case 2:
+                return 3.14159265358979311600e+00 - (z - PI_LO);/* ieee_atan(+,-) */
+            default: /* case 3 */
+                return (z - PI_LO) - 3.14159265358979311600e+00;/* ieee_atan(-,-) */
+        }
+    }
+
+    private static final int B1 = 715094163;
+    private static final int B2 = 696219795;
+    private static final double C = 5.42857142857142815906e-01;
+    private static final double D = -7.05306122448979611050e-01;
+    private static final double CBRTE = 1.41428571428571436819e+00;
+    private static final double F = 1.60714285714285720630e+00;
+    private static final double G = 3.57142857142857150787e-01;
 
     /**
      * Returns the closest double approximation of the cube root of the
@@ -207,11 +512,79 @@ public final class StrictMath {
      * <li>{@code cbrt(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose cube root has to be computed.
      * @return the cube root of the argument.
      */
-    public static native double cbrt(double d);
+    public static double cbrt(double x) {
+        if (x < 0) {
+            return -cbrt(-x);
+        }
+        int hx;
+        double r, s, w;
+        int sign; // caution: should be unsigned
+        long bits = Double.doubleToRawLongBits(x);
+
+        hx = (int) (bits >>> 32);
+        sign = hx & 0x80000000; /* sign= sign(x) */
+        hx ^= sign;
+        if (hx >= 0x7ff00000) {
+            return (x + x); /* ieee_cbrt(NaN,INF) is itself */
+        }
+
+        if ((hx | ((int) bits)) == 0) {
+            return x; /* ieee_cbrt(0) is itself */
+        }
+
+        // __HI(x) = hx; /* x <- |x| */
+        bits &= 0x00000000ffffffffL;
+        bits |= ((long) hx << 32);
+
+        long tBits = Double.doubleToRawLongBits(0.0) & 0x00000000ffffffffL;
+        double t = 0.0;
+        /* rough cbrt to 5 bits */
+        if (hx < 0x00100000) { /* subnormal number */
+            // __HI(t)=0x43500000; /*set t= 2**54*/
+            tBits |= 0x43500000L << 32;
+            t = Double.longBitsToDouble(tBits);
+            t *= x;
+
+            // __HI(t)=__HI(t)/3+B2;
+            tBits = Double.doubleToRawLongBits(t);
+            long tBitsHigh = tBits >> 32;
+            tBits &= 0x00000000ffffffffL;
+            tBits |= ((tBitsHigh / 3) + B2) << 32;
+            t = Double.longBitsToDouble(tBits);
+
+        } else {
+            // __HI(t)=hx/3+B1;
+            tBits |= ((long) ((hx / 3) + B1)) << 32;
+            t = Double.longBitsToDouble(tBits);
+        }
+
+        /* new cbrt to 23 bits, may be implemented in single precision */
+        r = t * t / x;
+        s = C + r * t;
+        t *= G + F / (s + CBRTE + D / s);
+
+        /* chopped to 20 bits and make it larger than ieee_cbrt(x) */
+        tBits = Double.doubleToRawLongBits(t);
+        tBits &= 0xFFFFFFFFL << 32;
+        tBits += 0x00000001L << 32;
+        t = Double.longBitsToDouble(tBits);
+
+        /* one step newton iteration to 53 bits with error less than 0.667 ulps */
+        s = t * t; /* t*t is exact */
+        r = x / s;
+        w = t + t;
+        r = (r - t) / (w + r); /* r-s is exact */
+        t = t + t * r;
+
+        /* retore the sign bit */
+        tBits = Double.doubleToRawLongBits(t);
+        tBits |= ((long) sign) << 32;
+        return Double.longBitsToDouble(tBits);
+    }
 
     /**
      * Returns the double conversion of the most negative (closest to negative
@@ -229,6 +602,8 @@ public final class StrictMath {
      */
     public static native double ceil(double d);
 
+    private static final long ONEBITS = Double.doubleToRawLongBits(1.00000000000000000000e+00)
+            & 0x00000000ffffffffL;
 
     /**
      * Returns the closest double approximation of the hyperbolic cosine of the
@@ -241,11 +616,54 @@ public final class StrictMath {
      * <li>{@code cosh(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose hyperbolic cosine has to be computed.
      * @return the hyperbolic cosine of the argument.
      */
-    public static native double cosh(double d);
+    public static double cosh(double x) {
+        double t, w;
+        int ix;
+        final long bits = Double.doubleToRawLongBits(x);
+        ix = (int) (bits >>> 32) & 0x7fffffff;
+
+        /* x is INF or NaN */
+        if (ix >= 0x7ff00000) {
+            return x * x;
+        }
+
+        /* |x| in [0,0.5*ln2], return 1+ieee_expm1(|x|)^2/(2*ieee_exp(|x|)) */
+        if (ix < 0x3fd62e43) {
+            t = expm1(Math.abs(x));
+            w = 1.00000000000000000000e+00 + t;
+            if (ix < 0x3c800000)
+                return w; /* ieee_cosh(tiny) = 1 */
+            return 1.00000000000000000000e+00 + (t * t) / (w + w);
+        }
+
+        /* |x| in [0.5*ln2,22], return (ieee_exp(|x|)+1/ieee_exp(|x|)/2; */
+        if (ix < 0x40360000) {
+            t = exp(Math.abs(x));
+            return 0.5 * t + 0.5 / t;
+        }
+
+        /* |x| in [22, ieee_log(maxdouble)] return half*ieee_exp(|x|) */
+        if (ix < 0x40862E42) {
+            return 0.5 * exp(Math.abs(x));
+        }
+
+        /* |x| in [log(maxdouble), overflowthresold] */
+        final long lx = ((ONEBITS >>> 29) + ((int) bits)) & 0x00000000ffffffffL;
+        // watch out: lx should be an unsigned int
+        // lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
+        if (ix < 0x408633CE || (ix == 0x408633ce) && (lx <= 0x8fb9f87dL)) {
+            w = exp(0.5 * Math.abs(x));
+            t = 0.5 * w;
+            return t * w;
+        }
+
+        /* |x| > overflowthresold, ieee_cosh(x) overflow */
+        return HUGE * HUGE;
+    }
 
     /**
      * Returns the closest double approximation of the cosine of the argument.
@@ -263,6 +681,19 @@ public final class StrictMath {
      */
     public static native double cos(double d);
 
+    private static final double TWON24 = 5.96046447753906250000e-08;
+    private static final double TWO54 = 1.80143985094819840000e+16,
+            TWOM54 = 5.55111512312578270212e-17;
+    private static final double TWOM1000 = 9.33263618503218878990e-302;
+    private static final double O_THRESHOLD = 7.09782712893383973096e+02;
+    private static final double U_THRESHOLD = -7.45133219101941108420e+02;
+    private static final double INVLN2 = 1.44269504088896338700e+00;
+    private static final double P1 = 1.66666666666666019037e-01;
+    private static final double P2 = -2.77777777770155933842e-03;
+    private static final double P3 = 6.61375632143793436117e-05;
+    private static final double P4 = -1.65339022054652515390e-06;
+    private static final double P5 = 4.13813679705723846039e-08;
+
     /**
      * Returns the closest double approximation of the raising "e" to the power
      * of the argument.
@@ -274,11 +705,88 @@ public final class StrictMath {
      * <li>{@code exp(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose exponential has to be computed.
      * @return the exponential of the argument.
      */
-    public static native double exp(double d);
+    public static double exp(double x) {
+        double y, c, t;
+        double hi = 0, lo = 0;
+        int k = 0, xsb;
+        int hx; // should be unsigned, be careful!
+        final long bits = Double.doubleToRawLongBits(x);
+        int lowBits = (int) bits;
+        int highBits = (int) (bits >>> 32);
+        hx = highBits & 0x7fffffff;
+        xsb = (highBits >>> 31) & 1;
+
+        /* filter out non-finite argument */
+        if (hx >= 0x40862E42) { /* if |x|>=709.78... */
+            if (hx >= 0x7ff00000) {
+                if (((hx & 0xfffff) | lowBits) != 0) {
+                    return x + x; /* NaN */
+                } else {
+                    return (xsb == 0) ? x : 0.0; /* ieee_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 - ((xsb == 0) ? 6.93147180369123816490e-01 :
+                        -6.93147180369123816490e-01); // LN2HI[xsb];
+                lo = (xsb == 0) ? 1.90821492927058770002e-10 :
+                        -1.90821492927058770002e-10; // LN2LO[xsb];
+                k = 1 - xsb - xsb;
+            } else {
+                k = (int) (INVLN2 * x + ((xsb == 0) ? 0.5 : -0.5 ));//halF[xsb]);
+                t = k;
+                hi = x - t * 6.93147180369123816490e-01; //ln2HI[0]; /* t*ln2HI is exact here */
+                lo = t * 1.90821492927058770002e-10; //ln2LO[0];
+            }
+            x = hi - lo;
+        } else if (hx < 0x3e300000) { /* when |x|<2**-28 */
+            if (HUGE + x > 1.00000000000000000000e+00)
+                return 1.00000000000000000000e+00 + x;/* trigger inexact */
+        } else {
+            k = 0;
+        }
+
+        /* x is now in primary range */
+        t = x * x;
+        c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
+        if (k == 0) {
+            return 1.00000000000000000000e+00 - ((x * c) / (c - 2.0) - x);
+        } else {
+            y = 1.00000000000000000000e+00 - ((lo - (x * c) / (2.0 - c)) - hi);
+        }
+        long yBits = Double.doubleToRawLongBits(y);
+        if (k >= -1021) {
+            yBits += ((long) (k << 20)) << 32; /* add k to y's exponent */
+            return Double.longBitsToDouble(yBits);
+        } else {
+            yBits += ((long) ((k + 1000) << 20)) << 32;/* add k to y's exponent */
+            return Double.longBitsToDouble(yBits) * TWOM1000;
+        }
+    }
+
+    private static final double TINY = 1.0e-300;
+    private static final double LN2_HI = 6.93147180369123816490e-01;
+    private static final double LN2_LO = 1.90821492927058770002e-10;
+    private static final double Q1 = -3.33333333333331316428e-02;
+    private static final double Q2 = 1.58730158725481460165e-03;
+    private static final double Q3 = -7.93650757867487942473e-05;
+    private static final double Q4 = 4.00821782732936239552e-06;
+    private static final double Q5 = -2.01099218183624371326e-07;
 
     /**
      * Returns the closest double approximation of <i>{@code e}</i><sup>
@@ -295,17 +803,124 @@ public final class StrictMath {
      * <li>{@code expm1(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value to compute the <i>{@code e}</i><sup>{@code d}</sup>
      *            {@code - 1} of.
-     * @return the <i>{@code e}</i><sup>{@code d}</sup>{@code - 1} value
-     *         of the argument.
+     * @return the <i>{@code e}</i><sup>{@code d}</sup>{@code - 1} value of the
+     *         argument.
      */
-    public static native double expm1(double d);
+    public static double expm1(double x) {
+        double y, hi, lo, t, e, hxs, hfx, r1, c = 0.0;
+        int k, xsb;
+        long yBits = 0;
+        final long bits = Double.doubleToRawLongBits(x);
+        int highBits = (int) (bits >>> 32);
+        int lowBits = (int) (bits);
+        int hx = highBits & 0x7fffffff; // caution: should be unsigned!
+        xsb = highBits & 0x80000000; /* sign bit of x */
+        y = xsb == 0 ? x : -x; /* y = |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) {
+                    if (((hx & 0xfffff) | lowBits) != 0) {
+                        return x + x; /* NaN */
+                    } else {
+                        return (xsb == 0) ? x : -1.0;/* ieee_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 - 1.00000000000000000000e+00; /* 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 = (int) (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;
+            }
+            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));
+            return x; // inexact flag is not set, but Java ignors this flag
+                      // anyway
+        } else {
+            k = 0;
+        }
+
+        /* x is now in primary range */
+        hfx = 0.5 * x;
+        hxs = x * hfx;
+        r1 = 1.00000000000000000000e+00 + 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 {
+            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 1.00000000000000000000e+00 + 2.0 * (x - e);
+                }
+            }
+
+            if (k <= -2 || k > 56) { /* suffice to return ieee_exp(x)-1 */
+                y = 1.00000000000000000000e+00 - (e - x);
+                yBits = Double.doubleToRawLongBits(y);
+                yBits += (((long) k) << 52); /* add k to y's exponent */
+                return Double.longBitsToDouble(yBits) - 1.00000000000000000000e+00;
+            }
+
+            long tBits = Double.doubleToRawLongBits(1.00000000000000000000e+00) & 0x00000000ffffffffL;
+
+            if (k < 20) {
+                tBits |= (((long) 0x3ff00000) - (0x200000 >> k)) << 32;
+                y = Double.longBitsToDouble(tBits) - (e - x);
+                yBits = Double.doubleToRawLongBits(y);
+                yBits += (((long) k) << 52); /* add k to y's exponent */
+                return Double.longBitsToDouble(yBits);
+            } else {
+                tBits |= ((((long) 0x3ff) - k) << 52); /* 2^-k */
+                y = x - (e + Double.longBitsToDouble(tBits));
+                y += 1.00000000000000000000e+00;
+                yBits = Double.doubleToRawLongBits(y);
+                yBits += (((long) k) << 52); /* add k to y's exponent */
+                return Double.longBitsToDouble(yBits);
+            }
+        }
+    }
 
     /**
-     * Returns the double conversion of the most positive (closest to
-     * positive infinity) integer less than or equal to the argument.
+     * Returns the double conversion of the most positive (closest to positive
+     * infinity) integer less than or equal to the argument.
      * <p>
      * Special cases:
      * <ul>
@@ -319,9 +934,9 @@ public final class StrictMath {
     public static native double floor(double d);
 
     /**
-     * Returns {@code sqrt(}<i>{@code x}</i><sup>{@code 2}</sup>{@code +}
-     * <i> {@code y}</i><sup>{@code 2}</sup>{@code )}. The final result is
-     * without medium underflow or overflow.
+     * Returns {@code sqrt(}<i>{@code x}</i><sup>{@code 2}</sup>{@code +} <i>
+     * {@code y}</i><sup>{@code 2}</sup>{@code )}. The final result is without
+     * medium underflow or overflow.
      * <p>
      * Special cases:
      * <ul>
@@ -369,6 +984,14 @@ public final class StrictMath {
      */
     public static native double IEEEremainder(double x, double y);
 
+    private static final double LG1 = 6.666666666666735130e-01;
+    private static final double LG2 = 3.999999999940941908e-01;
+    private static final double LG3 = 2.857142874366239149e-01;
+    private static final double LG4 = 2.222219843214978396e-01;
+    private static final double LG5 = 1.818357216161805012e-01;
+    private static final double LG6 = 1.531383769920937332e-01;
+    private static final double LG7 = 1.479819860511658591e-01;
+
     /**
      * Returns the closest double approximation of the natural logarithm of the
      * argument.
@@ -383,11 +1006,95 @@ public final class StrictMath {
      * <li>{@code log(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose log has to be computed.
      * @return the natural logarithm of the argument.
      */
-    public static native double log(double d);
+    public static double log(double x) {
+        double hfsq, f, s, z, R, w, t1, t2, dk;
+        int hx, i, j, k = 0;
+        int lx; // watch out, should be unsigned
+
+        long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >>> 32); /* high word of x */
+        lx = (int) bits; /* low word of x */
+
+        if (hx < 0x00100000) { /* x < 2**-1022 */
+            if (((hx & 0x7fffffff) | lx) == 0) {
+                return -TWO54 / 0.0; /* ieee_log(+-0)=-inf */
+            }
+
+            if (hx < 0) {
+                return (x - x) / 0.0; /* ieee_log(-#) = NaN */
+            }
+
+            k -= 54;
+            x *= TWO54; /* subnormal number, scale up x */
+            bits = Double.doubleToRawLongBits(x);
+            hx = (int) (bits >>> 32); /* high word of x */
+        }
+
+        if (hx >= 0x7ff00000) {
+            return x + x;
+        }
+
+        k += (hx >> 20) - 1023;
+        hx &= 0x000fffff;
+        bits &= 0x00000000ffffffffL;
+        i = (hx + 0x95f64) & 0x100000;
+        bits |= ((long) hx | (i ^ 0x3ff00000)) << 32; /* normalize x or x/2 */
+        x = Double.longBitsToDouble(bits);
+        k += (i >> 20);
+        f = x - 1.0;
+
+        if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */
+            if (f == 0.0) {
+                if (k == 0) {
+                    return 0.0;
+                } else {
+                    dk = k;
+                }
+                return dk * LN2_HI + dk * LN2_LO;
+            }
+
+            R = f * f * (0.5 - 0.33333333333333333 * f);
+            if (k == 0) {
+                return f - R;
+            } else {
+                dk = k;
+                return dk * LN2_HI - ((R - dk * LN2_LO) - f);
+            }
+        }
+        s = f / (2.0 + f);
+        dk = 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);
+            }
+        }
+    }
+
+    private static final double IVLN10 = 4.34294481903251816668e-01;
+    private static final double LOG10_2HI = 3.01029995663611771306e-01;
+    private static final double LOG10_2LO = 3.69423907715893078616e-13;
 
     /**
      * Returns the closest double approximation of the base 10 logarithm of the
@@ -403,11 +1110,54 @@ public final class StrictMath {
      * <li>{@code log10(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose base 10 log has to be computed.
-     * @return the natural logarithm of the argument.
+     * @return the the base 10 logarithm of x
      */
-    public static native double log10(double d);
+    public static double log10(double x) {
+        double y, z;
+        int i, k = 0, hx;
+        int lx; // careful: lx should be unsigned!
+        long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >> 32); /* high word of x */
+        lx = (int) bits; /* low word of x */
+        if (hx < 0x00100000) { /* x < 2**-1022 */
+            if (((hx & 0x7fffffff) | lx) == 0) {
+                return -TWO54 / 0.0; /* ieee_log(+-0)=-inf */
+            }
+
+            if (hx < 0) {
+                return (x - x) / 0.0; /* ieee_log(-#) = NaN */
+            }
+
+            k -= 54;
+            x *= TWO54; /* subnormal number, scale up x */
+            bits = Double.doubleToRawLongBits(x);
+            hx = (int) (bits >> 32); /* high word of x */
+        }
+
+        if (hx >= 0x7ff00000) {
+            return x + x;
+        }
+
+        k += (hx >> 20) - 1023;
+        i = (int) (((k & 0x00000000ffffffffL) & 0x80000000) >>> 31);
+        hx = (hx & 0x000fffff) | ((0x3ff - i) << 20);
+        y = k + i;
+        bits &= 0x00000000ffffffffL;
+        bits |= ((long) hx) << 32;
+        x = Double.longBitsToDouble(bits); // __HI(x) = hx;
+        z = y * LOG10_2LO + IVLN10 * log(x);
+        return z + y * LOG10_2HI;
+    }
+
+    private static final double LP1 = 6.666666666666735130e-01,
+            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 */
 
     /**
      * Returns the closest double approximation of the natural logarithm of the
@@ -426,11 +1176,107 @@ public final class StrictMath {
      * <li>{@code log1p(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value to compute the {@code ln(1+d)} of.
      * @return the natural logarithm of the sum of the argument and 1.
      */
-    public static native double log1p(double d);
+
+    public static double log1p(double x) {
+        double hfsq, f = 0.0, c = 0.0, s, z, R, u = 0.0;
+        int k, hx, hu = 0, ax;
+
+        final long bits = Double.doubleToRawLongBits(x);
+        hx = (int) (bits >>> 32); /* high word of x */
+        ax = hx & 0x7fffffff;
+
+        k = 1;
+        if (hx < 0x3FDA827A) { /* x < 0.41422 */
+            if (ax >= 0x3ff00000) { /* x <= -1.0 */
+                if (x == -1.0) {
+                    return -TWO54 / 0.0; /* ieee_log1p(-1)=+inf */
+                } else {
+                    return (x - x) / (x - x); /* ieee_log1p(x<-1)=NaN */
+                }
+            }
+            if (ax < 0x3e200000) {
+                if (TWO54 + x > 0.0 && ax < 0x3c900000) {
+                    return x;
+                } else {
+                    return x - x * x * 0.5;
+                }
+            }
+            if (hx > 0 || hx <= 0xbfd2bec3) {
+                k = 0;
+                f = x;
+                hu = 1;
+            } /* -0.2929<x<0.41422 */
+        }
+
+        if (hx >= 0x7ff00000) {
+            return x + x;
+        }
+
+        if (k != 0) {
+            long uBits;
+            if (hx < 0x43400000) {
+                u = 1.0 + x;
+                uBits = Double.doubleToRawLongBits(u);
+                hu = (int) (uBits >>> 32);
+                k = (hu >> 20) - 1023;
+                c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0);/* correction term */
+                c /= u;
+            } else {
+                uBits = Double.doubleToRawLongBits(x);
+                hu = (int) (uBits >>> 32);
+                k = (hu >> 20) - 1023;
+                c = 0;
+            }
+            hu &= 0x000fffff;
+            if (hu < 0x6a09e) {
+                // __HI(u) = hu|0x3ff00000; /* normalize u */
+                uBits &= 0x00000000ffffffffL;
+                uBits |= ((long) hu | 0x3ff00000) << 32;
+                u = Double.longBitsToDouble(uBits);
+            } else {
+                k += 1;
+                // __HI(u) = hu|0x3fe00000; /* normalize u/2 */
+                uBits &= 0xffffffffL;
+                uBits |= ((long) hu | 0x3fe00000) << 32;
+                u = Double.longBitsToDouble(uBits);
+                hu = (0x00100000 - hu) >> 2;
+            }
+            f = u - 1.0;
+        }
+        hfsq = 0.5 * f * f;
+        if (hu == 0) { /* |f| < 2**-20 */
+            if (f == 0.0) {
+                if (k == 0) {
+                    return 0.0;
+                } 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);
+        }
+    }
 
     /**
      * Returns the most positive (closest to positive infinity) of the two
@@ -453,8 +1299,8 @@ public final class StrictMath {
         if (d1 != d2)
             return Double.NaN;
         /* max( +0.0,-0.0) == +0.0 */
-        if (d1 == 0.0
-                && ((Double.doubleToLongBits(d1) & Double.doubleToLongBits(d2)) & 0x8000000000000000L) == 0)
+        if (d1 == 0.0 &&
+                ((Double.doubleToLongBits(d1) & Double.doubleToLongBits(d2)) & 0x8000000000000000L) == 0)
             return 0.0;
         return d1;
     }
@@ -480,8 +1326,8 @@ public final class StrictMath {
         if (f1 != f2)
             return Float.NaN;
         /* max( +0.0,-0.0) == +0.0 */
-        if (f1 == 0.0f
-                && ((Float.floatToIntBits(f1) & Float.floatToIntBits(f2)) & 0x80000000) == 0)
+        if (f1 == 0.0f &&
+                ((Float.floatToIntBits(f1) & Float.floatToIntBits(f2)) & 0x80000000) == 0)
             return 0.0f;
         return f1;
     }
@@ -523,8 +1369,8 @@ public final class StrictMath {
         if (d1 != d2)
             return Double.NaN;
         /* min( +0.0,-0.0) == -0.0 */
-        if (d1 == 0.0
-                && ((Double.doubleToLongBits(d1) | Double.doubleToLongBits(d2)) & 0x8000000000000000l) != 0)
+        if (d1 == 0.0 &&
+                ((Double.doubleToLongBits(d1) | Double.doubleToLongBits(d2)) & 0x8000000000000000l) != 0)
             return 0.0 * (-1.0);
         return d1;
     }
@@ -550,8 +1396,8 @@ public final class StrictMath {
         if (f1 != f2)
             return Float.NaN;
         /* min( +0.0,-0.0) == -0.0 */
-        if (f1 == 0.0f
-                && ((Float.floatToIntBits(f1) | Float.floatToIntBits(f2)) & 0x80000000) != 0)
+        if (f1 == 0.0f &&
+                ((Float.floatToIntBits(f1) | Float.floatToIntBits(f2)) & 0x80000000) != 0)
             return 0.0f * (-1.0f);
         return f1;
     }
@@ -706,7 +1552,7 @@ public final class StrictMath {
      *            the value whose signum has to be computed.
      * @return the value of the signum function.
      */
-    public static double signum(double d){
+    public static double signum(double d) {
         return Math.signum(d);
     }
 
@@ -729,10 +1575,12 @@ public final class StrictMath {
      *            the value whose signum has to be computed.
      * @return the value of the signum function.
      */
-    public static float signum(float f){
+    public static float signum(float f) {
         return Math.signum(f);
     }
 
+    private static final double shuge = 1.0e307;
+
     /**
      * Returns the closest double approximation of the hyperbolic sine of the
      * argument.
@@ -746,11 +1594,57 @@ public final class StrictMath {
      * <li>{@code sinh(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose hyperbolic sine has to be computed.
      * @return the hyperbolic sine of the argument.
      */
-    public static native double sinh(double d);
+    public static double sinh(double x) {
+        double t, w, h;
+        int ix, jx;
+        final long bits = Double.doubleToRawLongBits(x);
+
+        jx = (int) (bits >>> 32);
+        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 > 1.00000000000000000000e+00) {
+                    return x;/* ieee_sinh(tiny) = tiny with inexact */
+                }
+            t = expm1(Math.abs(x));
+            if (ix < 0x3ff00000)
+                return h * (2.0 * t - t * t / (t + 1.00000000000000000000e+00));
+            return h * (t + t / (t + 1.00000000000000000000e+00));
+        }
+
+        /* |x| in [22, ieee_log(maxdouble)] return 0.5*ieee_exp(|x|) */
+        if (ix < 0x40862E42) {
+            return h * exp(Math.abs(x));
+        }
+
+        /* |x| in [log(maxdouble), overflowthresold] */
+        final long lx = ((ONEBITS >>> 29) + ((int) bits)) & 0x00000000ffffffffL;
+        // lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
+        if (ix < 0x408633CE || (ix == 0x408633ce) && (lx <= 0x8fb9f87dL)) {
+            w = exp(0.5 * Math.abs(x));
+            t = h * w;
+            return t * w;
+        }
+
+        /* |x| > overflowthresold, ieee_sinh(x) overflow */
+        return x * shuge;
+    }
 
     /**
      * Returns the closest double approximation of the sine of the argument.
@@ -816,11 +1710,47 @@ public final class StrictMath {
      * <li>{@code tanh(NaN) = NaN}</li>
      * </ul>
      *
-     * @param d
+     * @param x
      *            the value whose hyperbolic tangent has to be computed.
      * @return the hyperbolic tangent of the argument
      */
-    public static native double tanh(double d);
+    public static double tanh(double x) {
+        double t, z;
+        int jx, ix;
+
+        final long bits = Double.doubleToRawLongBits(x);
+        /* High word of |x|. */
+        jx = (int) (bits >>> 32);
+        ix = jx & 0x7fffffff;
+
+        /* x is INF or NaN */
+        if (ix >= 0x7ff00000) {
+            if (jx >= 0) {
+                return 1.00000000000000000000e+00 / x + 1.00000000000000000000e+00; /* ieee_tanh(+-inf)=+-1 */
+            } else {
+                return 1.00000000000000000000e+00 / x - 1.00000000000000000000e+00; /* ieee_tanh(NaN) = NaN */
+            }
+        }
+
+        /* |x| < 22 */
+        if (ix < 0x40360000) { /* |x|<22 */
+            if (ix < 0x3c800000) { /* |x|<2**-55 */
+                return x * (1.00000000000000000000e+00 + x);/* ieee_tanh(small) = small */
+            }
+
+            if (ix >= 0x3ff00000) { /* |x|>=1 */
+                t = Math.expm1(2.0 * Math.abs(x));
+                z = 1.00000000000000000000e+00 - 2.0 / (t + 2.0);
+            } else {
+                t = Math.expm1(-2.0 * Math.abs(x));
+                z = -t / (t + 2.0);
+            }
+            /* |x| > 22, return +-1 */
+        } else {
+            z = 1.00000000000000000000e+00 - TINY; /* raised inexact flag */
+        }
+        return (jx >= 0) ? z : -z;
+    }
 
     /**
      * Returns the measure in degrees of the supplied radian angle. The result
@@ -922,6 +1852,7 @@ public final class StrictMath {
     /**
      * Returns a double with the given magnitude and the sign of {@code sign}.
      * If {@code sign} is NaN, the sign of the result is positive.
+     *
      * @since 1.6
      */
     public static double copySign(double magnitude, double sign) {
@@ -932,13 +1863,15 @@ public final class StrictMath {
         // (Tested on a Nexus One.)
         long magnitudeBits = Double.doubleToRawLongBits(magnitude);
         long signBits = Double.doubleToRawLongBits((sign != sign) ? 1.0 : sign);
-        magnitudeBits = (magnitudeBits & ~Double.SIGN_MASK) | (signBits & Double.SIGN_MASK);
+        magnitudeBits = (magnitudeBits & ~Double.SIGN_MASK)
+                | (signBits & Double.SIGN_MASK);
         return Double.longBitsToDouble(magnitudeBits);
     }
 
     /**
-     * Returns a float with the given magnitude and the sign of {@code sign}.
-     * If {@code sign} is NaN, the sign of the result is positive.
+     * Returns a float with the given magnitude and the sign of {@code sign}. If
+     * {@code sign} is NaN, the sign of the result is positive.
+     *
      * @since 1.6
      */
     public static float copySign(float magnitude, float sign) {
@@ -949,12 +1882,14 @@ public final class StrictMath {
         // (Tested on a Nexus One.)
         int magnitudeBits = Float.floatToRawIntBits(magnitude);
         int signBits = Float.floatToRawIntBits((sign != sign) ? 1.0f : sign);
-        magnitudeBits = (magnitudeBits & ~Float.SIGN_MASK) | (signBits & Float.SIGN_MASK);
+        magnitudeBits = (magnitudeBits & ~Float.SIGN_MASK)
+                | (signBits & Float.SIGN_MASK);
         return Float.intBitsToFloat(magnitudeBits);
     }
 
     /**
      * Returns the exponent of float {@code f}.
+     *
      * @since 1.6
      */
     public static int getExponent(float f) {
@@ -963,14 +1898,17 @@ public final class StrictMath {
 
     /**
      * Returns the exponent of double {@code d}.
+     *
      * @since 1.6
      */
-    public static int getExponent(double d){
+    public static int getExponent(double d) {
         return Math.getExponent(d);
     }
 
     /**
-     * Returns the next double after {@code start} in the given {@code direction}.
+     * Returns the next double after {@code start} in the given
+     * {@code direction}.
+     *
      * @since 1.6
      */
     public static double nextAfter(double start, double direction) {
@@ -981,7 +1919,9 @@ public final class StrictMath {
     }
 
     /**
-     * Returns the next float after {@code start} in the given {@code direction}.
+     * Returns the next float after {@code start} in the given {@code direction}
+     * .
+     *
      * @since 1.6
      */
     public static float nextAfter(float start, double direction) {
@@ -990,6 +1930,7 @@ public final class StrictMath {
 
     /**
      * Returns the next double larger than {@code d}.
+     *
      * @since 1.6
      */
     public static double nextUp(double d) {
@@ -998,6 +1939,7 @@ public final class StrictMath {
 
     /**
      * Returns the next float larger than {@code f}.
+     *
      * @since 1.6
      */
     public static float nextUp(float f) {
@@ -1006,6 +1948,7 @@ public final class StrictMath {
 
     /**
      * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded.
+     *
      * @since 1.6
      */
     public static double scalb(double d, int scaleFactor) {
@@ -1049,12 +1992,10 @@ public final class StrictMath {
         } else {
             if (Math.abs(d) >= Double.MIN_NORMAL) {
                 // common situation
-                result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS)
-                        | (bits & Double.MANTISSA_MASK);
+                result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS) | (bits & Double.MANTISSA_MASK);
             } else {
                 // origin d is sub-normal, change mantissa to normal style
-                result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS)
-                        | ((bits << (subNormalFactor + 1)) & Double.MANTISSA_MASK);
+                result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS) | ((bits << (subNormalFactor + 1)) & Double.MANTISSA_MASK);
             }
         }
         return Double.longBitsToDouble(result | sign);
@@ -1062,6 +2003,7 @@ public final class StrictMath {
 
     /**
      * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded.
+     *
      * @since 1.6
      */
     public static float scalb(float d, int scaleFactor) {
@@ -1073,8 +2015,7 @@ public final class StrictMath {
         int factor = ((bits & Float.EXPONENT_MASK) >> Float.MANTISSA_BITS)
                 - Float.EXPONENT_BIAS + scaleFactor;
         // calculates the factor of sub-normal values
-        int subNormalFactor = Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK)
-                - Float.EXPONENT_BITS;
+        int subNormalFactor = Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK) - Float.EXPONENT_BITS;
         if (subNormalFactor < 0) {
             // not sub-normal values
             subNormalFactor = 0;
@@ -1105,8 +2046,9 @@ public final class StrictMath {
                         | (bits & Float.MANTISSA_MASK);
             } else {
                 // origin d is sub-normal, change mantissa to normal style
-                result = ((factor + Float.EXPONENT_BIAS) << Float.MANTISSA_BITS)
-                        | ((bits << (subNormalFactor + 1)) & Float.MANTISSA_MASK);
+                result = ((factor + Float.EXPONENT_BIAS)
+                        << Float.MANTISSA_BITS) | (
+                        (bits << (subNormalFactor + 1)) & Float.MANTISSA_MASK);
             }
         }
         return Float.intBitsToFloat(result | sign);
@@ -1120,10 +2062,10 @@ public final class StrictMath {
         }
         // change it to positive
         int absDigits = -digits;
-        if (Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK) <= (32 - absDigits)) {
+        if (Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK)
+                <= (32 - absDigits)) {
             // some bits will remain after shifting, calculates its carry
-            if ((((bits >> (absDigits - 1)) & 0x1) == 0)
-                    || Integer.numberOfTrailingZeros(bits) == (absDigits - 1)) {
+            if ((((bits >> (absDigits - 1)) & 0x1) == 0) || Integer.numberOfTrailingZeros(bits) == (absDigits - 1)) {
                 return bits >> absDigits;
             }
             return ((bits >> absDigits) + 1);
@@ -1139,10 +2081,10 @@ public final class StrictMath {
         }
         // change it to positive
         long absDigits = -digits;
-        if (Long.numberOfLeadingZeros(bits & ~Double.SIGN_MASK) <= (64 - absDigits)) {
+        if (Long.numberOfLeadingZeros(bits & ~Double.SIGN_MASK)
+                <= (64 - absDigits)) {
             // some bits will remain after shifting, calculates its carry
-            if ((((bits >> (absDigits - 1)) & 0x1) == 0)
-                    || Long.numberOfTrailingZeros(bits) == (absDigits - 1)) {
+            if ((((bits >> (absDigits - 1)) & 0x1) == 0) || Long.numberOfTrailingZeros(bits) == (absDigits - 1)) {
                 return bits >> absDigits;
             }
             return ((bits >> absDigits) + 1);
diff --git a/luni/src/main/native/java_lang_StrictMath.cpp b/luni/src/main/native/java_lang_StrictMath.cpp
index cfe375e..e8c6dfb 100644
--- a/luni/src/main/native/java_lang_StrictMath.cpp
+++ b/luni/src/main/native/java_lang_StrictMath.cpp
@@ -34,26 +34,6 @@ static jdouble StrictMath_tan(JNIEnv*, jclass, jdouble a) {
     return ieee_tan(a);
 }
 
-static jdouble StrictMath_asin(JNIEnv*, jclass, jdouble a) {
-    return ieee_asin(a);
-}
-
-static jdouble StrictMath_acos(JNIEnv*, jclass, jdouble a) {
-    return ieee_acos(a);
-}
-
-static jdouble StrictMath_atan(JNIEnv*, jclass, jdouble a) {
-    return ieee_atan(a);
-}
-
-static jdouble StrictMath_exp(JNIEnv*, jclass, jdouble a) {
-    return ieee_exp(a);
-}
-
-static jdouble StrictMath_log(JNIEnv*, jclass, jdouble a) {
-    return ieee_log(a);
-}
-
 static jdouble StrictMath_sqrt(JNIEnv*, jclass, jdouble a) {
     return ieee_sqrt(a);
 }
@@ -74,75 +54,30 @@ static jdouble StrictMath_rint(JNIEnv*, jclass, jdouble a) {
     return ieee_rint(a);
 }
 
-static jdouble StrictMath_atan2(JNIEnv*, jclass, jdouble a, jdouble b) {
-    return ieee_atan2(a, b);
-}
-
 static jdouble StrictMath_pow(JNIEnv*, jclass, jdouble a, jdouble b) {
     return ieee_pow(a,b);
 }
 
-static jdouble StrictMath_sinh(JNIEnv*, jclass, jdouble a) {
-    return ieee_sinh(a);
-}
-
-static jdouble StrictMath_tanh(JNIEnv*, jclass, jdouble a) {
-    return ieee_tanh(a);
-}
-
-static jdouble StrictMath_cosh(JNIEnv*, jclass, jdouble a) {
-    return ieee_cosh(a);
-}
-
-static jdouble StrictMath_log10(JNIEnv*, jclass, jdouble a) {
-    return ieee_log10(a);
-}
-
-static jdouble StrictMath_cbrt(JNIEnv*, jclass, jdouble a) {
-    return ieee_cbrt(a);
-}
-
-static jdouble StrictMath_expm1(JNIEnv*, jclass, jdouble a) {
-    return ieee_expm1(a);
-}
-
 static jdouble StrictMath_hypot(JNIEnv*, jclass, jdouble a, jdouble b) {
     return ieee_hypot(a, b);
 }
 
-static jdouble StrictMath_log1p(JNIEnv*, jclass, jdouble a) {
-    return ieee_log1p(a);
-}
-
 static jdouble StrictMath_nextafter(JNIEnv*, jclass, jdouble a, jdouble b) {
     return ieee_nextafter(a, b);
 }
 
 static JNINativeMethod gMethods[] = {
     NATIVE_METHOD(StrictMath, IEEEremainder, "!(DD)D"),
-    NATIVE_METHOD(StrictMath, acos, "!(D)D"),
-    NATIVE_METHOD(StrictMath, asin, "!(D)D"),
-    NATIVE_METHOD(StrictMath, atan, "!(D)D"),
-    NATIVE_METHOD(StrictMath, atan2, "!(DD)D"),
-    NATIVE_METHOD(StrictMath, cbrt, "!(D)D"),
     NATIVE_METHOD(StrictMath, ceil, "!(D)D"),
     NATIVE_METHOD(StrictMath, cos, "!(D)D"),
-    NATIVE_METHOD(StrictMath, cosh, "!(D)D"),
-    NATIVE_METHOD(StrictMath, exp, "!(D)D"),
-    NATIVE_METHOD(StrictMath, expm1, "!(D)D"),
     NATIVE_METHOD(StrictMath, floor, "!(D)D"),
     NATIVE_METHOD(StrictMath, hypot, "!(DD)D"),
-    NATIVE_METHOD(StrictMath, log, "!(D)D"),
-    NATIVE_METHOD(StrictMath, log10, "!(D)D"),
-    NATIVE_METHOD(StrictMath, log1p, "!(D)D"),
     NATIVE_METHOD(StrictMath, nextafter, "!(DD)D"),
     NATIVE_METHOD(StrictMath, pow, "!(DD)D"),
     NATIVE_METHOD(StrictMath, rint, "!(D)D"),
     NATIVE_METHOD(StrictMath, sin, "!(D)D"),
-    NATIVE_METHOD(StrictMath, sinh, "!(D)D"),
     NATIVE_METHOD(StrictMath, sqrt, "!(D)D"),
     NATIVE_METHOD(StrictMath, tan, "!(D)D"),
-    NATIVE_METHOD(StrictMath, tanh, "!(D)D"),
 };
 void register_java_lang_StrictMath(JNIEnv* env) {
     jniRegisterNativeMethods(env, "java/lang/StrictMath", gMethods, NELEM(gMethods));