/*
 * -------------------------------------------------------------------------
 *        $Id: JMath.java,v 1.2 2009/02/16 05:19:12 jfrijters Exp $
 * -------------------------------------------------------------------------
 *        Copyright (c) 1999 Visual Numerics Inc. All Rights Reserved.
 *
 * Permission to use, copy, modify, and distribute this software is freely
 * granted by Visual Numerics, Inc., provided that the copyright notice
 * above and the following warranty disclaimer are preserved in human
 * readable form.
 *
 * Because this software is licenses free of charge, it is provided
 * "AS IS", with NO WARRANTY.  TO THE EXTENT PERMITTED BY LAW, VNI
 * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
 * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 * VNI WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
 * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
 * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
 *
 *
 * This Java code is based on C code in the package fdlibm,
 * which can be obtained from www.netlib.org.
 * The original fdlibm C code contains the following notice.
 *
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 *
 *--------------------------------------------------------------------------
 */

package ikvm.internal;
import java.util.Random;


/*
 *        Pure Java implementation of the standard java.lang.Math class.
 *        This Java code is based on C code in the package fdlibm,
 *        which can be obtained from www.netlib.org.
 *
 *        @author Sun Microsystems (original C code in fdlibm)
 *        @author        John F. Brophy (translated from C to Java)
 */
@ikvm.lang.Internal
public final class JMath {
    static public final double PI = 0x1.921fb54442d18p1; /* 3.14159265358979323846 */
    static public final double E = 2.7182818284590452354;
    static private Random random;

    /**
     *        Returns the absolute value of its argument.
     *        @param        x        The argument, an integer.
     *        @return        Returns |x|.
     */
    strictfp public static int abs(int x) {
        return ((x < 0) ? (-x) : x);
    }

    /**
     *        Returns the absolute value of its argument.
     *        @param        x        The argument, a long.
     *        @return        Returns |x|.
     */
    strictfp public static long abs(long x) {
        return ((x < 0L) ? (-x) : x);
    }

    /**
     *        Returns the absolute value of its argument.
     *        @param        x        The argument, a float.
     *        @return        Returns |x|.
     */
    strictfp public static float abs(float x) {
        return ((x <= 0.0f) ? (0.0f - x) : x);
    }

    /**
     *        Returns the absolute value of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns |x|.
     */
    strictfp public static double abs(double x) {
        return ((x <= 0.0) ? (0.0 - x) : x);
    }

    /**
     *        Returns the smaller of its two arguments.
     *        @param        x        The first argument, an integer.
     *        @param        y        The second argument, an integer.
     *        @return        Returns the smaller of x and y.
     */
    strictfp public static int min(int x, int y) {
        return ((x < y) ? x : y);
    }

    /**
     *        Returns the smaller of its two arguments.
     *        @param        x        The first argument, a long.
     *        @param        y        The second argument, a long.
     *        @return        Returns the smaller of x and y.
     */
    strictfp public static long min(long x, long y) {
        return ((x < y) ? x : y);
    }

    /**
     *        Returns the smaller of its two arguments.
     *        @param        x        The first argument, a float.
     *        @param        y        The second argument, a float.
     *        @return        Returns the smaller of x and y.
     *                        This function considers -0.0f to
     *                        be less than 0.0f.
     */
    strictfp public static float min(float x, float y) {
        if (Float.isNaN(x)) {
            return x;
        }
        float ans = ((x <= y) ? x : y);
        if ((ans == 0.0f) && (Float.floatToIntBits(y) == 0x80000000)) {
            ans = y;
        }
        return ans;
    }

    /**
     *        Returns the smaller of its two arguments.
     *        @param        x        The first argument, a double.
     *        @param        y        The second argument, a double.
     *        @return        Returns the smaller of x and y.
     *                        This function considers -0.0 to
     *                        be less than 0.0.
     */
    strictfp public static double min(double x, double y) {
        if (Double.isNaN(x)) {
            return x;
        }
        double ans = ((x <= y) ? x : y);
        if ((x == 0.0) && (y == 0.0) &&
                (Double.doubleToLongBits(y) == 0x8000000000000000L)) {
            ans = y;
        }
        return ans;
    }

    /**
     *        Returns the larger of its two arguments.
     *        @param        x        The first argument, an integer.
     *        @param        y        The second argument, an integer.
     *        @return        Returns the larger of x and y.
     */
    strictfp public static int max(int x, int y) {
        return ((x > y) ? x : y);
    }

    /**
     *        Returns the larger of its two arguments.
     *        @param        x        The first argument, a long.
     *        @param        y        The second argument, a long.
     *        @return        Returns the larger of x and y.
     */
    strictfp public static long max(long x, long y) {
        return ((x > y) ? x : y);
    }

    /**
     *        Returns the larger of its two arguments.
     *        @param        x        The first argument, a float.
     *        @param        y        The second argument, a float.
     *        @return        Returns the larger of x and y.
     *                        This function considers -0.0f to
     *                        be less than 0.0f.
     */
    strictfp public static float max(float x, float y) {
        if (Float.isNaN(x)) {
            return x;
        }
        float ans = ((x >= y) ? x : y);
        if ((ans == 0.0f) && (Float.floatToIntBits(x) == 0x80000000)) {
            ans = y;
        }
        return ans;
    }

    /**
     *        Returns the larger of its two arguments.
     *        @param        x        The first argument, a double.
     *        @param        y        The second argument, a double.
     *        @return        Returns the larger of x and y.
     *                        This function considers -0.0 to
     *                        be less than 0.0.
     */
    strictfp public static double max(double x, double y) {
        if (Double.isNaN(x)) {
            return x;
        }
        double ans = ((x >= y) ? x : y);
        if ((x == 0.0) && (y == 0.0) &&
                (Double.doubleToLongBits(x) == 0x8000000000000000L)) {
            ans = y;
        }
        return ans;
    }

    /**
     *        Returns the integer closest to the arguments.
     *        @param        x        The argument, a float.
     *        @return        Returns the integer closest to x.
     */
    strictfp public static int round(float x) {
        return (int) floor(x + 0.5f);
    }

    /**
     *        Returns the long closest to the arguments.
     *        @param        x        The argument, a double.
     *        @return        Returns the long closest to x.
     */
    strictfp public static long round(double x) {
        return (long) floor(x + 0.5);
    }

    /**
     *        Returns the random number.
     *        @return        Returns a random number from a uniform distribution.
     */
    synchronized strictfp public static double random() {
        if (random == null) {
            random = new Random();
        }
        return random.nextDouble();
    }

    /*
     *        This following code is derived from fdlibm, which contained
     *        the following notice.
     * ====================================================
     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     *
     * Developed at SunSoft, a Sun Microsystems, Inc. business.
     * Permission to use, copy, modify, and distribute this
     * software is freely granted, provided that this notice
     * is preserved.
     * ====================================================
     *
     * Constants:
     * The hexadecimal values are the intended ones for the following
     * constants. The decimal values may be used, provided that the
     * compiler will convert from decimal to binary accurately enough
     * to produce the hexadecimal values shown.
     */
    static private final double huge = 1.0e+300;
    static private final double tiny = 1.0e-300;

    /**
     *        Returns the value of its argument rounded toward
     *        positive infinity to an integral value.
     *        @param        x        The argument, a double.
     *        @return        Returns the smallest double, not less than x,
     *                        that is an integral value.
     */
    static public double ceil(double x) {
        int exp;
        int sign;
        long ix;

        if (x == 0) {
            return x;
        }

        ix = Double.doubleToLongBits(x);
        sign = (int) ((ix >> 63) & 1);
        exp = ((int) (ix >> 52) & 0x7ff) - 0x3ff;

        if (exp < 0) {
            if (x < 0.0) {
                return NEGATIVE_ZERO;
            } else if (x == 0.0) {
                return x;
            } else {
                return 1.0;
            }
        } else if (exp < 53) {
            long mask = (0x000fffffffffffffL >>> exp);
            if ((mask & ix) == 0) {
                return x; // x is integral
            }
            if (x > 0.0) {
                ix += (0x0010000000000000L >> exp);
            }
            ix = ix & (~mask);
        } else if (exp == 1024) { // infinity
            return x;
        }

        return Double.longBitsToDouble(ix);
    }

    /**
     *        Returns the value of its argument rounded toward
     *        negative infinity to an integral value.
     *        @param        x        The argument, a double.
     *        @return        Returns the smallest double, not greater than x,
     *                        that is an integral value.
     */
    static public double floor(double x) {
        int exp;
        int sign;
        long ix;

        if (x == 0) {
            return x;
        }

        ix = Double.doubleToLongBits(x);
        sign = (int) ((ix >> 63) & 1);
        exp = ((int) (ix >> 52) & 0x7ff) - 0x3ff;

        if (exp < 0) {
            if (x < 0.0) {
                return -1.0;
            } else if (x == 0.0) {
                return x;
            } else {
                return 0.0;
            }
        } else if (exp < 53) {
            long mask = (0x000fffffffffffffL >>> exp);
            if ((mask & ix) == 0) {
                return x; // x is integral
            }
            if (x < 0.0) {
                ix += (0x0010000000000000L >> exp);
            }
            ix = ix & (~mask);
        } else if (exp == 1024) { // infinity
            return x;
        }

        return Double.longBitsToDouble(ix);
    }

    static private final double[] TWO52 = {
            0x1.0p52, /*  4.50359962737049600000e+15 */
            -0x1.0p52 /* -4.50359962737049600000e+15 */
        };
    static private final double NEGATIVE_ZERO = -0x0.0p0;

    /**
     *        Returns the value of its argument rounded toward
     *        the closest integral value.
     *        @param        x        The argument, a double.
     *        @return        Returns the double closest to x
     *                        that is an integral value.
     */
    static public double rint(double x) {
        int exp;
        int sign;
        long ix;
        double w;

        if (x == 0) {
            return x;
        }

        ix = Double.doubleToLongBits(x);
        sign = (int) ((ix >> 63) & 1);
        exp = ((int) (ix >> 52) & 0x7ff) - 0x3ff;

        if (exp < 0) {
            if (x < -0.5) {
                return -1.0;
            } else if (x > 0.5) {
                return 1.0;
            } else if (sign == 0) {
                return 0.0;
            } else {
                return NEGATIVE_ZERO;
            }
        } else if (exp < 53) {
            long mask = (0x000fffffffffffffL >>> exp);
            if ((mask & ix) == 0) {
                return x; // x is integral
            }
        } else if (exp == 1024) { // infinity
            return x;
        }

        x = Double.longBitsToDouble(ix);
        w = TWO52[sign] + x;
        return w - TWO52[sign];
    }

