path: root/luni/src/main/java/java/lang/StrictMath.java

/*
 *  Licensed to the Apache Software Foundation (ASF) under one or more
 *  contributor license agreements.  See the NOTICE file distributed with
 *  this work for additional information regarding copyright ownership.
 *  The ASF licenses this file to You under the Apache License, Version 2.0
 *  (the "License"); you may not use this file except in compliance with
 *  the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  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;

/**
 * Class StrictMath provides basic math constants and operations such as
 * trigonometric functions, hyperbolic functions, exponential, logarithms, etc.
 * <p>
 * In contrast to class {@link Math}, the methods in this class return exactly
 * the same results on all platforms. Algorithms based on these methods thus
 * behave the same (e.g. regarding numerical convergence) on all platforms,
 * complying with the slogan "write once, run everywhere". On the other side,
 * the implementation of class StrictMath may be less efficient than that of
 * class Math, as class StrictMath cannot utilize platform specific features
 * such as an extended precision math co-processors.
 * <p>
 * The methods in this class are specified using the "Freely Distributable Math
 * Library" (fdlibm), version 5.3.
 * <p>
 * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a>
 */
public final class StrictMath {
    /**
     * The double value closest to e, the base of the natural logarithm.
     */
    public static final double E = Math.E;

    /**
     * The double value closest to pi, the ratio of a circle's circumference to
     * its diameter.
     */
    public static final double PI = Math.PI;

    /**
     * Prevents this class from being instantiated.
     */
    private StrictMath() {
    }

    /**
     * Returns the absolute value of the argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code abs(-0.0) = +0.0}</li>
     * <li>{@code abs(+infinity) = +infinity}</li>
     * <li>{@code abs(-infinity) = +infinity}</li>
     * <li>{@code abs(NaN) = NaN}</li>
     * </ul>
     */
    public static double abs(double d) {
        return Math.abs(d);
    }

    /**
     * Returns the absolute value of the argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code abs(-0.0) = +0.0}</li>
     * <li>{@code abs(+infinity) = +infinity}</li>
     * <li>{@code abs(-infinity) = +infinity}</li>
     * <li>{@code abs(NaN) = NaN}</li>
     * </ul>
     */
    public static float abs(float f) {
        return Math.abs(f);
    }

    /**
     * Returns the absolute value of the argument.
     * <p>
     * If the argument is {@code Integer.MIN_VALUE}, {@code Integer.MIN_VALUE}
     * is returned.
     */
    public static int abs(int i) {
        return Math.abs(i);
    }

