cc.redberry.core.number.Exponentiation.java Source code

Java tutorial

Introduction

Here is the source code for cc.redberry.core.number.Exponentiation.java

Source

/*
 * Redberry: symbolic tensor computations.
 *
 * Copyright (c) 2010-2013:
 *   Stanislav Poslavsky   <stvlpos@mail.ru>
 *   Bolotin Dmitriy       <bolotin.dmitriy@gmail.com>
 *
 * This file is part of Redberry.
 *
 * Redberry is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Redberry is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Redberry. If not, see <http://www.gnu.org/licenses/>.
 */
package cc.redberry.core.number;

import org.apache.commons.math3.util.ArithmeticUtils;

import java.math.BigInteger;

/**
 * This class contains methods for complex numbers exponentiation.
 */
public final class Exponentiation {
    public static Real exponentiateIfPossible(Real base, Real power) {
        if (base.isZero())
            if (power.isInfinite())
                return Numeric.NaN;
            else
                return Rational.ZERO;

        if (base.isNumeric() || power.isNumeric()) { // Bease or power are numeric
            return new Numeric(
                    Math.pow(base.getNumericValue().doubleValue(), power.getNumericValue().doubleValue()));
        }

        //<-- Power and Base are rational

        if (power.isInteger())
            return new Rational(((Rational) base).getBigFraction().pow(((Rational) power).getNumerator())); //Using BigFraction pow method.

        //<-- Power is not integer

        BigInteger powerNum = ((Rational) power).getNumerator();
        BigInteger powerDen = ((Rational) power).getDenominator();

        BigInteger baseNum = ((Rational) base).getNumerator();
        BigInteger baseDen = ((Rational) base).getDenominator();

        baseNum = findIntegerRoot(baseNum, powerDen);
        baseDen = findIntegerRoot(baseDen, powerDen);

        if (baseNum == null || baseDen == null) //Result is irrational
            return null;

        return exponentiateIfPossible(new Rational(baseNum, baseDen), new Rational(powerNum));
    }

    private static BigInteger BI_MINUS_ONE = BigInteger.ONE.negate();

    static BigInteger findIntegerRoot(BigInteger base, BigInteger power) {
        BigInteger maxBits = BigInteger.valueOf(base.bitLength() + 1); // base < 2 ^ (maxBits + 1)
        // => base ^ ( 1 / power ) < 2 ^ ( (maxBits + 1) / power )

        BigInteger[] divResult = maxBits.divideAndRemainder(power);
        if (divResult[1].signum() == 0) // i.e. divResult[1] == 0
            maxBits = divResult[0];
        else
            maxBits = divResult[0].add(BigInteger.ONE);

        if (maxBits.bitLength() > 31)
            throw new RuntimeException("Too many bits...");

        int targetBitsNumber = maxBits.intValue();
        int resultLengthM1 = targetBitsNumber / 8 + 1; //resultLength minus one
        byte[] result = new byte[resultLengthM1];
        resultLengthM1--;

        int bitNumber = targetBitsNumber;

        int cValue;
        BigInteger testValue;

        while ((--bitNumber) >= 0) {
            //setting bit
            result[resultLengthM1 - (bitNumber >> 3)] |= 1 << (bitNumber & 0x7);

            //Testing
            testValue = new BigInteger(result);
            cValue = ArithmeticUtils.pow(testValue, power).compareTo(base);
            if (cValue == 0)
                return testValue;
            if (cValue > 0)
                result[resultLengthM1 - (bitNumber >> 3)] &= ~(1 << (bitNumber & 0x7));
        }

        return null;
    }

    public static Complex exponentiateIfPossible(Complex base, Complex power) {
        //Partially copied from PowerFactory
        if (base.isInfinite())
            if (power.isZero())
                return Complex.ComplexNaN;
            else
                return base;

        if (base.isOne())
            if (power.isInfinite())
                return power.multiply(base);
            else
                return base;

        if (power.isOne())
            return base;

        if (base.isZero()) {
            if (power.getReal().signum() <= 0)
                return Complex.ComplexNaN;
            return base;
        }

        if (power.isZero())
            return Complex.ONE;

        if (base.isNumeric() || power.isNumeric())
            return base.powNumeric(power);

        // <-- base and power are rational

        if (power.isReal()) {
            Rational pp = (Rational) power.getReal();

            if (base.isReal()) {
                Real value = exponentiateIfPossible(base.getReal(), pp);
                if (value == null)
                    return null;
                return new Complex(value);
            }

            if (pp.isInteger())
                return base.pow(pp.getNumerator());
            else {
                Complex root = findIntegerRoot(base, pp.getDenominator());
                if (root == null)
                    return null;
                return root.pow(pp.getNumerator());
            }
        }

        return null;
    }

    public static Complex findIntegerRoot(Complex base, BigInteger power) {
        BigInteger rDenominator = ((Rational) base.getReal()).getDenominator();
        BigInteger iDenominator = ((Rational) base.getImaginary()).getDenominator();

        BigInteger lcm = rDenominator.gcd(iDenominator);
        lcm = rDenominator.divide(lcm);
        lcm = lcm.multiply(iDenominator);

        BigInteger lcmRoot = findIntegerRoot(lcm, power);

        if (lcm == null)
            return null;

        base = base.multiply(lcm);

        Complex numericValue = base.pow(1.0 / power.doubleValue());
        double real = numericValue.getReal().doubleValue();
        double imaginary = numericValue.getImaginary().doubleValue();

        int ceilReal = (int) Math.ceil(real), floorReal = (int) Math.floor(real),
                ceilImaginary = (int) Math.ceil(imaginary), floorImaginary = (int) Math.floor(imaginary);

        Complex candidate;
        if ((candidate = new Complex(ceilReal, ceilImaginary)).pow(power).equals(base))
            return candidate.divide(lcmRoot);
        if ((candidate = new Complex(floorReal, ceilImaginary)).pow(power).equals(base))
            return candidate.divide(lcmRoot);
        if ((candidate = new Complex(ceilReal, floorImaginary)).pow(power).equals(base))
            return candidate.divide(lcmRoot);
        if ((candidate = new Complex(floorReal, floorImaginary)).pow(power).equals(base))
            return candidate.divide(lcmRoot);
        return null;
    }
}