    /**
     *        Returns  x REM p =  x - [x/p]*p as if in infinite
     *         precise arithmetic, where [x/p] is the (infinite bit)
     *        integer nearest x/p (in half way case choose the even one).
     *        @param        x        The dividend.
     *        @param        y        The divisor.
     *        @return        The remainder computed according to the IEEE 754 standard.
     */
    static public double IEEEremainder(double x, double p) {
        int hx;
        int hp;
        int sx; // unsigned
        int lx; // unsigned
        int lp; // unsigned
        double p_half;

        hx = __HI(x); /* high word of x */
        lx = __LO(x); /* low  word of x */
        hp = __HI(p); /* high word of p */
        lp = __LO(p); /* low  word of p */
        sx = hx & 0x80000000;
        hp &= 0x7fffffff;
        hx &= 0x7fffffff;

        /* purge off exception values */
        if ((hp | lp) == 0) {
            return (x * p) / (x * p); /* p = 0 */
        }

        if ((hx >= 0x7ff00000) || /* x not finite */
                ((hp >= 0x7ff00000) && /* p is NaN */
                (((hp - 0x7ff00000) | lp) != 0))) {
            return (x * p) / (x * p);
        }

        if (hp <= 0x7fdfffff) {
            x = x % (p + p); /* now x < 2p */
        }

        if (((hx - hp) | (lx - lp)) == 0) {
            return zero * x;
        }

        x = abs(x);
        p = abs(p);
        if (hp < 0x00200000) {
            if ((x + x) > p) {
                x -= p;
                if ((x + x) >= p) {
                    x -= p;
                }
            }
        } else {
            p_half = 0.5 * p;
            if (x > p_half) {
                x -= p;
                if (x >= p_half) {
                    x -= p;
                }
            }
        }
        lx = __HI(x);
        lx ^= sx;
        return setHI(x, lx);
    }

    /* sqrt(x)
     * Return correctly rounded sqrt.
     *   ------------------------------------------
     *   |  Use the hardware sqrt if you have one |
     *   ------------------------------------------
     * Method:
     *   Bit by bit method using integer arithmetic. (Slow, but portable)
     *   1. Normalization
     *  Scale x to y in [1,4) with even powers of 2:
     *  find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
     *      sqrt(x) = 2^k * sqrt(y)
     *   2. Bit by bit computation
     *  Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
     *       i                                                   0
     *                             i+1         2
     *      s  = 2*q , and  y  =  2   * ( y - q  ).     (1)
     *       i      i            i                 i
     *
     *  To compute q    from q , one checks whether
     *              i+1       i
     *
     *                -(i+1) 2
     *          (q + 2      ) <= y.         (2)
     *                i
     *                         -(i+1)
     *  If (2) is false, then q   = q ; otherwise q   = q  + 2      .
     *                         i+1   i             i+1   i
     *
     *  With some algebric manipulation, it is not difficult to see
     *  that (2) is equivalent to
     *                             -(i+1)
     *          s  +  2       <= y          (3)
     *           i                i
     *
     *  The advantage of (3) is that s  and y  can be computed by
     *                    i      i
     *  the following recurrence formula:
     *      if (3) is false
     *
     *      s     =  s  ,   y    = y   ;            (4)
     *       i+1      i      i+1    i
     *
     *      otherwise,
     *                     -i                      -(i+1)
     *      s     =  s  + 2  ,  y    = y  -  s  - 2         (5)
     *       i+1      i          i+1    i     i
     *
     *  One may easily use induction to prove (4) and (5).
     *  Note. Since the left hand side of (3) contain only i+2 bits,
     *        it does not necessary to do a full (53-bit) comparison
     *        in (3).
     *   3. Final rounding
     *  After generating the 53 bits result, we compute one more bit.
     *  Together with the remainder, we can decide whether the
     *  result is exact, bigger than 1/2ulp, or less than 1/2ulp
     *  (it will never equal to 1/2ulp).
     *  The rounding mode can be detected by checking whether
     *  huge + tiny is equal to huge, and whether huge - tiny is
     *  equal to huge for some floating point number "huge" and "tiny".
     *
     *
     * Special cases:
     *  sqrt(+-0) = +-0     ... exact
     *  sqrt(inf) = inf
     *  sqrt(-ve) = NaN     ... with invalid signal
     *  sqrt(NaN) = NaN     ... with invalid signal for signaling NaN
     */

    /**
     *        Returns the square root of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns the square root of x.
     */
    static public double sqrt(double x) {
        long ix = Double.doubleToLongBits(x);

        /* take care of Inf and NaN */
        if ((ix & 0x7ff0000000000000L) == 0x7ff0000000000000L) {

            /* sqrt(NaN)=NaN, sqrt(+inf)=+inf  sqrt(-inf)=sNaN */
            return (x * x) + x;
        }

        /* take care of zero */
        if (x < 0.0) {
            return Double.NaN;
        } else if (x == 0.0) {
            return x; /* sqrt(+-0) = +-0 */
        }

        /* normalize x */
        long m = (ix >> 52);
        ix &= 0x000fffffffffffffL;

        /* add implicit bit, if not sub-normal */
        if (m != 0) {
            ix |= 0x0010000000000000L;
        }

        m -= 1023L; /* unbias exponent */
        if ((m & 1) != 0) { /* odd m, double x to make it even */
            ix += ix;
        }
        m >>= 1; /* m = [m/2] */
        m += 1023L;

        /* generate sqrt(x) bit by bit */
        ix += ix;
        long q = 0L; /* q = sqrt(x) */
        long s = 0L;
        long r = 0x0020000000000000L; /* r = moving bit from right to left */

        while (r != 0) {
            long t = s + r;
            if (t <= ix) {
                s = t + r;
                ix -= t;
                q += r;
            }
            ix += ix;
            r >>= 1;
        }

        /* round */
        if (ix != 0) {
            q += (q & 1L);
        }

        /* assemble result */
        ix = (m << 52) | (0x000fffffffffffffL & (q >> 1));
        return Double.longBitsToDouble(ix);
    }

    static private final double[] halF = { 0.5, -0.5 };
    static private final double twom1000 = 0x1.0p-1000; /* 2**-1000=9.33263618503218878990e-302 */
    static private final double o_threshold = 0x1.62e42fefa39efp9; /* 7.09782712893383973096e+02 */
    static private final double u_threshold = -0x1.74910d52d3051p9; /* -7.45133219101941108420e+02 */
    static private final double[] ln2HI = {
            0x1.62e42feep-1, /*  6.93147180369123816490e-01 */
            -0x1.62e42feep-1
        }; /* -6.93147180369123816490e-01 */
    static private final double[] ln2LO = {
            0x1.a39ef35793c76p-33, /*  1.90821492927058770002e-10 */
            -0x1.a39ef35793c76p-33
        }; /* -1.90821492927058770002e-10 */
    static private final double invln2 = 0x1.71547652b82fep0; /* 1.44269504088896338700e+00 */
    static private final double P1 = 0x1.555555555553ep-3; /*  1.66666666666666019037e-01 */
    static private final double P2 = -0x1.6c16c16bebd93p-9; /* -2.77777777770155933842e-03 */
    static private final double P3 = 0x1.1566aaf25de2cp-14; /*  6.61375632143793436117e-05 */
    static private final double P4 = -0x1.bbd41c5d26bf1p-20; /* -1.65339022054652515390e-06 */
    static private final double P5 = 0x1.6376972bea4dp-25; /*  4.13813679705723846039e-08 */

    /* exp(x)
     * Returns the exponential of x.
     *
     * Method
     *   1. Argument reduction:
     *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
     *        Given x, find r and integer k such that
     *
     *               x = k*ln2 + r,  |r| <= 0.5*ln2.
     *
     *      Here r will be represented as r = hi-lo for better
     *         accuracy.
     *
     *   2. Approximation of exp(r) by a special rational function on
     *        the interval [0,0.34658]:
     *        Write
     *            R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
     *      We use a special Reme algorithm on [0,0.34658] to generate
     *         a polynomial of degree 5 to approximate R. The maximum error
     *          of this polynomial approximation is bounded by 2**-59. In
     *        other words,
     *              R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
     *         (where z=r*r, and the values of P1 to P5 are listed below)
     *        and
     *            |                  5          |     -59
     *             | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
     *            |                             |
     *        The computation of exp(r) thus becomes
     *                             2*r
     *                exp(r) = 1 + -------
     *                              R - r
     *                          r*R1(r)
     *                       = 1 + r + ----------- (for better accuracy)
     *                                  2 - R1(r)
     *         where
     *                                   2       4             10
     *                R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
     *
     *   3. Scale back to obtain exp(x):
     *        From step 1, we have
     *           exp(x) = 2^k * exp(r)
     *
     * Special cases:
     *        exp(INF) is INF, exp(NaN) is NaN;
     *        exp(-INF) is 0, and
     *        for finite argument, only exp(0)=1 is exact.
     *
     * Accuracy:
     *        according to an error analysis, the error is always less than
     *        1 ulp (unit in the last place).
     *
     * Misc. info.
     *        For IEEE double
     *            if x >  7.09782712893383973096e+02 then exp(x) overflow
     *            if x < -7.45133219101941108420e+02 then exp(x) underflow
     */

    /**
     *        Returns the exponential of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns e to the power x.
     */
    static public double exp(double x) {
        double y;
        double hi = 0;
        double lo = 0;
        double c;
        double t;
        int k = 0;
        int xsb;
        int hx;

        hx = __HI(x); /* high word of x */
        xsb = (hx >>> 31) & 1; /* sign bit of x */
        hx &= 0x7fffffff; /* high word of |x| */

        /* filter out non-finite argument */
        if (hx >= 0x40862E42) { /* if |x|>=709.78... */
            if (hx >= 0x7ff00000) {
                if (((hx & 0xfffff) | __LO(x)) != 0) {
                    return x + x; /* NaN */
                } else {
                    return ((xsb == 0) ? x : 0.0); /* exp(+-inf)={inf,0} */
                }
            }
            if (x > o_threshold) {
                return huge * huge; /* overflow */
            }
            if (x < u_threshold) {
                return twom1000 * twom1000; /* underflow */
            }
        }

        /* argument reduction */
        if (hx > 0x3fd62e42) { /* if  |x| > 0.5 ln2 */
            if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
                hi = x - ln2HI[xsb];
                lo = ln2LO[xsb];
                k = 1 - xsb - xsb;
            } else {
                k = (int) ((invln2 * x) + halF[xsb]);
                t = k;
                hi = x - (t * ln2HI[0]); /* t*ln2HI is exact here */
                lo = t * ln2LO[0];
            }
            x = hi - lo;
        } else if (hx < 0x3e300000) { /* when |x|<2**-28 */
            if ((huge + x) > one) {
                return one + x; /* trigger inexact */
            }
        } else {
            k = 0;
        }

        /* x is now in primary range */
        t = x * x;
        c = x - (t * (P1 + (t * (P2 + (t * (P3 + (t * (P4 + (t * P5)))))))));
        if (k == 0) {
            return one - (((x * c) / (c - 2.0)) - x);
        } else {
            y = one - ((lo - ((x * c) / (2.0 - c))) - hi);
        }

