public static double pow()

in commons-math-core/src/main/java/org/apache/commons/math4/core/jdkmath/AccurateMath.java [1442:1572]


    public static double pow(final double x, final double y) {

        if (y == 0) {
            // y = -0 or y = +0
            return 1.0;
        } else {

            final long yBits        = Double.doubleToRawLongBits(y);
            final int  yRawExp      = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
            final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
            final long xBits        = Double.doubleToRawLongBits(x);
            final int  xRawExp      = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
            final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;

            if (yRawExp > 1085) {
                // y is either a very large integral value that does not fit in a long or it is a special number

                if ((yRawExp == 2047 && yRawMantissa != 0) ||
                    (xRawExp == 2047 && xRawMantissa != 0)) {
                    // NaN
                    return Double.NaN;
                } else if (xRawExp == 1023 && xRawMantissa == 0) {
                    // x = -1.0 or x = +1.0
                    if (yRawExp == 2047) {
                        // y is infinite
                        return Double.NaN;
                    } else {
                        // y is a large even integer
                        return 1.0;
                    }
                } else {
                    // the absolute value of x is either greater or smaller than 1.0

                    // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
                    // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
                    // accuracy, at this magnitude it behaves just like infinity with regards to x
                    if ((y > 0) ^ (xRawExp < 1023)) {
                        // either y = +infinity (or large engouh) and abs(x) > 1.0
                        // or     y = -infinity (or large engouh) and abs(x) < 1.0
                        return Double.POSITIVE_INFINITY;
                    } else {
                        // either y = +infinity (or large engouh) and abs(x) < 1.0
                        // or     y = -infinity (or large engouh) and abs(x) > 1.0
                        return +0.0;
                    }
                }
            } else {
                // y is a regular non-zero number

                if (yRawExp >= 1023) {
                    // y may be an integral value, which should be handled specifically
                    final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
                    if (yRawExp < 1075) {
                        // normal number with negative shift that may have a fractional part
                        final long integralMask = (-1L) << (1075 - yRawExp);
                        if ((yFullMantissa & integralMask) == yFullMantissa) {
                            // all fractional bits are 0, the number is really integral
                            final long l = yFullMantissa >> (1075 - yRawExp);
                            return AccurateMath.pow(x, (y < 0) ? -l : l);
                        }
                    } else {
                        // normal number with positive shift, always an integral value
                        // we know it fits in a primitive long because yRawExp > 1085 has been handled above
                        final long l =  yFullMantissa << (yRawExp - 1075);
                        return AccurateMath.pow(x, (y < 0) ? -l : l);
                    }
                }

                // y is a non-integral value

                if (x == 0) {
                    // x = -0 or x = +0
                    // the integer powers have already been handled above
                    return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
                } else if (xRawExp == 2047) {
                    if (xRawMantissa == 0) {
                        // x = -infinity or x = +infinity
                        return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
                    } else {
                        // NaN
                        return Double.NaN;
                    }
                } else if (x < 0) {
                    // the integer powers have already been handled above
                    return Double.NaN;
                } else {

                    // this is the general case, for regular fractional numbers x and y

                    // Split y into ya and yb such that y = ya+yb
                    final double tmp = y * HEX_40000000;
                    final double ya = (y + tmp) - tmp;
                    final double yb = y - ya;

                    /* Compute ln(x) */
                    final double[] lns = new double[2];
                    final double lores = log(x, lns);
                    if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
                        return lores;
                    }

                    double lna = lns[0];
                    double lnb = lns[1];

                    /* resplit lns */
                    final double tmp1 = lna * HEX_40000000;
                    final double tmp2 = (lna + tmp1) - tmp1;
                    lnb += lna - tmp2;
                    lna = tmp2;

                    // y*ln(x) = (aa+ab)
                    final double aa = lna * ya;
                    final double ab = lna * yb + lnb * ya + lnb * yb;

                    lna = aa + ab;
                    lnb = -(lna - aa - ab);

                    double z = 1.0 / 120.0;
                    z = z * lnb + (1.0 / 24.0);
                    z = z * lnb + (1.0 / 6.0);
                    z = z * lnb + 0.5;
                    z = z * lnb + 1.0;
                    z *= lnb;

                    final double result = exp(lna, z, null);
                    //result = result + result * z;
                    return result;
                }
            }
        }
    }