Implements some math functions for faster performance

Changes cosh, exp, expm1, log, log10, sinh, cbrt,
asin, acos, atan, tanh, log1p
These functions are now implemnted in java to
save the overhead of a jni call

Change-Id: I4d3b70dd1693378af5b04cbbe0fecd95d97c1cff

2 files changed, 1006 insertions, 129 deletions

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));