    /**
     * Returns the absolute value of the argument.
     * <p>
     * If the argument is {@code Long.MIN_VALUE}, {@code Long.MIN_VALUE} is
     * returned.
     */
    public static long abs(long l) {
        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]}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code acos((anything > 1) = NaN}</li>
     * <li>{@code acos((anything < -1) = NaN}</li>
     * <li>{@code acos(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value to compute arc cosine of.
     * @return the arc cosine of the argument.
     */
    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
     * within the range {@code [-pi/2..pi/2]}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code asin((anything > 1)) = NaN}</li>
     * <li>{@code asin((anything < -1)) = NaN}</li>
     * <li>{@code asin(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose arc sine has to be computed.
     * @return the arc sine of the argument.
     */
    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
     * argument within the range {@code [-pi/2..pi/2]}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code atan(+0.0) = +0.0}</li>
     * <li>{@code atan(-0.0) = -0.0}</li>
     * <li>{@code atan(+infinity) = +pi/2}</li>
     * <li>{@code atan(-infinity) = -pi/2}</li>
     * <li>{@code atan(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose arc tangent has to be computed.
     * @return the arc tangent of the argument.
     */
    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
     * {@code y/x} within the range {@code [-pi..pi]}. This is the angle of the
     * polar representation of the rectangular coordinates (x,y).
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code atan2((anything), NaN ) = NaN;}</li>
     * <li>{@code atan2(NaN , (anything) ) = NaN;}</li>
     * <li>{@code atan2(+0.0, +(anything but NaN)) = +0.0}</li>
     * <li>{@code atan2(-0.0, +(anything but NaN)) = -0.0}</li>
     * <li>{@code atan2(+0.0, -(anything but NaN)) = +pi}</li>
     * <li>{@code atan2(-0.0, -(anything but NaN)) = -pi}</li>
     * <li>{@code atan2(+(anything but 0 and NaN), 0) = +pi/2}</li>
     * <li>{@code atan2(-(anything but 0 and NaN), 0) = -pi/2}</li>
     * <li>{@code atan2(+(anything but infinity and NaN), +infinity)} {@code =}
     * {@code +0.0}</li>
     * <li>{@code atan2(-(anything but infinity and NaN), +infinity)} {@code =}
     * {@code -0.0}</li>
     * <li>{@code atan2(+(anything but infinity and NaN), -infinity) = +pi}</li>
     * <li>{@code atan2(-(anything but infinity and NaN), -infinity) = -pi}</li>
     * <li>{@code atan2(+infinity, +infinity ) = +pi/4}</li>
     * <li>{@code atan2(-infinity, +infinity ) = -pi/4}</li>
     * <li>{@code atan2(+infinity, -infinity ) = +3pi/4}</li>
     * <li>{@code atan2(-infinity, -infinity ) = -3pi/4}</li>
     * <li>{@code atan2(+infinity, (anything but,0, NaN, and infinity))}
     * {@code =} {@code +pi/2}</li>
     * <li>{@code atan2(-infinity, (anything but,0, NaN, and infinity))}
     * {@code =} {@code -pi/2}</li>
     * </ul>
     *
     * @param y
     *            the numerator of the value whose atan has to be computed.
     * @param x
     *            the denominator of the value whose atan has to be computed.
     * @return the arc tangent of {@code y/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
     * argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code cbrt(+0.0) = +0.0}</li>
     * <li>{@code cbrt(-0.0) = -0.0}</li>
     * <li>{@code cbrt(+infinity) = +infinity}</li>
     * <li>{@code cbrt(-infinity) = -infinity}</li>
     * <li>{@code cbrt(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose cube root has to be computed.
     * @return the cube root of the argument.
     */
    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
     * infinity) integer value greater than or equal to the argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code ceil(+0.0) = +0.0}</li>
     * <li>{@code ceil(-0.0) = -0.0}</li>
     * <li>{@code ceil((anything in range (-1,0)) = -0.0}</li>
     * <li>{@code ceil(+infinity) = +infinity}</li>
     * <li>{@code ceil(-infinity) = -infinity}</li>
     * <li>{@code ceil(NaN) = NaN}</li>
     * </ul>
     */
    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
     * argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code cosh(+infinity) = +infinity}</li>
     * <li>{@code cosh(-infinity) = +infinity}</li>
     * <li>{@code cosh(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose hyperbolic cosine has to be computed.
     * @return the hyperbolic cosine of the argument.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code cos(+infinity) = NaN}</li>
     * <li>{@code cos(-infinity) = NaN}</li>
     * <li>{@code cos(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the angle whose cosine has to be computed, in radians.
     * @return the cosine of the argument.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code exp(+infinity) = +infinity}</li>
     * <li>{@code exp(-infinity) = +0.0}</li>
     * <li>{@code exp(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose exponential has to be computed.
     * @return the exponential of the argument.
     */
    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>
     * {@code d}</sup>{@code - 1}. If the argument is very close to 0, it is
     * much more accurate to use {@code expm1(d)+1} than {@code exp(d)} (due to
     * cancellation of significant digits).
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code expm1(+0.0) = +0.0}</li>
     * <li>{@code expm1(-0.0) = -0.0}</li>
     * <li>{@code expm1(+infinity) = +infinity}</li>
     * <li>{@code expm1(-infinity) = -1.0}</li>
     * <li>{@code expm1(NaN) = NaN}</li>
     * </ul>
     *
     * @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.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code floor(+0.0) = +0.0}</li>
     * <li>{@code floor(-0.0) = -0.0}</li>
     * <li>{@code floor(+infinity) = +infinity}</li>
     * <li>{@code floor(-infinity) = -infinity}</li>
     * <li>{@code floor(NaN) = NaN}</li>
     * </ul>
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code hypot(+infinity, (anything including NaN)) = +infinity}</li>
     * <li>{@code hypot(-infinity, (anything including NaN)) = +infinity}</li>
     * <li>{@code hypot((anything including NaN), +infinity) = +infinity}</li>
     * <li>{@code hypot((anything including NaN), -infinity) = +infinity}</li>
     * <li>{@code hypot(NaN, NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            a double number.
     * @param y
     *            a double number.
     * @return the {@code sqrt(}<i>{@code x}</i><sup>{@code 2}</sup>{@code +}
     *         <i> {@code y}</i><sup>{@code 2}</sup>{@code )} value of the
     *         arguments.
     */
    public static native double hypot(double x, double y);

    /**
     * Returns the remainder of dividing {@code x} by {@code y} using the IEEE
     * 754 rules. The result is {@code x-round(x/p)*p} where {@code round(x/p)}
     * is the nearest integer (rounded to even), but without numerical
     * cancellation problems.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code IEEEremainder((anything), 0) = NaN}</li>
     * <li>{@code IEEEremainder(+infinity, (anything)) = NaN}</li>
     * <li>{@code IEEEremainder(-infinity, (anything)) = NaN}</li>
     * <li>{@code IEEEremainder(NaN, (anything)) = NaN}</li>
     * <li>{@code IEEEremainder((anything), NaN) = NaN}</li>
     * <li>{@code IEEEremainder(x, +infinity) = x } where x is anything but
     * +/-infinity</li>
     * <li>{@code IEEEremainder(x, -infinity) = x } where x is anything but
     * +/-infinity</li>
     * </ul>
     *
     * @param x
     *            the numerator of the operation.
     * @param y
     *            the denominator of the operation.
     * @return the IEEE754 floating point reminder of of {@code x/y}.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code log(+0.0) = -infinity}</li>
     * <li>{@code log(-0.0) = -infinity}</li>
     * <li>{@code log((anything < 0) = NaN}</li>
     * <li>{@code log(+infinity) = +infinity}</li>
     * <li>{@code log(-infinity) = NaN}</li>
     * <li>{@code log(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose log has to be computed.
     * @return the natural logarithm of the argument.
     */
    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
     * argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code log10(+0.0) = -infinity}</li>
     * <li>{@code log10(-0.0) = -infinity}</li>
     * <li>{@code log10((anything < 0) = NaN}</li>
     * <li>{@code log10(+infinity) = +infinity}</li>
     * <li>{@code log10(-infinity) = NaN}</li>
     * <li>{@code log10(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose base 10 log has to be computed.
     * @return the the base 10 logarithm of x
     */
    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
     * sum of the argument and 1. If the argument is very close to 0, it is much
     * more accurate to use {@code log1p(d)} than {@code log(1.0+d)} (due to
     * numerical cancellation).
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code log1p(+0.0) = +0.0}</li>
     * <li>{@code log1p(-0.0) = -0.0}</li>
     * <li>{@code log1p((anything < 1)) = NaN}</li>
     * <li>{@code log1p(-1.0) = -infinity}</li>
     * <li>{@code log1p(+infinity) = +infinity}</li>
     * <li>{@code log1p(-infinity) = NaN}</li>
     * <li>{@code log1p(NaN) = NaN}</li>
     * </ul>
     *
     * @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 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
     * arguments.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code max(NaN, (anything)) = NaN}</li>
     * <li>{@code max((anything), NaN) = NaN}</li>
     * <li>{@code max(+0.0, -0.0) = +0.0}</li>
     * <li>{@code max(-0.0, +0.0) = +0.0}</li>
     * </ul>
     */
    public static double max(double d1, double d2) {
        if (d1 > d2)
            return d1;
        if (d1 < d2)
            return d2;
        /* if either arg is NaN, return NaN */
        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)
            return 0.0;
        return d1;
    }

    /**
     * Returns the most positive (closest to positive infinity) of the two
     * arguments.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code max(NaN, (anything)) = NaN}</li>
     * <li>{@code max((anything), NaN) = NaN}</li>
     * <li>{@code max(+0.0, -0.0) = +0.0}</li>
     * <li>{@code max(-0.0, +0.0) = +0.0}</li>
     * </ul>
     */
    public static float max(float f1, float f2) {
        if (f1 > f2)
            return f1;
        if (f1 < f2)
            return f2;
        /* if either arg is NaN, return NaN */
        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)
            return 0.0f;
        return f1;
    }

    /**
     * Returns the most positive (closest to positive infinity) of the two
     * arguments.
     */
    public static int max(int i1, int i2) {
        return Math.max(i1, i2);
    }

    /**
     * Returns the most positive (closest to positive infinity) of the two
     * arguments.
     */
    public static long max(long l1, long l2) {
        return l1 > l2 ? l1 : l2;
    }

    /**
     * Returns the most negative (closest to negative infinity) of the two
     * arguments.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code min(NaN, (anything)) = NaN}</li>
     * <li>{@code min((anything), NaN) = NaN}</li>
     * <li>{@code min(+0.0, -0.0) = -0.0}</li>
     * <li>{@code min(-0.0, +0.0) = -0.0}</li>
     * </ul>
     */
    public static double min(double d1, double d2) {
        if (d1 > d2)
            return d2;
        if (d1 < d2)
            return d1;
        /* if either arg is NaN, return NaN */
        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)
            return 0.0 * (-1.0);
        return d1;
    }

    /**
     * Returns the most negative (closest to negative infinity) of the two
     * arguments.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code min(NaN, (anything)) = NaN}</li>
     * <li>{@code min((anything), NaN) = NaN}</li>
     * <li>{@code min(+0.0, -0.0) = -0.0}</li>
     * <li>{@code min(-0.0, +0.0) = -0.0}</li>
     * </ul>
     */
    public static float min(float f1, float f2) {
        if (f1 > f2)
            return f2;
        if (f1 < f2)
            return f1;
        /* if either arg is NaN, return NaN */
        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)
            return 0.0f * (-1.0f);
        return f1;
    }

    /**
     * Returns the most negative (closest to negative infinity) of the two
     * arguments.
     */
    public static int min(int i1, int i2) {
        return Math.min(i1, i2);
    }

    /**
     * Returns the most negative (closest to negative infinity) of the two
     * arguments.
     */
    public static long min(long l1, long l2) {
        return l1 < l2 ? l1 : l2;
    }

    /**
     * Returns the closest double approximation of the result of raising
     * {@code x} to the power of {@code y}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code pow((anything), +0.0) = 1.0}</li>
     * <li>{@code pow((anything), -0.0) = 1.0}</li>
     * <li>{@code pow(x, 1.0) = x}</li>
     * <li>{@code pow((anything), NaN) = NaN}</li>
     * <li>{@code pow(NaN, (anything except 0)) = NaN}</li>
     * <li>{@code pow(+/-(|x| > 1), +infinity) = +infinity}</li>
     * <li>{@code pow(+/-(|x| > 1), -infinity) = +0.0}</li>
     * <li>{@code pow(+/-(|x| < 1), +infinity) = +0.0}</li>
     * <li>{@code pow(+/-(|x| < 1), -infinity) = +infinity}</li>
     * <li>{@code pow(+/-1.0 , +infinity) = NaN}</li>
     * <li>{@code pow(+/-1.0 , -infinity) = NaN}</li>
     * <li>{@code pow(+0.0, (+anything except 0, NaN)) = +0.0}</li>
     * <li>{@code pow(-0.0, (+anything except 0, NaN, odd integer)) = +0.0}</li>
     * <li>{@code pow(+0.0, (-anything except 0, NaN)) = +infinity}</li>
     * <li>{@code pow(-0.0, (-anything except 0, NAN, odd integer))} {@code =}
     * {@code +infinity}</li>
     * <li>{@code pow(-0.0, (odd integer)) = -pow( +0 , (odd integer) )}</li>
     * <li>{@code pow(+infinity, (+anything except 0, NaN)) = +infinity}</li>
     * <li>{@code pow(+infinity, (-anything except 0, NaN)) = +0.0}</li>
     * <li>{@code pow(-infinity, (anything)) = -pow(0, (-anything))}</li>
     * <li>{@code pow((-anything), (integer))} {@code =}
     * {@code pow(-1,(integer))*pow(+anything,integer)}</li>
     * <li>{@code pow((-anything except 0 and infinity), (non-integer))}
     * {@code =} {@code NAN}</li>
     * </ul>
     *
     * @param x
     *            the base of the operation.
     * @param y
     *            the exponent of the operation.
     * @return {@code x} to the power of {@code y}.
     */
    public static native double pow(double x, double y);

    /**
     * Returns a pseudo-random number between 0.0 (inclusive) and 1.0
     * (exclusive).
     *
     * @return a pseudo-random number.
     */
    public static double random() {
        return Math.random();
    }

    /**
     * Returns the double conversion of the result of rounding the argument to
     * an integer. Tie breaks are rounded towards even.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code rint(+0.0) = +0.0}</li>
     * <li>{@code rint(-0.0) = -0.0}</li>
     * <li>{@code rint(+infinity) = +infinity}</li>
     * <li>{@code rint(-infinity) = -infinity}</li>
     * <li>{@code rint(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the value to be rounded.
     * @return the closest integer to the argument (as a double).
     */
    public static native double rint(double d);

    /**
     * Returns the result of rounding the argument to an integer. The result is
     * equivalent to {@code (long) Math.floor(d+0.5)}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code round(+0.0) = +0.0}</li>
     * <li>{@code round(-0.0) = +0.0}</li>
     * <li>{@code round((anything > Long.MAX_VALUE) = Long.MAX_VALUE}</li>
     * <li>{@code round((anything < Long.MIN_VALUE) = Long.MIN_VALUE}</li>
     * <li>{@code round(+infinity) = Long.MAX_VALUE}</li>
     * <li>{@code round(-infinity) = Long.MIN_VALUE}</li>
     * <li>{@code round(NaN) = +0.0}</li>
     * </ul>
     *
     * @param d
     *            the value to be rounded.
     * @return the closest integer to the argument.
     */
    public static long round(double d) {
        return Math.round(d);
    }

    /**
     * Returns the result of rounding the argument to an integer. The result is
     * equivalent to {@code (int) Math.floor(f+0.5)}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code round(+0.0) = +0.0}</li>
     * <li>{@code round(-0.0) = +0.0}</li>
     * <li>{@code round((anything > Integer.MAX_VALUE) = Integer.MAX_VALUE}</li>
     * <li>{@code round((anything < Integer.MIN_VALUE) = Integer.MIN_VALUE}</li>
     * <li>{@code round(+infinity) = Integer.MAX_VALUE}</li>
     * <li>{@code round(-infinity) = Integer.MIN_VALUE}</li>
     * <li>{@code round(NaN) = +0.0}</li>
     * </ul>
     *
     * @param f
     *            the value to be rounded.
     * @return the closest integer to the argument.
     */
    public static int round(float f) {
        return Math.round(f);
    }

    /**
     * Returns the signum function of the argument. If the argument is less than
     * zero, it returns -1.0. If the argument is greater than zero, 1.0 is
     * returned. If the argument is either positive or negative zero, the
     * argument is returned as result.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code signum(+0.0) = +0.0}</li>
     * <li>{@code signum(-0.0) = -0.0}</li>
     * <li>{@code signum(+infinity) = +1.0}</li>
     * <li>{@code signum(-infinity) = -1.0}</li>
     * <li>{@code signum(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the value whose signum has to be computed.
     * @return the value of the signum function.
     */
    public static double signum(double d) {
        return Math.signum(d);
    }

    /**
     * Returns the signum function of the argument. If the argument is less than
     * zero, it returns -1.0. If the argument is greater than zero, 1.0 is
     * returned. If the argument is either positive or negative zero, the
     * argument is returned as result.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code signum(+0.0) = +0.0}</li>
     * <li>{@code signum(-0.0) = -0.0}</li>
     * <li>{@code signum(+infinity) = +1.0}</li>
     * <li>{@code signum(-infinity) = -1.0}</li>
     * <li>{@code signum(NaN) = NaN}</li>
     * </ul>
     *
     * @param f
     *            the value whose signum has to be computed.
     * @return the value of the signum function.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code sinh(+0.0) = +0.0}</li>
     * <li>{@code sinh(-0.0) = -0.0}</li>
     * <li>{@code sinh(+infinity) = +infinity}</li>
     * <li>{@code sinh(-infinity) = -infinity}</li>
     * <li>{@code sinh(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose hyperbolic sine has to be computed.
     * @return the hyperbolic sine of the argument.
     */
    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.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code sin(+0.0) = +0.0}</li>
     * <li>{@code sin(-0.0) = -0.0}</li>
     * <li>{@code sin(+infinity) = NaN}</li>
     * <li>{@code sin(-infinity) = NaN}</li>
     * <li>{@code sin(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the angle whose sin has to be computed, in radians.
     * @return the sine of the argument.
     */
    public static native double sin(double d);

    /**
     * Returns the closest double approximation of the square root of the
     * argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code sqrt(+0.0) = +0.0}</li>
     * <li>{@code sqrt(-0.0) = -0.0}</li>
     * <li>{@code sqrt( (anything < 0) ) = NaN}</li>
     * <li>{@code sqrt(+infinity) = +infinity}</li>
     * <li>{@code sqrt(NaN) = NaN}</li>
     * </ul>
     */
    public static native double sqrt(double d);

    /**
     * Returns the closest double approximation of the tangent of the argument.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code tan(+0.0) = +0.0}</li>
     * <li>{@code tan(-0.0) = -0.0}</li>
     * <li>{@code tan(+infinity) = NaN}</li>
     * <li>{@code tan(-infinity) = NaN}</li>
     * <li>{@code tan(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the angle whose tangent has to be computed, in radians.
     * @return the tangent of the argument.
     */
    public static native double tan(double d);

    /**
     * Returns the closest double approximation of the hyperbolic tangent of the
     * argument. The absolute value is always less than 1.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code tanh(+0.0) = +0.0}</li>
     * <li>{@code tanh(-0.0) = -0.0}</li>
     * <li>{@code tanh(+infinity) = +1.0}</li>
     * <li>{@code tanh(-infinity) = -1.0}</li>
     * <li>{@code tanh(NaN) = NaN}</li>
     * </ul>
     *
     * @param x
     *            the value whose hyperbolic tangent has to be computed.
     * @return the hyperbolic tangent of the argument
     */
    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
     * is {@code angrad * 180 / pi}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code toDegrees(+0.0) = +0.0}</li>
     * <li>{@code toDegrees(-0.0) = -0.0}</li>
     * <li>{@code toDegrees(+infinity) = +infinity}</li>
     * <li>{@code toDegrees(-infinity) = -infinity}</li>
     * <li>{@code toDegrees(NaN) = NaN}</li>
     * </ul>
     *
     * @param angrad
     *            an angle in radians.
     * @return the degree measure of the angle.
     */
    public static double toDegrees(double angrad) {
        return Math.toDegrees(angrad);
    }

    /**
     * Returns the measure in radians of the supplied degree angle. The result
     * is {@code angdeg / 180 * pi}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code toRadians(+0.0) = +0.0}</li>
     * <li>{@code toRadians(-0.0) = -0.0}</li>
     * <li>{@code toRadians(+infinity) = +infinity}</li>
     * <li>{@code toRadians(-infinity) = -infinity}</li>
     * <li>{@code toRadians(NaN) = NaN}</li>
     * </ul>
     *
     * @param angdeg
     *            an angle in degrees.
     * @return the radian measure of the angle.
     */
    public static double toRadians(double angdeg) {
        return Math.toRadians(angdeg);
    }

    /**
     * Returns the argument's ulp (unit in the last place). The size of a ulp of
     * a double value is the positive distance between this value and the double
     * value next larger in magnitude. For non-NaN {@code x},
     * {@code ulp(-x) == ulp(x)}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code ulp(+0.0) = Double.MIN_VALUE}</li>
     * <li>{@code ulp(-0.0) = Double.MIN_VALUE}</li>
     * <li>{@code ulp(+infinity) = infinity}</li>
     * <li>{@code ulp(-infinity) = infinity}</li>
     * <li>{@code ulp(NaN) = NaN}</li>
     * </ul>
     *
     * @param d
     *            the floating-point value to compute ulp of.
     * @return the size of a ulp of the argument.
     */
    public static double ulp(double d) {
        // special cases
        if (Double.isInfinite(d)) {
            return Double.POSITIVE_INFINITY;
        } else if (d == Double.MAX_VALUE || d == -Double.MAX_VALUE) {
            return pow(2, 971);
        }
        d = Math.abs(d);
        return nextafter(d, Double.MAX_VALUE) - d;
    }

    /**
     * Returns the argument's ulp (unit in the last place). The size of a ulp of
     * a float value is the positive distance between this value and the float
     * value next larger in magnitude. For non-NaN {@code x},
     * {@code ulp(-x) == ulp(x)}.
     * <p>
     * Special cases:
     * <ul>
     * <li>{@code ulp(+0.0) = Float.MIN_VALUE}</li>
     * <li>{@code ulp(-0.0) = Float.MIN_VALUE}</li>
     * <li>{@code ulp(+infinity) = infinity}</li>
     * <li>{@code ulp(-infinity) = infinity}</li>
     * <li>{@code ulp(NaN) = NaN}</li>
     * </ul>
     *
     * @param f
     *            the floating-point value to compute ulp of.
     * @return the size of a ulp of the argument.
     */
    public static float ulp(float f) {
        return Math.ulp(f);
    }

    private static native double nextafter(double x, double y);

    /**
     * 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) {
        // We manually inline Double.isNaN here because the JIT can't do it yet.
        // With Double.isNaN: 236.3ns
        // With manual inline: 141.2ns
        // With no check (i.e. Math's behavior): 110.0ns
        // (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);
        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.
     *
     * @since 1.6
     */
    public static float copySign(float magnitude, float sign) {
        // We manually inline Float.isNaN here because the JIT can't do it yet.
        // With Float.isNaN: 214.7ns
        // With manual inline: 112.3ns
        // With no check (i.e. Math's behavior): 93.1ns
        // (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);
        return Float.intBitsToFloat(magnitudeBits);
    }

    /**
     * Returns the exponent of float {@code f}.
     *
     * @since 1.6
     */
    public static int getExponent(float f) {
        return Math.getExponent(f);
    }

    /**
     * Returns the exponent of double {@code d}.
     *
     * @since 1.6
     */
    public static int getExponent(double d) {
        return Math.getExponent(d);
    }

    /**
     * Returns the next double after {@code start} in the given
     * {@code direction}.
     *
     * @since 1.6
     */
    public static double nextAfter(double start, double direction) {
        if (start == 0 && direction == 0) {
            return direction;
        }
        return nextafter(start, direction);
    }

    /**
     * Returns the next float after {@code start} in the given {@code direction}
     * .
     *
     * @since 1.6
     */
    public static float nextAfter(float start, double direction) {
        return Math.nextAfter(start, direction);
    }

    /**
     * Returns the next double larger than {@code d}.
     *
     * @since 1.6
     */
    public static double nextUp(double d) {
        return Math.nextUp(d);
    }

    /**
     * Returns the next float larger than {@code f}.
     *
     * @since 1.6
     */
    public static float nextUp(float f) {
        return Math.nextUp(f);
    }

    /**
     * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded.
     *
     * @since 1.6
     */
    public static double scalb(double d, int scaleFactor) {
        if (Double.isNaN(d) || Double.isInfinite(d) || d == 0) {
            return d;
        }
        // change double to long for calculation
        long bits = Double.doubleToLongBits(d);
        // the sign of the results must be the same of given d
        long sign = bits & Double.SIGN_MASK;
        // calculates the factor of the result
        long factor = (int) ((bits & Double.EXPONENT_MASK) >> Double.MANTISSA_BITS)
                - Double.EXPONENT_BIAS + scaleFactor;

        // calculates the factor of sub-normal values
        int subNormalFactor = Long.numberOfLeadingZeros(bits & ~Double.SIGN_MASK)
                - Double.EXPONENT_BITS;
        if (subNormalFactor < 0) {
            // not sub-normal values
            subNormalFactor = 0;
        }
        if (Math.abs(d) < Double.MIN_NORMAL) {
            factor = factor - subNormalFactor;
        }
        if (factor > Double.MAX_EXPONENT) {
            return (d > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
        }

        long result;
        // if result is a sub-normal
        if (factor < -Double.EXPONENT_BIAS) {
            // the number of digits that shifts
            long digits = factor + Double.EXPONENT_BIAS + subNormalFactor;
            if (Math.abs(d) < Double.MIN_NORMAL) {
                // origin d is already sub-normal
                result = shiftLongBits(bits & Double.MANTISSA_MASK, digits);
            } else {
                // origin d is not sub-normal, change mantissa to sub-normal
                result = shiftLongBits(bits & Double.MANTISSA_MASK | 0x0010000000000000L, digits - 1);
            }
        } else {
            if (Math.abs(d) >= Double.MIN_NORMAL) {
                // common situation
                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);
            }
        }
        return Double.longBitsToDouble(result | sign);
    }

    /**
     * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded.
     *
     * @since 1.6
     */
    public static float scalb(float d, int scaleFactor) {
        if (Float.isNaN(d) || Float.isInfinite(d) || d == 0) {
            return d;
        }
        int bits = Float.floatToIntBits(d);
        int sign = bits & Float.SIGN_MASK;
        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;
        if (subNormalFactor < 0) {
            // not sub-normal values
            subNormalFactor = 0;
        }
        if (Math.abs(d) < Float.MIN_NORMAL) {
            factor = factor - subNormalFactor;
        }
        if (factor > Float.MAX_EXPONENT) {
            return (d > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY);
        }

        int result;
        // if result is a sub-normal
        if (factor < -Float.EXPONENT_BIAS) {
            // the number of digits that shifts
            int digits = factor + Float.EXPONENT_BIAS + subNormalFactor;
            if (Math.abs(d) < Float.MIN_NORMAL) {
                // origin d is already sub-normal
                result = shiftIntBits(bits & Float.MANTISSA_MASK, digits);
            } else {
                // origin d is not sub-normal, change mantissa to sub-normal
                result = shiftIntBits(bits & Float.MANTISSA_MASK | 0x00800000, digits - 1);
            }
        } else {
            if (Math.abs(d) >= Float.MIN_NORMAL) {
                // common situation
                result = ((factor + Float.EXPONENT_BIAS) << Float.MANTISSA_BITS)
                        | (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);
            }
        }
        return Float.intBitsToFloat(result | sign);
    }

    // Shifts integer bits as float, if the digits is positive, left-shift; if
    // not, shift to right and calculate its carry.
    private static int shiftIntBits(int bits, int digits) {
        if (digits > 0) {
            return bits << digits;
        }
        // change it to positive
        int absDigits = -digits;
        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)) {
                return bits >> absDigits;
            }
            return ((bits >> absDigits) + 1);
        }
        return 0;
    }

    // Shifts long bits as double, if the digits is positive, left-shift; if
    // not, shift to right and calculate its carry.
    private static long shiftLongBits(long bits, long digits) {
        if (digits > 0) {
            return bits << digits;
        }
        // change it to positive
        long absDigits = -digits;
        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)) {
                return bits >> absDigits;
            }
            return ((bits >> absDigits) + 1);
        }
        return 0;
    }
}