        long iy = Double.doubleToLongBits(y);
        if (k >= -1021) {
            iy += ((long) k << 52);
        } else {
            iy += ((k + 1000L) << 52);
        }
        return Double.longBitsToDouble(iy);
    }

    static private final double ln2_hi = 0x1.62e42feep-1; /* 6.93147180369123816490e-01 */
    static private final double ln2_lo = 0x1.a39ef35793c76p-33; /* 1.90821492927058770002e-10 */
    static private final double Lg1 = 0x1.5555555555593p-1; /* 6.666666666666735130e-01 */
    static private final double Lg2 = 0x1.999999997fa04p-2; /* 3.999999999940941908e-01 */
    static private final double Lg3 = 0x1.2492494229359p-2; /* 2.857142874366239149e-01 */
    static private final double Lg4 = 0x1.c71c51d8e78afp-3; /* 2.222219843214978396e-01 */
    static private final double Lg5 = 0x1.7466496cb03dep-3; /* 1.818357216161805012e-01 */
    static private final double Lg6 = 0x1.39a09d078c69fp-3; /* 1.531383769920937332e-01 */
    static private final double Lg7 = 0x1.2f112df3e5244p-3; /* 1.479819860511658591e-01 */

    /*
     * Return the logrithm of x
     *
     * Method :
     *   1. Argument Reduction: find k and f such that
     *                        x = 2^k * (1+f),
     *           where  sqrt(2)/2 < 1+f < sqrt(2) .
     *
     *   2. Approximation of log(1+f).
     *        Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
     *                 = 2s + 2/3 s**3 + 2/5 s**5 + .....,
     *                      = 2s + s*R
     *      We use a special Reme algorithm on [0,0.1716] to generate
     *         a polynomial of degree 14 to approximate R The maximum error
     *        of this polynomial approximation is bounded by 2**-58.45. In
     *        other words,
     *                        2      4      6      8      10      12      14
     *            R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
     *          (the values of Lg1 to Lg7 are listed in the program)
     *        and
     *            |      2          14          |     -58.45
     *            | Lg1*s +...+Lg7*s    -  R(z) | <= 2
     *            |                             |
     *        Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
     *        In order to guarantee error in log below 1ulp, we compute log
     *        by
     *                log(1+f) = f - s*(f - R)        (if f is not too large)
     *                log(1+f) = f - (hfsq - s*(hfsq+R)).        (better accuracy)
     *
     *        3. Finally,  log(x) = k*ln2 + log(1+f).
     *                            = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
     *           Here ln2 is split into two floating point number:
     *                        ln2_hi + ln2_lo,
     *           where n*ln2_hi is always exact for |n| < 2000.
     *
     * Special cases:
     *        log(x) is NaN with signal if x < 0 (including -INF) ;
     *        log(+INF) is +INF; log(0) is -INF with signal;
     *        log(NaN) is that NaN with no signal.
     *
     * Accuracy:
     *        according to an error analysis, the error is always less than
     *        1 ulp (unit in the last place).
     */

    /**
     *        Returns the natural logarithm of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns the natural (base e) logarithm of x.
     */
    static public double log(double x) {
        double hfsq;
        double f;
        double s;
        double z;
        double R;
        double w;
        double t1;
        double t2;
        double dk;
        int k;
        int hx;
        int i;
        int j;
        int lx;

        hx = __HI(x); /* high word of x */
        lx = __LO(x); /* low  word of x */

        k = 0;
        if (hx < 0x00100000) { /* x < 2**-1022  */
            if (((hx & 0x7fffffff) | lx) == 0) {
                return -two54 / zero; /* log(+-0)=-inf */
            }
            if (hx < 0) {
                return (x - x) / zero; /* log(-#) = NaN */
            }
            k -= 54;
            x *= two54; /* subnormal number, scale up x */
            hx = __HI(x); /* high word of x */
        }
        if (hx >= 0x7ff00000) {
            return x + x;
        }
        k += ((hx >> 20) - 1023);
        hx &= 0x000fffff;
        i = (hx + 0x95f64) & 0x100000;
        x = setHI(x, hx | (i ^ 0x3ff00000)); /* normalize x or x/2 */
        k += (i >> 20);
        f = x - 1.0;
        if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */
            if (f == zero) {
                if (k == 0) {
                    return zero;
                } else {
                    dk = (double) k;
                }
                return (dk * ln2_hi) + (dk * ln2_lo);
            }
            R = f * f * (0.5 - (0.33333333333333333 * f));
            if (k == 0) {
                return f - R;
            } else {
                dk = (double) k;
                return (dk * ln2_hi) - ((R - (dk * ln2_lo)) - f);
            }
        }
        s = f / (2.0 + f);
        dk = (double) k;
        z = s * s;
        i = hx - 0x6147a;
        w = z * z;
        j = 0x6b851 - hx;
        t1 = w * (Lg2 + (w * (Lg4 + (w * Lg6))));
        t2 = z * (Lg1 + (w * (Lg3 + (w * (Lg5 + (w * Lg7))))));
        i |= j;
        R = t2 + t1;
        if (i > 0) {
            hfsq = 0.5 * f * f;
            if (k == 0) {
                return f - (hfsq - (s * (hfsq + R)));
            } else {
                return (dk * ln2_hi) -
                ((hfsq - ((s * (hfsq + R)) + (dk * ln2_lo))) - f);
            }
        } else {
            if (k == 0) {
                return f - (s * (f - R));
            } else {
                return (dk * ln2_hi) - (((s * (f - R)) - (dk * ln2_lo)) - f);
            }
        }
    }

    /**
     *        Returns the sine of its argument.
     *        @param        x        The argument, a double, assumed to be in radians.
     *        @return        Returns the sine of x.
     */
    static public double sin(double x) {
        double z = 0.0;
        int n;
        int ix = __HI(x);

        ix &= 0x7fffffff; /* |x| ~< pi/4 */

        if (ix <= 0x3fe921fb) {
            return __kernel_sin(x, z, 0);
        } else if (ix >= 0x7ff00000) {

            /* sin(Inf or NaN) is NaN */
            return x - x;
        } else {

            /* argument reduction needed */
            double[] y = new double[2];
            n = __ieee754_rem_pio2(x, y);
            switch (n & 3) {
            case 0:
                return __kernel_sin(y[0], y[1], 1);
            case 1:
                return __kernel_cos(y[0], y[1]);
            case 2:
                return -__kernel_sin(y[0], y[1], 1);
            default:
                return -__kernel_cos(y[0], y[1]);
            }
        }
    }

    static private final double S1 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */
    static private final double S2 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */
    static private final double S3 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */
    static private final double S4 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */
    static private final double S5 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */
    static private final double S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */

    /*
     * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
     * Input x is assumed to be bounded by ~pi/4 in magnitude.
     * Input y is the tail of x. * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
     *
     * Algorithm
     *        1. Since sin(-x) = -sin(x), we need only to consider positive x.
     *        2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
     *        3. sin(x) is approximated by a polynomial of degree 13 on
     *           [0,pi/4]
     *                                   3            13
     *                   sin(x) ~ x + S1*x + ... + S6*x
     *           where
     *
     *         |sin(x)         2     4     6     8     10     12  |     -58
     *         |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x  +S6*x   )| <= 2
     *         |  x                                                    |
     *
     *        4. sin(x+y) = sin(x) + sin'(x')*y
     *                    ~ sin(x) + (1-x*x/2)*y
     *           For better accuracy, let
     *                     3      2      2      2      2
     *                r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
     *           then                   3    2
     *                sin(x) = x + (S1*x + (x *(r-y/2)+y))
     */
    static double __kernel_sin(double x, double y, int iy) {
        double z;
        double r;
        double v;
        int ix;
        ix = __HI(x) & 0x7fffffff; /* high word of x */
        if (ix < 0x3e400000) { /* |x| < 2**-27 */
            if ((int) x == 0) {
                return x; /* generate inexact */
            }
        }
        z = x * x;
        v = z * x;
        r = S2 + (z * (S3 + (z * (S4 + (z * (S5 + (z * S6)))))));
        if (iy == 0) {
            return x + (v * (S1 + (z * r)));
        } else {
            return x - (((z * ((half * y) - (v * r))) - y) - (v * S1));
        }
    }

    /**
     *        Returns the cosine of its argument.
     *        @param        x        The argument, a double, assumed to be in radians.
     *        @return        Returns the cosine of x.
     */
    static public double cos(double x) {
        double z = 0.0;
        int n;
        int ix;

        /* High word of x. */
        ix = __HI(x);

        /* |x| ~< pi/4 */
        ix &= 0x7fffffff;
        if (ix <= 0x3fe921fb) {
            return __kernel_cos(x, z);

            /* cos(Inf or NaN) is NaN */
        } else if (ix >= 0x7ff00000) {
            return x - x;

            /* argument reduction needed */
        } else {
            double[] y = new double[2];
            n = __ieee754_rem_pio2(x, y);
            switch (n & 3) {
            case 0:
                return __kernel_cos(y[0], y[1]);
            case 1:
                return -__kernel_sin(y[0], y[1], 1);
            case 2:
                return -__kernel_cos(y[0], y[1]);
            default:
                return __kernel_sin(y[0], y[1], 1);
            }
        }
    }

    static private final double one = 0x1.0p0; /*  1.00000000000000000000e+00 */
    static private final double C1 = 0x1.555555555554cp-5; /*  4.16666666666666019037e-02 */
    static private final double C2 = -0x1.6c16c16c15177p-10; /* -1.38888888888741095749e-03 */
    static private final double C3 = 0x1.a01a019cb159p-16; /*  2.48015872894767294178e-05 */
    static private final double C4 = -0x1.27e4f809c52adp-22; /* -2.75573143513906633035e-07 */
    static private final double C5 = 0x1.1ee9ebdb4b1c4p-29; /*  2.08757232129817482790e-09 */
    static private final double C6 = -0x1.8fae9be8838d4p-37; /* -1.13596475577881948265e-11 */

    /*
     * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
     * Input x is assumed to be bounded by ~pi/4 in magnitude.
     * Input y is the tail of x.
     *
     * Algorithm
     *        1. Since cos(-x) = cos(x), we need only to consider positive x.
     *        2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
     *        3. cos(x) is approximated by a polynomial of degree 14 on
     *           [0,pi/4]
     *                               4            14
     *                   cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
     *           where the remez error is
     *
     *         |              2     4     6     8     10    12     14 |     -58
     *         |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  )| <= 2
     *         |                                                           |
     *
     *                        4     6     8     10    12     14
     *        4. let r = C1*x +C2*x +C3*x +C4*x +C5*x  +C6*x  , then
     *               cos(x) = 1 - x*x/2 + r
     *           since cos(x+y) ~ cos(x) - sin(x)*y
     *                          ~ cos(x) - x*y,
     *           a correction term is necessary in cos(x) and hence
     *                cos(x+y) = 1 - (x*x/2 - (r - x*y))
     *           For better accuracy when x > 0.3, let qx = |x|/4 with
     *           the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
     *           Then
     *                cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
     *           Note that 1-qx and (x*x/2-qx) is EXACT here, and the
     *           magnitude of the latter is at least a quarter of x*x/2,
     *           thus, reducing the rounding error in the subtraction.
     */
    static private double __kernel_cos(double x, double y) {
        double a;
        double hz;
        double z;
        double r;
        double qx = zero;
        int ix;
        ix = __HI(x) & 0x7fffffff; /* ix = |x|'s high word*/
        if (ix < 0x3e400000) {

            /* if x < 2**27 */
            if (((int) x) == 0) {
                return one; /* generate inexact */
            }
        }
        z = x * x;
        r = z * (C1 +
            (z * (C2 + (z * (C3 + (z * (C4 + (z * (C5 + (z * C6))))))))));
        if (ix < 0x3FD33333) {

            /* if |x| < 0.3 */
            return one - ((0.5 * z) - ((z * r) - (x * y)));
        } else {
            if (ix > 0x3fe90000) { /* x > 0.78125 */
                qx = 0.28125;
            } else {
                qx = set(ix - 0x00200000, 0); /* x/4 */
            }
            hz = (0.5 * z) - qx;
            a = one - qx;
            return a - (hz - ((z * r) - (x * y)));
        }
    }

    /**
     *        Returns the tangent of its argument.
     *        @param        x        The argument, a double, assumed to be in radians.
     *        @return        Returns the tangent of x.
     */
    static public double tan(double x) {
        double z = zero;

        int n;
        int ix = __HI(x);

        /* |x| ~< pi/4 */
        ix &= 0x7fffffff;
        if (ix <= 0x3fe921fb) {
            return __kernel_tan(x, z, 1);
        } else if (ix >= 0x7ff00000) {

            /* tan(Inf or NaN) is NaN */
            return x - x; /* NaN */
        } else {

            /* argument reduction needed */
            double[] y = new double[2];
            n = __ieee754_rem_pio2(x, y);

            /*   1 -- n even -1 -- n odd */
            return __kernel_tan(y[0], y[1], 1 - ((n & 1) << 1));
        }
    }

    static private final double pio4 = 0x1.921fb54442d18p-1; /* 7.85398163397448278999e-01 */
    static private final double pio4lo = 0x1.1a62633145c07p-55; /* 3.06161699786838301793e-17 */
    static private final double[] T = {
            0x1.5555555555563p-2, /* 3.33333333333334091986e-01 */
            0x1.111111110fe7ap-3, /* 1.33333333333201242699e-01 */
            0x1.ba1ba1bb341fep-5, /* 5.39682539762260521377e-02 */
            0x1.664f48406d637p-6, /* 2.18694882948595424599e-02 */
            0x1.226e3e96e8493p-7, /* 8.86323982359930005737e-03 */
            0x1.d6d22c9560328p-9, /* 3.59207910759131235356e-03 */
            0x1.7dbc8fee08315p-10, /* 1.45620945432529025516e-03 */
            0x1.344d8f2f26501p-11, /* 5.88041240820264096874e-04 */
            0x1.026f71a8d1068p-12, /* 2.46463134818469906812e-04 */
            0x1.47e88a03792a6p-14, /* 7.81794442939557092300e-05 */
            0x1.2b80f32f0a7e9p-14, /* 7.14072491382608190305e-05 */
            -0x1.375cbdb605373p-16, /* -1.85586374855275456654e-05 */
            0x1.b2a7074bf7ad4p-16 /* 2.59073051863633712884e-05 */
        };

    /*
     *        __kernel_tan( x, y, k )
     * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
     * Input x is assumed to be bounded by ~pi/4 in magnitude.
     * Input y is the tail of x. * Input k indicates whether tan (if k=1) or
     * -1/tan (if k= -1) is returned. * * Algorithm
     *        1. Since tan(-x) = -tan(x), we need only to consider positive x.
     *        2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
     *        3. tan(x) is approximated by a odd polynomial of degree 27 on
     *           [0,0.67434]
     *                       3             27
     *                   tan(x) ~ x + T1*x + ... + T13*x
     *           where
     *
     *                 |tan(x)         2     4            26   |     -59.2
     *                 |----- - (1+T1*x +T2*x +.... +T13*x    )| <= 2
     *          |  x                                    |
     *
     *           Note: tan(x+y) = tan(x) + tan'(x)*y
     *                          ~ tan(x) + (1+x*x)*y
     *           Therefore, for better accuracy in computing tan(x+y), let
     *                     3      2      2       2       2
     *                r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
     *           then
     *                                     3    2
     *                tan(x+y) = x + (T1*x + (x *(r+y)+y)) *
     *   4. For x in [0.67434,pi/4],  let y = pi/4 - x, then
     *                tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
     *                       = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
     */
    static private double __kernel_tan(double x, double y, int iy) {
        double z;
        double r;
        double v;
        double w;
        double s;
        int ix;
        int hx;

        hx = __HI(x); /* high word of x */
        ix = hx & 0x7fffffff; /* high word of |x| */
        if (ix < 0x3e300000) { /* x < 2**-28 */
            if ((int) x == 0) { /* generate inexact */
                if (((ix | __LO(x)) | (iy + 1)) == 0) {
                    return one / abs(x);
                } else {
                    return (iy == 1) ? x : (-one / x);
                }
            }
        }
        if (ix >= 0x3FE59428) {

            /* |x|>=0.6744 */
            if (hx < 0) {
                x = -x;
                y = -y;
            }
            z = pio4 - x;
            w = pio4lo - y;
            x = z + w;
            y = 0.0;
        }
        z = x * x;
        w = z * z;

        /*
         *        Break x^5*(T[1]+x^2*T[2]+...) into
         *          x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
         *          x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
         */
        r = T[1] +
            (w * (T[3] +
            (w * (T[5] + (w * (T[7] + (w * (T[9] + (w * T[11])))))))));
        v = z * (T[2] +
            (w * (T[4] +
            (w * (T[6] + (w * (T[8] + (w * (T[10] + (w * T[12]))))))))));
        s = z * x;
        r = y + (z * ((s * (r + v)) + y));
        r += (T[0] * s);
        w = x + r;
        if (ix >= 0x3FE59428) {
            v = (double) iy;
            return (double) (1 - ((hx >> 30) & 2)) * (v -
            (2.0 * (x - (((w * w) / (w + v)) - r))));
        }
        if (iy == 1) {
            return w;
        } else {

            /* if allow error up to 2 ulp, simply return -1.0/(x+r) here */
            /*  compute -1.0/(x+r) accurately */
            double a;

            /* if allow error up to 2 ulp, simply return -1.0/(x+r) here */
            /*  compute -1.0/(x+r) accurately */
            double t;
            z = w;
            z = setLO(z, 0);
            v = r - (z - x);

            /* z+v = r+x */
            t = a = -1.0 / w;

            /* a = -1.0/w */
            t = setLO(t, 0);
            s = 1.0 + (t * z);
            return t + (a * (s + (t * v)));
        }
    }

    static private final double pio2_hi = 0x1.921fb54442d18p0; /* 1.57079632679489655800e+00 */
    static private final double pio2_lo = 0x1.1a62633145c07p-54; /* 6.12323399573676603587e-17 */
    static private final double pio4_hi = 0x1.921fb54442d18p-1; /* 7.85398163397448278999e-01 */

    /* coefficient for R(x^2) */
    static private final double pS0 = 0x1.5555555555555p-3; /*  1.66666666666666657415e-01 */
    static private final double pS1 = -0x1.4d61203eb6f7dp-2; /* -3.25565818622400915405e-01 */
    static private final double pS2 = 0x1.9c1550e884455p-3; /*  2.01212532134862925881e-01 */
    static private final double pS3 = -0x1.48228b5688f3bp-5; /* -4.00555345006794114027e-02 */
    static private final double pS4 = 0x1.9efe07501b288p-11; /*  7.91534994289814532176e-04 */
    static private final double pS5 = 0x1.23de10dfdf709p-15; /*  3.47933107596021167570e-05 */
    static private final double qS1 = -0x1.33a271c8a2d4bp1; /* -2.40339491173441421878e+00 */
    static private final double qS2 = 0x1.02ae59c598ac8p1; /*  2.02094576023350569471e+00 */
    static private final double qS3 = -0x1.6066c1b8d0159p-1; /* -6.88283971605453293030e-01 */
    static private final double qS4 = 0x1.3b8c5b12e9282p-4; /*  7.70381505559019352791e-02 */

    /*
     *        asin(x)
     * Method :
     *        Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
     *        we approximate asin(x) on [0,0.5] by
     *                asin(x) = x + x*x^2*R(x^2)
     *        where
     *                R(x^2) is a rational approximation of (asin(x)-x)/x^3
     *        and its remez error is bounded by
     *                |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
     *
     *        For x in [0.5,1]
     *                asin(x) = pi/2-2*asin(sqrt((1-x)/2))
     *        Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
     *        then for x>0.98
     *                asin(x) = pi/2 - 2*(s+s*z*R(z))
     *                        = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
     *        For x<=0.98, let pio4_hi = pio2_hi/2, then
     *                f = hi part of s;
     *                c = sqrt(z) - f = (z-f*f)/(s+f)         ...f+c=sqrt(z)
     *        and
     *                asin(x) = pi/2 - 2*(s+s*z*R(z))
     *                        = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
     *                        = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
     *
     * Special cases:
     *        if x is NaN, return x itself;
     *        if |x|>1, return NaN with invalid signal.
     *
     */

    /**
     *        Returns the inverse (arc) sine of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns the angle, in radians, whose sine is x.
     *                        It is in the range [-pi/2,pi/2].
     */
    static public double asin(double x) {
        double t = zero;
        double w;
        double p;
        double q;
        double c;
        double r;
        double s;
        int hx;
        int ix;
        hx = __HI(x);

        ix = hx & 0x7fffffff;
        if (ix >= 0x3ff00000) { /* |x|>= 1 */
            if (((ix - 0x3ff00000) | __LO(x)) == 0) {

                /* asin(1)=+-pi/2 with inexact */
                return (x * pio2_hi) + (x * pio2_lo);
            }
            return (x - x) / (x - x); /* asin(|x|>1) is NaN */
        } else if (ix < 0x3fe00000) { /* |x|<0.5 */
            if (ix < 0x3e400000) { /* if |x| < 2**-27 */
                if ((huge + x) > one) {
                    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 = one + (t * (qS1 + (t * (qS2 + (t * (qS3 + (t * qS4)))))));
            w = p / q;
            return x + (x * w);
        }

        /* 1> |x|>= 0.5 */
        w = one - abs(x);
        t = w * 0.5;
        p = t * (pS0 +
            (t * (pS1 + (t * (pS2 + (t * (pS3 + (t * (pS4 + (t * pS5))))))))));
        q = one + (t * (qS1 + (t * (qS2 + (t * (qS3 + (t * qS4)))))));
        s = sqrt(t);
        if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
            w = p / q;
            t = pio2_hi - ((2.0 * (s + (s * w))) - pio2_lo);
        } else {
            w = s;
            w = setLO(w, 0);
            c = (t - (w * w)) / (s + w);
            r = p / q;
            p = (2.0 * s * r) - (pio2_lo - (2.0 * c));
            q = pio4_hi - (2.0 * w);
            t = pio4_hi - (p - q);
        }
        return ((hx > 0) ? t : (-t));
    }

    /*
     * Method :
     *        acos(x)  = pi/2 - asin(x)
     *        acos(-x) = pi/2 + asin(x)
     * For |x|<=0.5
     *        acos(x) = pi/2 - (x + x*x^2*R(x^2))        (see asin.c)
     * For x>0.5
     *         acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
     *                = 2asin(sqrt((1-x)/2))
     *                = 2s + 2s*z*R(z)         ...z=(1-x)/2, s=sqrt(z)
     *                = 2f + (2c + 2s*z*R(z))
     *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
     *     for f so that f+c ~ sqrt(z).
     * For x<-0.5
     *        acos(x) = pi - 2asin(sqrt((1-|x|)/2))
     *          = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
     *
     * Special cases:
     *        if x is NaN, return x itself;
     *        if |x|>1, return NaN with invalid signal.
     *
     * Function needed: sqrt
     */

    /**
     *        Returns the inverse (arc) cosine of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns the angle, in radians, whose cosine is x.
     *                        It is in the range [0,pi].
     */
    static public double acos(double x) {
        double z;
        double p;
        double q;
        double r;
        double w;
        double s;
        double c;
        double df;
        int hx;
        int ix;
        hx = __HI(x);
        ix = hx & 0x7fffffff;
        if (ix >= 0x3ff00000) { /* |x| >= 1 */
            if (((ix - 0x3ff00000) | __LO(x)) == 0) { /* |x|==1 */
                if (hx > 0) {
                    return 0.0; /* acos(1) = 0  */
                } else {
                    return PI + (2.0 * pio2_lo); /* acos(-1)= pi */
                }
            }
            return (x - x) / (x - x); /* acos(|x|>1) is NaN */
        }
        if (ix < 0x3fe00000) { /* |x| < 0.5 */
            if (ix <= 0x3c600000) {
                return pio2_hi + pio2_lo; /*if|x|<2**-57*/
            }
            z = x * x;
            p = z * (pS0 +
                (z * (pS1 +
                (z * (pS2 + (z * (pS3 + (z * (pS4 + (z * pS5))))))))));
            q = one + (z * (qS1 + (z * (qS2 + (z * (qS3 + (z * qS4)))))));
            r = p / q;
            return pio2_hi - (x - (pio2_lo - (x * r)));
        } else if (hx < 0) { /* x < -0.5 */
            z = (one + x) * 0.5;
            p = z * (pS0 +
                (z * (pS1 +
                (z * (pS2 + (z * (pS3 + (z * (pS4 + (z * pS5))))))))));
            q = one + (z * (qS1 + (z * (qS2 + (z * (qS3 + (z * qS4)))))));
            s = sqrt(z);
            r = p / q;
            w = (r * s) - pio2_lo;
            return PI - (2.0 * (s + w));
        } else { /* x > 0.5 */
            z = (one - x) * 0.5;
            s = sqrt(z);
            df = s;
            df = setLO(df, 0);
            c = (z - (df * df)) / (s + df);
            p = z * (pS0 +
                (z * (pS1 +
                (z * (pS2 + (z * (pS3 + (z * (pS4 + (z * pS5))))))))));
            q = one + (z * (qS1 + (z * (qS2 + (z * (qS3 + (z * qS4)))))));
            r = p / q;
            w = (r * s) + c;
            return 2.0 * (df + w);
        }
    }

    static private final double[] atanhi = {
            0x1.dac670561bb4fp-2, /* 4.63647609000806093515e-01 atan(0.5)hi  */
            0x1.921fb54442d18p-1, /* 7.85398163397448278999e-01 atan(1.0)hi  */
            0x1.f730bd281f69bp-1, /* 9.82793723247329054082e-01 atan(1.5)hi  */
            0x1.921fb54442d18p0   /* 1.57079632679489655800e+00 atan(inf)hi  */
        };
    static private final double[] atanlo = {
            0x1.a2b7f222f65e2p-56, /* 2.26987774529616870924e-17 atan(0.5)lo  */
            0x1.1a62633145c07p-55, /* 3.06161699786838301793e-17 atan(1.0)lo  */
            0x1.007887af0cbbdp-56, /* 1.39033110312309984516e-17 atan(1.5)lo  */
            0x1.1a62633145c07p-54  /* 6.12323399573676603587e-17 atan(inf)lo  */
        };
    static private final double[] aT = {
            0x1.555555555550dp-2,  /*  3.33333333333329318027e-01 */
            -0x1.999999998ebc4p-3, /* -1.99999999998764832476e-01 */
            0x1.24924920083ffp-3,  /*  1.42857142725034663711e-01 */
            -0x1.c71c6fe231671p-4, /* -1.11111104054623557880e-01 */
            0x1.745cdc54c206ep-4,  /*  9.09088713343650656196e-02 */
            -0x1.3b0f2af749a6dp-4, /* -7.69187620504482999495e-02 */
            0x1.10d66a0d03d51p-4,  /*  6.66107313738753120669e-02 */
            -0x1.dde2d52defd9ap-5, /* -5.83357013379057348645e-02 */
            0x1.97b4b24760debp-5,  /*  4.97687799461593236017e-02 */
            -0x1.2b4442c6a6c2fp-5, /* -3.65315727442169155270e-02 */
            0x1.0ad3ae322da11p-6   /*  1.62858201153657823623e-02 */
        };

    /*
     * Method
     *   1. Reduce x to positive by atan(x) = -atan(-x).
     *   2. According to the integer k=4t+0.25 chopped, t=x, the argument
     *      is further reduced to one of the following intervals and the
     *      arctangent of t is evaluated by the corresponding formula:
     *
     *      [0,7/16]      atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
     *      [7/16,11/16]  atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
     *      [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
     *      [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
     *      [39/16,INF]   atan(x) = atan(INF) + atan( -1/t )
     *
     * Constants:
     * The hexadecimal values are the intended ones for the following
     * constants. The decimal values may be used, provided that the
     * compiler will convert from decimal to binary accurately enough
     * to produce the hexadecimal values shown.
     */

    /**
     *        Returns the inverse (arc) tangent of its argument.
     *        @param        x        The argument, a double.
     *        @return        Returns the angle, in radians, whose tangent is x.
     *                        It is in the range [-pi/2,pi/2].
     */
    static public double atan(double x) {
        double w;
        double s1;
        double s2;
        double z;
        int ix;
        int hx;
        int id;

        hx = __HI(x);
        ix = hx & 0x7fffffff;
        if (ix >= 0x44100000) { /* if |x| >= 2^66 */
            if ((ix > 0x7ff00000) || ((ix == 0x7ff00000) && (__LO(x) != 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) > one) {
                    return x; /* raise inexact */
                }
            }
            id = -1;
        } else {
            x = abs(x);
            if (ix < 0x3ff30000) { /* |x| < 1.1875 */
                if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
                    id = 0;
                    x = ((2.0 * x) - one) / (2.0 + x);
                } else { /* 11/16<=|x|< 19/16 */
                    id = 1;
                    x = (x - one) / (x + one);
                }
            } else {
                if (ix < 0x40038000) { /* |x| < 2.4375 */
                    id = 2;
                    x = (x - 1.5) / (one + (1.5 * x));
                } else { /* 2.4375 <= |x| < 2^66 */
                    id = 3;
                    x = -1.0 / x;
                }
            }
        }

        /* end of argument reduction */
        z = x * x;
        w = z * z;

        /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
        s1 = z * (aT[0] +
            (w * (aT[2] +
            (w * (aT[4] + (w * (aT[6] + (w * (aT[8] + (w * aT[10]))))))))));
        s2 = w * (aT[1] +
            (w * (aT[3] + (w * (aT[5] + (w * (aT[7] + (w * aT[9]))))))));
        if (id < 0) {
            return x - (x * (s1 + s2));
        } else {
            z = atanhi[id] - (((x * (s1 + s2)) - atanlo[id]) - x);
            return (hx < 0) ? (-z) : z;
        }
    }

    static private final double pi_o_4 = 0x1.921fb54442d18p-1; /* 7.8539816339744827900e-01 */
    static private final double pi_o_2 = 0x1.921fb54442d18p0; /* 1.5707963267948965580e+00 */
    static private final double pi_lo = 0x1.1a62633145c07p-53; /* 1.2246467991473531772e-16 */

    /*
     * Method :
     *        1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
     *        2. Reduce x to positive by (if x and y are unexceptional):
     *                ARG (x+iy) = arctan(y/x)              ... if x > 0,
     *                ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
     *
     * Special cases:
     *
     *        ATAN2((anything), NaN ) is NaN;
     *        ATAN2(NAN , (anything) ) is NaN;
     *        ATAN2(+-0, +(anything but NaN)) is +-0  ;
     *        ATAN2(+-0, -(anything but NaN)) is +-pi ;
     *        ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
     *        ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
     *        ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
     *        ATAN2(+-INF,+INF ) is +-pi/4 ;
     *        ATAN2(+-INF,-INF ) is +-3pi/4;
     *        ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
     *
     * Constants:
     * The hexadecimal values are the intended ones for the following
     * constants. The decimal values may be used, provided that the
     * compiler will convert from decimal to binary accurately enough
     * to produce the hexadecimal values shown.
     */

    /**
     *        Returns angle corresponding to a Cartesian point.
     *        @param        x        The first argument, a double.
     *        @param        y        The second argument, a double.
     *        @return        Returns the angle, in radians, the the line
     *                        from (0,0) to (x,y) makes with the x-axis.
     *                        It is in the range [-pi,pi].
     */
    static public double atan2(double y, double x) {
        double z;
        int k;
        int m;
        int hx;
        int hy;
        int ix;
        int iy;
        int lx;
        int ly;

        hx = __HI(x);
        ix = hx & 0x7fffffff;
        lx = __LO(x);
        hy = __HI(y);
        iy = hy & 0x7fffffff;
        ly = __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 atan(y); /* x=1.0 */
        }
        m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */

        /* when y = 0 */
        if ((iy | ly) == 0) {
            switch (m) {
            case 0:
            case 1:
                return y; /* atan(+-0,+anything)=+-0 */
            case 2:
                return PI + tiny; /* atan(+0,-anything) = pi */
            case 3:
                return -PI - tiny; /* atan(-0,-anything) =-pi */
            }
        }

        /* when x = 0 */
        if ((ix | lx) == 0) {
            return ((hy < 0) ? (-pi_o_2 - tiny) : (pi_o_2 + tiny));
        }

        /* when x is INF */
        if (ix == 0x7ff00000) {
            if (iy == 0x7ff00000) {
                switch (m) {
                case 0:
                    return pi_o_4 + tiny; /* atan(+INF,+INF) */
                case 1:
                    return -pi_o_4 - tiny; /* atan(-INF,+INF) */
                case 2:
                    return (3.0 * pi_o_4) + tiny; /*atan(+INF,-INF)*/
                case 3:
                    return (-3.0 * pi_o_4) - tiny; /*atan(-INF,-INF)*/
                }
            } else {
                switch (m) {
                case 0:
                    return zero; /* atan(+...,+INF) */
                case 1:
                    return -zero; /* atan(-...,+INF) */
                case 2:
                    return PI + tiny; /* atan(+...,-INF) */
                case 3:
                    return -PI - tiny; /* atan(-...,-INF) */
                }
            }
        }

        /* when y is INF */
        if (iy == 0x7ff00000) {
            return (hy < 0) ? (-pi_o_2 - tiny) : (pi_o_2 + tiny);
        }

        /* compute y/x */
        k = (iy - ix) >> 20;
        if (k > 60) {

            /* |y/x| >  2**60 */
            z = pi_o_2 + (0.5 * pi_lo);
        } else if ((hx < 0) && (k < -60)) {

            /* |y|/x < -2**60 */
            z = 0.0;
        } else {

            /* safe to do y/x */
            z = atan(abs(y / x));
        }

        switch (m) {
        case 0:
            return z; /* atan(+,+) */
        case 1:
            return setHI(z, __HI(z) ^ 0x80000000); /* atan(-,+) */
        case 2:
            return PI - (z - pi_lo); /* atan(+,-) */
        default: /* case 3 */
            return (z - pi_lo) - PI; /* atan(-,-) */
        }
    }

    /*
     * kernel function:
     *        __kernel_sin                ... sine function on [-pi/4,pi/4]
     *        __ieee754_rem_pio2        ... argument reduction routine
     *
     * Method.
     *      Let S,C and T denote the sin, cos and tan respectively on
     *        [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
     *        in [-pi/4 , +pi/4], and let n = k mod 4.
     *        We have
     *
     *          n        sin(x)      cos(x)        tan(x)
     *     ----------------------------------------------------------
     *            0               S           C                 T
     *            1               C          -S                -1/T
     *            2              -S          -C                 T
     *            3              -C           S                -1/T
     *     ----------------------------------------------------------
     *
     * Special cases:
     *      Let trig be any of sin, cos, or tan.
     *      trig(+-INF)  is NaN, with signals;
     *      trig(NaN)    is that NaN;
     *
     * Accuracy:
     *        TRIG(x) returns trig(x) nearly rounded
     */
    /*
     * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
     */
    static private final int[] two_over_pi = {
            0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, 0x95993c,
            0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, 0x424dd2, 0xe00649,
            0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, 0xa73ee8, 0x8235f5, 0x2ebb44,
            0x84e99c, 0x7026b4, 0x5f7e41, 0x3991d6, 0x398353, 0x39f49c, 0x845f8b,
            0xbdf928, 0x3b1ff8, 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d,
            0x367ecf, 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
            0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, 0x560330,
            0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, 0x91615e, 0xe61b08,
            0x659985, 0x5f14a0, 0x68408d, 0xffd880, 0x4d7327, 0x310606, 0x1556ca,
            0x73a8c9, 0x60e27b, 0xc08c6b
        };
    static private final int[] npio2_hw = {
            0x3ff921fb, 0x400921fb, 0x4012d97c, 0x401921fb, 0x401f6a7a,
            0x4022d97c, 0x4025fdbb, 0x402921fb, 0x402c463a, 0x402f6a7a,
            0x4031475c, 0x4032d97c, 0x40346b9c, 0x4035fdbb, 0x40378fdb,
            0x403921fb, 0x403ab41b, 0x403c463a, 0x403dd85a, 0x403f6a7a,
            0x40407e4c, 0x4041475c, 0x4042106c, 0x4042d97c, 0x4043a28c,
            0x40446b9c, 0x404534ac, 0x4045fdbb, 0x4046c6cb, 0x40478fdb,
            0x404858eb, 0x404921fb
        };
    static private final double zero = 0.00000000000000000000e+00; // 0x0000000000000000 
    static private final double half = 0x1.0p-1; /* 5.00000000000000000000e-01 */
    static private final double two24 = 0x1.0p24; /* 1.67772160000000000000e+07 */
    static private final double invpio2 = 0x1.45f306dc9c883p-1; /* 6.36619772367581382433e-01 53 bits of 2/pi */
    static private final double pio2_1 = 0x1.921fb544p0; /* 1.57079632673412561417e+00 first  33 bit of pi/2 */
    static private final double pio2_1t = 0x1.0b4611a626331p-34; /* 6.07710050650619224932e-11 pi/2 - pio2_1 */
    static private final double pio2_2 = 0x1.0b4611a6p-34; /* 6.07710050630396597660e-11 second 33 bit of pi/2 */
    static private final double pio2_2t = 0x1.3198a2e037073p-69; /* 2.02226624879595063154e-21 pi/2 - (pio2_1+pio2_2) */
    static private final double pio2_3 = 0x1.3198a2ep-69; /* 2.02226624871116645580e-21 third  33 bit of pi/2 */
    static private final double pio2_3t = 0x1.b839a252049c1p-104; /* 8.47842766036889956997e-32 pi/2 - (pio2_1+pio2_2+pio2_3) */

    /*
     *        Return the remainder of x % pi/2 in y[0]+y[1]
     */
    static private int __ieee754_rem_pio2(double x, double[] y) {
        double z = zero;
        double w;
        double t;
        double r;
        double fn;
        int i;
        int j;
        int nx;
        int n;
        int ix;
        int hx;

        hx = __HI(x); /* high word of x */
        ix = hx & 0x7fffffff;
        if (ix <= 0x3fe921fb) {

            /* |x| ~<= pi/4 , no need for reduction */
            y[0] = x;
            y[1] = 0;
            return 0;
        }

        if (ix < 0x4002d97c) {

            /* |x| < 3pi/4, special case with n=+-1 */
            if (hx > 0) {
                z = x - pio2_1;
                if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
                    y[0] = z - pio2_1t;
                    y[1] = (z - y[0]) - pio2_1t;
                } else { /* near pi/2, use 33+33+53 bit pi */
                    z -= pio2_2;
                    y[0] = z - pio2_2t;
                    y[1] = (z - y[0]) - pio2_2t;
                }
                return 1;
            } else { /* negative x */
                z = x + pio2_1;
                if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
                    y[0] = z + pio2_1t;
                    y[1] = (z - y[0]) + pio2_1t;
                } else { /* near pi/2, use 33+33+53 bit pi */
                    z += pio2_2;
                    y[0] = z + pio2_2t;
                    y[1] = (z - y[0]) + pio2_2t;
                }
                return -1;
            }
        }

        if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
            t = abs(x);
            n = (int) ((t * invpio2) + half);
            fn = (double) n;
            r = t - (fn * pio2_1);
            w = fn * pio2_1t; /* 1st round good to 85 bit */
            if ((n < 32) && (ix != npio2_hw[n - 1])) {
                y[0] = r - w; /* quick check no cancellation */
            } else {
                j = ix >> 20;
                y[0] = r - w;
                i = j - (((__HI(y[0])) >> 20) & 0x7ff);
                if (i > 16) { /* 2nd iteration needed, good to 118 */
                    t = r;
                    w = fn * pio2_2;
                    r = t - w;
                    w = (fn * pio2_2t) - ((t - r) - w);
                    y[0] = r - w;
                    i = j - (((__HI(y[0])) >> 20) & 0x7ff);
                    if (i > 49) { /* 3rd iteration need, 151 bits acc */
                        t = r; /* will cover all possible cases */
                        w = fn * pio2_3;
                        r = t - w;
                        w = (fn * pio2_3t) - ((t - r) - w);
                        y[0] = r - w;
                    }
                }
            }
            y[1] = (r - y[0]) - w;
            if (hx < 0) {
                y[0] = -y[0];
                y[1] = -y[1];
                return -n;
            } else {
                return n;
            }
        }

        /*
         * all other (large) arguments
         */
        if (ix >= 0x7ff00000) {

            /* x is inf or NaN */
            y[0] = y[1] = x - x;
            return 0;
        }

        /* set z = scalbn(|x|,ilogb(x)-23) */
        double[] tx = new double[3];
        long lx = Double.doubleToLongBits(x);
        long exp = (0x7ff0000000000000L & lx) >> 52;
        exp -= 1046;
        lx -= (exp << 52);
        lx &= 0x7fffffffffffffffL;
        z = Double.longBitsToDouble(lx);
        for (i = 0; i < 2; i++) {
            tx[i] = (double) ((int) (z));
            z = (z - tx[i]) * two24;
        }
        tx[2] = z;
        nx = 3;
        while (tx[nx - 1] == zero)
            nx--; /* skip zero term */
        n = __kernel_rem_pio2(tx, y, (int) exp, nx);
        //System.out.println("KERNEL");
        //System.out.println("tx "+tx[0]+"  "+tx[1]+"  "+tx[2]);
        //System.out.println("y "+y[0]+"  "+y[1]);
        //System.out.println("exp "+(int)exp+"  nx "+nx+"   n "+n);
        if (hx < 0) {
            y[0] = -y[0];
            y[1] = -y[1];
            return -n;
        }
        return n;
    }

    /*
     * __kernel_rem_pio2(x,y,e0,nx)
     * double x[],y[]; int e0,nx; int two_over_pi[];
     *
     * __kernel_rem_pio2 return the last three digits of N with
     *                y = x - N*pi/2
     * so that |y| < pi/2.
     *
     * The method is to compute the integer (mod 8) and fraction parts of
     * (2/pi)*x without doing the full multiplication. In general we
     * skip the part of the product that are known to be a huge integer
     * (more accurately, = 0 mod 8 ). Thus the number of operations are
     * independent of the exponent of the input.
     *
     * (2/pi) is represented by an array of 24-bit integers in two_over_pi[].
     *
     * Input parameters:
     *         x[]        The input value (must be positive) is broken into nx
     *                pieces of 24-bit integers in double precision format.
     *                x[i] will be the i-th 24 bit of x. The scaled exponent
     *                of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
     *                match x's up to 24 bits.
     *
     *                Example of breaking a double positive z into x[0]+x[1]+x[2]:
     *                        e0 = ilogb(z)-23
     *                        z  = scalbn(z,-e0)
     *                for i = 0,1,2
     *                        x[i] = floor(z)
     *                        z    = (z-x[i])*2**24
     *
     *
     *        y[]        ouput result in an array of double precision numbers.
     *                The dimension of y[] is:
     *                        24-bit  precision        1
     *                        53-bit  precision        2
     *                        64-bit  precision        2
     *                        113-bit precision        3
     *                The actual value is the sum of them. Thus for 113-bit
     *                precison, one may have to do something like:
     *
     *                long double t,w,r_head, r_tail;
     *                t = (long double)y[2] + (long double)y[1];
     *                w = (long double)y[0];
     *                r_head = t+w;
     *                r_tail = w - (r_head - t);
     *
     *        e0        The exponent of x[0]
     *
     *        nx        dimension of x[]
     *
     *          prec        an integer indicating the precision:
     *                        0        24  bits (single)
     *                        1        53  bits (double)
     *                        2        64  bits (extended)
     *                        3        113 bits (quad)
     *
     *        two_over_pi[]
     *                integer array, contains the (24*i)-th to (24*i+23)-th
     *                bit of 2/pi after binary point. The corresponding
     *                floating value is
     *
     *                        two_over_pi[i] * 2^(-24(i+1)).
     *
     * External function:
     *        double scalbn(), floor();
     *
     *
     * Here is the description of some local variables:
     *
     *         jk        jk+1 is the initial number of terms of two_over_pi[] needed
     *                in the computation. The recommended value is 2,3,4,
     *                6 for single, double, extended,and quad.
     *
     *         jz        local integer variable indicating the number of
     *                terms of two_over_pi[] used.
     *
     *        jx        nx - 1
     *
     *        jv        index for pointing to the suitable two_over_pi[] for the
     *                computation. In general, we want
     *                        ( 2^e0*x[0] * two_over_pi[jv-1]*2^(-24jv) )/8
     *                is an integer. Thus
     *                        e0-3-24*jv >= 0 or (e0-3)/24 >= jv
     *                Hence jv = max(0,(e0-3)/24).
     *
     *        jp        jp+1 is the number of terms in PIo2[] needed, jp = jk.
     *
     *         q[]        double array with integral value, representing the
     *                24-bits chunk of the product of x and 2/pi.
     *
     *        q0        the corresponding exponent of q[0]. Note that the
     *                exponent for q[i] would be q0-24*i.
     *
     *        PIo2[]        double precision array, obtained by cutting pi/2
     *                into 24 bits chunks.
     *
     *        f[]        two_over_pi[] in floating point
     *
     *        iq[]        integer array by breaking up q[] in 24-bits chunk.
     *
     *        fq[]        final product of x*(2/pi) in fq[0],..,fq[jk]
     *
     *        ih        integer. If >0 it indicates q[] is >= 0.5, hence
     *                it also indicates the *sign* of the result.
     *
     */
    /*
     * Constants:
     * The hexadecimal values are the intended ones for the following
     * constants. The decimal values may be used, provided that the
     * compiler will convert from decimal to binary accurately enough
     * to produce the hexadecimal values shown.
     */
    static final private double[] PIo2 = {
            0x1.921fb4p0, /* 1.57079625129699707031e+00 */
            0x1.4442dp-24, /* 7.54978941586159635335e-08 */
            0x1.846988p-48, /* 5.39030252995776476554e-15 */
            0x1.8cc516p-72, /* 3.28200341580791294123e-22 */
            0x1.01b838p-96, /* 1.27065575308067607349e-29 */
            0x1.a25204p-120, /* 1.22933308981111328932e-36 */
            0x1.382228p-145, /* 2.73370053816464559624e-44 */
            0x1.9f31dp-169 /* 2.16741683877804819444e-51 */
        };
    static final private double twon24 = 0x1.0p-24; /* 5.96046447753906250000e-08 */

    static private int __kernel_rem_pio2(double[] x, double[] y, int e0, int nx) {
        int jz;
        int jx;
        int jv;
        int jp;
        int jk;
        int carry;
        int n;
        int i;
        int j;
        int k;
        int m;
        int q0;
        int ih;
        double z;
        double fw;
        double[] f = new double[20];
        double[] q = new double[20];
        double[] fq = new double[20];
        int[] iq = new int[20];

        /* initialize jk*/
        jk = 4;
        jp = jk;

        /* determine jx,jv,q0, note that 3>q0 */
        jx = nx - 1;
        jv = (e0 - 3) / 24;
        if (jv < 0) {
            jv = 0;
        }
        q0 = e0 - (24 * (jv + 1));

        /* set up f[0] to f[jx+jk] where f[jx+jk] = two_over_pi[jv+jk] */
        j = jv - jx;
        m = jx + jk;
        for (i = 0; i <= m; i++, j++)
            f[i] = ((j < 0) ? zero : (double) two_over_pi[j]);

        /* compute q[0],q[1],...q[jk] */
        for (i = 0; i <= jk; i++) {
            for (j = 0, fw = 0.0; j <= jx; j++)
                fw += (x[j] * f[(jx + i) - j]);
            q[i] = fw;
        }

        jz = jk;

        while (true) { // recompute:

            /* distill q[] into iq[] reversingly */
            for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
                fw = (double) ((int) (twon24 * z));
                iq[i] = (int) (z - (two24 * fw));
                z = q[j - 1] + fw;
            }

            /* compute n */
            z = scalbn(z, q0); /* actual value of z */
            z -= (8.0 * floor(z * 0.125)); /* trim off integer >= 8 */
            n = (int) z;
            z -= (double) n;
            ih = 0;
            if (q0 > 0) { /* need iq[jz-1] to determine n */
                i = (iq[jz - 1] >> (24 - q0));
                n += i;
                iq[jz - 1] -= (i << (24 - q0));
                ih = iq[jz - 1] >> (23 - q0);
            } else if (q0 == 0) {
                ih = iq[jz - 1] >> 23;
            } else if (z >= 0.5) {
                ih = 2;
            }
            if (ih > 0) { /* q > 0.5 */
                n += 1;
                carry = 0;
                for (i = 0; i < jz; i++) { /* compute 1-q */
                    j = iq[i];
                    if (carry == 0) {
                        if (j != 0) {
                            carry = 1;
                            iq[i] = 0x1000000 - j;
                        }
                    } else {
                        iq[i] = 0xffffff - j;
                    }
                }
                if (q0 > 0) { /* rare case: chance is 1 in 12 */
                    switch (q0) {
                    case 1:
                        iq[jz - 1] &= 0x7fffff;
                        break;
                    case 2:
                        iq[jz - 1] &= 0x3fffff;
                        break;
                    }
                }
                if (ih == 2) {
                    z = one - z;
                    if (carry != 0) {
                        z -= scalbn(one, q0);
                    }
                }
            }

            /* check if recomputation is needed */
            if (z == zero) {
                j = 0;
                for (i = jz - 1; i >= jk; i--)
                    j |= iq[i];
                if (j == 0) { /* need recomputation */
                    for (k = 1; iq[jk - k] == 0; k++)
                        ; /* k = no. of terms needed */for (i = jz + 1;
                            i <= (jz + k); i++) { /* add q[jz+1] to q[jz+k] */
                        f[jx + i] = (double) two_over_pi[jv + i];
                        for (j = 0, fw = 0.0; j <= jx; j++)
                            fw += (x[j] * f[(jx + i) - j]);
                        q[i] = fw;
                    }
                    jz += k;
                    continue; //goto recompute;
                }
            }
            break;
        }

        /* chop off zero terms */
        if (z == 0.0) {
            jz--;
            q0 -= 24;
            while (iq[jz] == 0) {
                jz--;
                q0 -= 24;
            }
        } else { /* break z into 24-bit if necessary */
            z = scalbn(z, -q0);
            if (z >= two24) {
                fw = (double) ((int) (twon24 * z));
                iq[jz] = (int) (z - (two24 * fw));
                jz++;
                q0 += 24;
                iq[jz] = (int) fw;
            } else {
                iq[jz] = (int) z;
            }
        }

        /* convert integer "bit" chunk to floating-point value */
        fw = scalbn(one, q0);
        for (i = jz; i >= 0; i--) {
            q[i] = fw * (double) iq[i];
            fw *= twon24;
        }

        /* compute PIo2[0,...,jp]*q[jz,...,0] */
        for (i = jz; i >= 0; i--) {
            for (fw = 0.0, k = 0; (k <= jp) && (k <= (jz - i)); k++)
                fw += (PIo2[k] * q[i + k]);
            fq[jz - i] = fw;
        }

        /* compress fq[] into y[] */
        fw = 0.0;
        for (i = jz; i >= 0; i--)
            fw += fq[i];
        y[0] = (ih == 0) ? fw : (-fw);
        fw = fq[0] - fw;
        for (i = 1; i <= jz; i++)
            fw += fq[i];
        y[1] = ((ih == 0) ? fw : (-fw));
        return n & 7;
    }

    static final private double[] bp = { 1.0, 1.5, };
    static final private double[] dp_h = {
            0.0, 0x1.2b8034p-1
        }; /* 5.84962487220764160156e-01 */
    static final private double[] dp_l = {
            0.0, 0x1.cfdeb43cfd006p-27
        }; /* 1.35003920212974897128e-08 */
    static final private double two53 = 0x1.0p53; /* 9007199254740992.0 */

    /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
    static final private double L1 = 0x1.3333333333303p-1; /* 5.99999999999994648725e-01 */
    static final private double L2 = 0x1.b6db6db6fabffp-2; /*  4.28571428578550184252e-01 */
    static final private double L3 = 0x1.55555518f264dp-2; /*  3.33333329818377432918e-01 */
    static final private double L4 = 0x1.17460a91d4101p-2; /*  2.72728123808534006489e-01 */
    static final private double L5 = 0x1.d864a93c9db65p-3; /*  2.30660745775561754067e-01 */
    static final private double L6 = 0x1.a7e284a454eefp-3; /*  2.06975017800338417784e-01 */
    static final private double lg2 = 0x1.62e42fefa39efp-1; /*  6.93147180559945286227e-01 */
    static final private double lg2_h = 0x1.62e43p-1; /*  6.93147182464599609375e-01 */
    static final private double lg2_l = -1.90465429995776804525e-09; /* 0xbe205c610ca86c39 */
    static final private double ovt = 8.0085662595372944372e-17; /* -(1024-log2(ovfl+.5ulp)) */
    static final private double cp = 0x1.ec709dc3a03fdp-1; /*  9.61796693925975554329e-01 = 2/(3ln2) */
    static final private double cp_h = 0x1.ec709ep-1; /*  9.61796700954437255859e-01 = (float)cp */
    static final private double cp_l = -0x1.e2fe0145b01f5p-28; /* -7.02846165095275826516e-09 = tail of cp_h*/
    static final private double ivln2 = 0x1.71547652b82fep0; /*  1.44269504088896338700e+00 = 1/ln2 */
    static final private double ivln2_h = 0x1.715476p0; /*  1.44269502162933349609e+00 = 24b 1/ln2*/
    static final private double ivln2_l = 0x1.4ae0bf85ddf44p-26; /*  1.92596299112661746887e-08 = 1/ln2 tail*/

    /*
     *        Returns x to the power y
     *
     *                      n
     * Method:  Let x =  2   * (1+f)
     *        1. Compute and return log2(x) in two pieces:
     *                log2(x) = w1 + w2,
     *           where w1 has 53-24 = 29 bit trailing zeros.
     *        2. Perform y*log2(x) = n+y' by simulating muti-precision
     *           arithmetic, where |y'|<=0.5.
     *        3. Return x**y = 2**n*exp(y'*log2)
     *
     * Special cases:
     *        1.  (anything) ** 0  is 1
     *        2.  (anything) ** 1  is itself
     *        3.  (anything) ** NAN is NAN
     *        4.  NAN ** (anything except 0) is NAN
     *        5.  +-(|x| > 1) **  +INF is +INF
     *        6.  +-(|x| > 1) **  -INF is +0
     *        7.  +-(|x| < 1) **  +INF is +0
     *        8.  +-(|x| < 1) **  -INF is +INF
     *        9.  +-1         ** +-INF is NAN
     *        10. +0 ** (+anything except 0, NAN)               is +0
     *        11. -0 ** (+anything except 0, NAN, odd integer)  is +0
     *        12. +0 ** (-anything except 0, NAN)               is +INF
     *        13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
     *        14. -0 ** (odd integer) = -( +0 ** (odd integer) )
     *        15. +INF ** (+anything except 0,NAN) is +INF
     *        16. +INF ** (-anything except 0,NAN) is +0
     *        17. -INF ** (anything)  = -0 ** (-anything)
     *        18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
     *        19. (-anything except 0 and inf) ** (non-integer) is NAN
     *
     * Accuracy:
     *        pow(x,y) returns x**y nearly rounded. In particular
     *                        pow(integer,integer)
     *        always returns the correct integer provided it is
     *        representable.
     */

    /**
     *        Returns x to the power y.
     *        @param        x        The base.
     *        @param        y        The exponent.
     *        @return        x to the power y.
     */
    static public double pow(double x, double y) {
        double z;
        double ax;
        double z_h;
        double z_l;
        double p_h;
        double p_l;
        double y1;
        double t1;
        double t2;
        double r;
        double s;
        double t;
        double u;
        double v;
        double w;
        int i;
        int j;
        int k;
        int yisint;
        int n;
        int hx;
        int hy;
        int ix;
        int iy;
        int lx;
        int ly;

        hx = __HI(x);
        lx = __LO(x);
        hy = __HI(y);
        ly = __LO(y);
        ix = hx & 0x7fffffff;
        iy = hy & 0x7fffffff;

        /* y==zero: x**0 = 1 */
        if ((iy | ly) == 0) {
            return one;
        }

        /* +-NaN return x+y */
        if ((ix > 0x7ff00000) || ((ix == 0x7ff00000) && (lx != 0)) ||
                (iy > 0x7ff00000) || ((iy == 0x7ff00000) && (ly != 0))) {
            return x + y;
        }

        /* determine if y is an odd int when x < 0
         * yisint = 0        ... y is not an integer
         * yisint = 1        ... y is an odd int
         * yisint = 2        ... y is an even int
         */
        yisint = 0;
        if (hx < 0) {
            if (iy >= 0x43400000) {
                yisint = 2; /* even integer y */
            } else if (iy >= 0x3ff00000) {
                k = (iy >> 20) - 0x3ff; /* exponent */
                if (k > 20) {
                    j = ly >>> (52 - k);
                    if ((j << (52 - k)) == ly) {
                        yisint = 2 - (j & 1);
                    }
                } else if (ly == 0) {
                    j = iy >> (20 - k);
                    if ((j << (20 - k)) == iy) {
                        yisint = 2 - (j & 1);
                    }
                }
            }
        }

        /* special value of y */
        if (ly == 0) {
            if (iy == 0x7ff00000) { /* y is +-inf */
                if (((ix - 0x3ff00000) | lx) == 0) {
                    return y - y; /* inf**+-1 is NaN */
                } else if (ix >= 0x3ff00000) { /* (|x|>1)**+-inf = inf,0 */
                    return (hy >= 0) ? y : zero;
                } else { /* (|x|<1)**-,+inf = inf,0 */
                    return (hy < 0) ? (-y) : zero;
                }
            }
            if (iy == 0x3ff00000) { /* y is  +-1 */
                if (hy < 0) {
                    return one / x;
                } else {
                    return x;
                }
            }
            if (hy == 0x40000000) {
                return x * x; /* y is  2 */
            }
            if (hy == 0x3fe00000) { /* y is  0.5 */
                if (hx >= 0) { /* x >= +0 */
                    return sqrt(x);
                }
            }
        }

        ax = abs(x);

        /* special value of x */
        if (lx == 0) {
            if ((ix == 0x7ff00000) || (ix == 0) || (ix == 0x3ff00000)) {
                z = ax; /*x is +-0,+-inf,+-1*/
                if (hy < 0) {
                    z = one / z; /* z = (1/|x|) */
                }
                if (hx < 0) {
                    if (((ix - 0x3ff00000) | yisint) == 0) {
                        z = (z - z) / (z - z); /* (-1)**non-int is NaN */
                    } else if (yisint == 1) {
                        z = -z; /* (x<0)**odd = -(|x|**odd) */
                    }
                }
                return z;
            }
        }

        /* (x<0)**(non-int) is NaN */
        if ((((hx >> 31) + 1) | yisint) == 0) {
            return (x - x) / (x - x);
        }

        /* |y| is huge */
        if (iy > 0x41e00000) { /* if |y| > 2**31 */
            if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */
                if (ix <= 0x3fefffff) {
                    return ((hy < 0) ? (huge * huge) : (tiny * tiny));
                }
                if (ix >= 0x3ff00000) {
                    return ((hy > 0) ? (huge * huge) : (tiny * tiny));
                }
            }

            /* over/underflow if x is not close to one */
            if (ix < 0x3fefffff) {
                return ((hy < 0) ? (huge * huge) : (tiny * tiny));
            }
            if (ix > 0x3ff00000) {
                return ((hy > 0) ? (huge * huge) : (tiny * tiny));
            }

            /* now |1-x| is tiny <= 2**-20, suffice to compute
               log(x) by x-x^2/2+x^3/3-x^4/4 */
            t = x - 1; /* t has 20 trailing zeros */
            w = (t * t) * (0.5 - (t * (0.3333333333333333333333 - (t * 0.25))));
            u = ivln2_h * t; /* ivln2_h has 21 sig. bits */
            v = (t * ivln2_l) - (w * ivln2);
            t1 = u + v;
            t1 = setLO(t1, 0);
            t2 = v - (t1 - u);
        } else {
            double s2;
            double s_h;
            double s_l;
            double t_h;
            double t_l;
            n = 0;

            /* take care subnormal number */
            if (ix < 0x00100000) {
                ax *= two53;
                n -= 53;
                ix = __HI(ax);
            }
            n += (((ix) >> 20) - 0x3ff);
            j = ix & 0x000fffff;

            /* determine interval */
            ix = j | 0x3ff00000; /* normalize ix */
            if (j <= 0x3988E) {
                k = 0; /* |x|<sqrt(3/2) */
            } else if (j < 0xBB67A) {
                k = 1; /* |x|<sqrt(3)   */
            } else {
                k = 0;
                n += 1;
                ix -= 0x00100000;
            }
            ax = setHI(ax, ix);

            /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
            u = ax - bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
            v = one / (ax + bp[k]);
            s = u * v;
            s_h = s;
            s_h = setLO(s_h, 0);

            /* t_h=ax+bp[k] High */
            t_h = zero;
            t_h = setHI(t_h, ((ix >> 1) | 0x20000000) + 0x00080000 + (k << 18));
            t_l = ax - (t_h - bp[k]);
            s_l = v * ((u - (s_h * t_h)) - (s_h * t_l));

            /* compute log(ax) */
            s2 = s * s;
            r = s2 * s2 * (L1 +
                (s2 * (L2 +
                (s2 * (L3 + (s2 * (L4 + (s2 * (L5 + (s2 * L6))))))))));
            r += (s_l * (s_h + s));
            s2 = s_h * s_h;
            t_h = 3.0 + s2 + r;
            t_h = setLO(t_h, 0);
            t_l = r - ((t_h - 3.0) - s2);

            /* u+v = s*(1+...) */
            u = s_h * t_h;
            v = (s_l * t_h) + (t_l * s);

            /* 2/(3log2)*(s+...) */
            p_h = u + v;
            p_h = setLO(p_h, 0);
            p_l = v - (p_h - u);
            z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */
            z_l = (cp_l * p_h) + (p_l * cp) + dp_l[k];

            /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
            t = (double) n;
            t1 = (((z_h + z_l) + dp_h[k]) + t);
            t1 = setLO(t1, 0);
            t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
        }

        s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
        if ((((hx >> 31) + 1) | (yisint - 1)) == 0) {
            s = -one; /* (-ve)**(odd int) */
        }

        /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
        y1 = y;
        y1 = setLO(y1, 0);
        p_l = ((y - y1) * t1) + (y * t2);
        p_h = y1 * t1;
        z = p_l + p_h;
        j = __HI(z);
        i = __LO(z);
        if (j >= 0x40900000) { /* z >= 1024 */
            if (((j - 0x40900000) | i) != 0) { /* if z > 1024 */
                return s * huge * huge; /* overflow */
            } else {
                if ((p_l + ovt) > (z - p_h)) {
                    return s * huge * huge; /* overflow */
                }
            }
        } else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */
            if (((j - 0xc090cc00) | i) != 0) { /* z < -1075 */
                return s * tiny * tiny; /* underflow */
            } else {
                if (p_l <= (z - p_h)) {
                    return s * tiny * tiny; /* underflow */
                }
            }
        }

        /*
         * compute 2**(p_h+p_l)
         */
        i = j & 0x7fffffff;
        k = (i >> 20) - 0x3ff;
        n = 0;
        if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
            n = j + (0x00100000 >> (k + 1));
            k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */
            t = zero;
            t = setHI(t, (n & ~(0x000fffff >> k)));
            n = ((n & 0x000fffff) | 0x00100000) >> (20 - k);
            if (j < 0) {
                n = -n;
            }
            p_h -= t;
        }
        t = p_l + p_h;
        t = setLO(t, 0);
        u = t * lg2_h;
        v = ((p_l - (t - p_h)) * lg2) + (t * lg2_l);
        z = u + v;
        w = v - (z - u);
        t = z * z;
        t1 = z - (t * (P1 + (t * (P2 + (t * (P3 + (t * (P4 + (t * P5)))))))));
        r = ((z * t1) / (t1 - 2.0)) - (w + (z * w));
        z = one - (r - z);
        j = __HI(z);
        j += (n << 20);
        if ((j >> 20) <= 0) {
            z = scalbn(z, n); /* subnormal output */
        } else {
            i = __HI(z);
            i += (n << 20);
            z = setHI(z, i);
        }
        return s * z;
    }

    /*
     * copysign(double x, double y)
     * copysign(x,y) returns a value with the magnitude of x and
     * with the sign bit of y.
     */
    static private double copysign(double x, double y) {
        long ix = Double.doubleToRawLongBits(x);
        long iy = Double.doubleToRawLongBits(y);
        ix = (0x7fffffffffffffffL & ix) | (0x8000000000000000L & iy);
        return Double.longBitsToDouble(ix);
    }

    static private final double two54 = 0x1.0p54; /*  1.80143985094819840000e+16 */
    static private final double twom54 = 0x1.0p-54; /*  5.55111512312578270212e-17 */

    /*
     * scalbn (double x, int n)
     * scalbn(x,n) returns x* 2**n  computed by  exponent
     * manipulation rather than by actually performing an
     * exponentiation or a multiplication.
     */
    static private double scalbn(double x, int n) {
        int k;
        int hx;
        int lx;
        hx = __HI(x);
        lx = __LO(x);
        k = (hx & 0x7ff00000) >> 20; /* extract exponent */
        if (k == 0) { /* 0 or subnormal x */
            if ((lx | (hx & 0x7fffffff)) == 0) {
                return x; /* +-0 */
            }
            x *= two54;
            hx = __HI(x);
            k = ((hx & 0x7ff00000) >> 20) - 54;
            if (n < -50000) {
                return tiny * x; /*underflow*/
            }
        }
        if (k == 0x7ff) {
            return x + x; /* NaN or Inf */
        }
        k = k + n;
        if (k > 0x7fe) {
            return huge * copysign(huge, x); /* overflow  */
        }
        if (k > 0) {

            /* normal result */
            return setHI(x, (hx & 0x800fffff) | (k << 20));
        }
        if (k <= -54) {
            if (n > 50000) { /* in case integer overflow in n+k */
                return huge * copysign(huge, x); /*overflow*/
            }
        } else {
            return tiny * copysign(tiny, x); /*underflow*/
        }
        k += 54; /* subnormal result */
        return twom54 * setHI(x, (hx & 0x800fffff) | (k << 20));
    }

    static private double set(int newHiPart, int newLowPart) {
        return Double.longBitsToDouble((((long) newHiPart) << 32) | newLowPart);
    }

    static private double setLO(double x, int newLowPart) {
        long lx = Double.doubleToRawLongBits(x);
        lx &= 0xFFFFFFFF00000000L;
        lx |= newLowPart;
        return Double.longBitsToDouble(lx);
    }

    static private double setHI(double x, int newHiPart) {
        long lx = Double.doubleToRawLongBits(x);
        lx &= 0x00000000FFFFFFFFL;
        lx |= (((long) newHiPart) << 32);
        return Double.longBitsToDouble(lx);
    }

    static private int __HI(double x) {
        return (int) (0xFFFFFFFF & (Double.doubleToRawLongBits(x) >> 32));
    }

    static private int __LO(double x) {
        return (int) (0xFFFFFFFF & Double.doubleToRawLongBits(x));
    }
} 