android/calculator2/UnifiedReal.java

/*
 * Copyright (C) 2016 The Android Open Source Project
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package com.android.calculator2;

import java.math.BigInteger;
import com.hp.creals.CR;
import com.hp.creals.UnaryCRFunction;

/**
 * Computable real numbers, represented so that we can get exact decidable comparisons
 * for a number of interesting special cases, including rational computations.
 *
 * A real number is represented as the product of two numbers with different representations:
 * A) A BoundedRational that can only represent a subset of the rationals, but supports
 *    exact computable comparisons.
 * B) A lazily evaluated "constructive real number" that provides operations to evaluate
 *    itself to any requested number of digits.
 * Whenever possible, we choose (B) to be one of a small set of known constants about which we
 * know more.  For example, whenever we can, we represent rationals such that (B) is 1.
 * This scheme allows us to do some very limited symbolic computation on numbers when both
 * have the same (B) value, as well as in some other situations.  We try to maximize that
 * possibility.
 *
 * Arithmetic operations and operations that produce finite approximations may throw unchecked
 * exceptions produced by the underlying CR and BoundedRational packages, including
 * CR.PrecisionOverflowException and CR.AbortedException.
 */
public class UnifiedReal {

    private final BoundedRational mRatFactor;
    private final CR mCrFactor;
    // TODO: It would be helpful to add flags to indicate whether the result is known
    // irrational, etc.  This sometimes happens even if mCrFactor is not one of the known ones.
    // And exact comparisons between rationals and known irrationals are decidable.

    /**
     * Perform some nontrivial consistency checks.
     * @hide
     */
    public static boolean enableChecks = true;

    private static void check(boolean b) {
        if (!b) {
            throw new AssertionError();
        }
    }

    private UnifiedReal(BoundedRational rat, CR cr) {
        if (rat == null) {
            throw new ArithmeticException("Building UnifiedReal from null");
        }
        // We don't normally traffic in null CRs, and hence don't test explicitly.
        mCrFactor = cr;
        mRatFactor = rat;
    }

    public UnifiedReal(CR cr) {
        this(BoundedRational.ONE, cr);
    }

    public UnifiedReal(BoundedRational rat) {
        this(rat, CR_ONE);
    }

    public UnifiedReal(BigInteger n) {
        this(new BoundedRational(n));
    }

    public UnifiedReal(long n) {
        this(new BoundedRational(n));
    }

    public static UnifiedReal valueOf(double x) {
        if (x == 0.0 || x == 1.0) {
            return valueOf((long) x);
        }
        return new UnifiedReal(BoundedRational.valueOf(x));
    }

    public static UnifiedReal valueOf(long x) {
        if (x == 0) {
            return UnifiedReal.ZERO;
        } else if (x == 1) {
            return UnifiedReal.ONE;
        } else {
            return new UnifiedReal(BoundedRational.valueOf(x));
        }
    }

    // Various helpful constants
    private final static BigInteger BIG_24 = BigInteger.valueOf(24);
    private final static int DEFAULT_COMPARE_TOLERANCE = -1000;

    // Well-known CR constants we try to use in the mCrFactor position:
    private final static CR CR_ONE = CR.ONE;
    private final static CR CR_PI = CR.PI;
    private final static CR CR_E = CR.ONE.exp();
    private final static CR CR_SQRT2 = CR.valueOf(2).sqrt();
    private final static CR CR_SQRT3 = CR.valueOf(3).sqrt();
    private final static CR CR_LN2 = CR.valueOf(2).ln();
    private final static CR CR_LN3 = CR.valueOf(3).ln();
    private final static CR CR_LN5 = CR.valueOf(5).ln();
    private final static CR CR_LN6 = CR.valueOf(6).ln();
    private final static CR CR_LN7 = CR.valueOf(7).ln();
    private final static CR CR_LN10 = CR.valueOf(10).ln();

    // Square roots that we try to recognize.
    // We currently recognize only a small fixed collection, since the sqrt() function needs to
    // identify numbers of the form <SQRT[i]>*n^2, and we don't otherwise know of a good
    // algorithm for that.
    private final static CR sSqrts[] = {
            null,
            CR.ONE,
            CR_SQRT2,
            CR_SQRT3,
            null,
            CR.valueOf(5).sqrt(),
            CR.valueOf(6).sqrt(),
            CR.valueOf(7).sqrt(),
            null,
            null,
            CR.valueOf(10).sqrt() };

    // Natural logs of small integers that we try to recognize.
    private final static CR sLogs[] = {
            null,
            null,
            CR_LN2,
            CR_LN3,
            null,
            CR_LN5,
            CR_LN6,
            CR_LN7,
            null,
            null,
            CR_LN10 };


    // Some convenient UnifiedReal constants.
    public static final UnifiedReal PI = new UnifiedReal(CR_PI);
    public static final UnifiedReal E = new UnifiedReal(CR_E);
    public static final UnifiedReal ZERO = new UnifiedReal(BoundedRational.ZERO);
    public static final UnifiedReal ONE = new UnifiedReal(BoundedRational.ONE);
    public static final UnifiedReal MINUS_ONE = new UnifiedReal(BoundedRational.MINUS_ONE);
    public static final UnifiedReal TWO = new UnifiedReal(BoundedRational.TWO);
    public static final UnifiedReal MINUS_TWO = new UnifiedReal(BoundedRational.MINUS_TWO);
    public static final UnifiedReal HALF = new UnifiedReal(BoundedRational.HALF);
    public static final UnifiedReal MINUS_HALF = new UnifiedReal(BoundedRational.MINUS_HALF);
    public static final UnifiedReal TEN = new UnifiedReal(BoundedRational.TEN);
    public static final UnifiedReal RADIANS_PER_DEGREE
            = new UnifiedReal(new BoundedRational(1, 180), CR_PI);
    private static final UnifiedReal SIX = new UnifiedReal(6);
    private static final UnifiedReal HALF_SQRT2 = new UnifiedReal(BoundedRational.HALF, CR_SQRT2);
    private static final UnifiedReal SQRT3 = new UnifiedReal(CR_SQRT3);
    private static final UnifiedReal HALF_SQRT3 = new UnifiedReal(BoundedRational.HALF, CR_SQRT3);
    private static final UnifiedReal THIRD_SQRT3 = new UnifiedReal(BoundedRational.THIRD, CR_SQRT3);
    private static final UnifiedReal PI_OVER_2 = new UnifiedReal(BoundedRational.HALF, CR_PI);
    private static final UnifiedReal PI_OVER_3 = new UnifiedReal(BoundedRational.THIRD, CR_PI);
    private static final UnifiedReal PI_OVER_4 = new UnifiedReal(BoundedRational.QUARTER, CR_PI);
    private static final UnifiedReal PI_OVER_6 = new UnifiedReal(BoundedRational.SIXTH, CR_PI);


    /**
     * Given a constructive real cr, try to determine whether cr is the square root of
     * a small integer.  If so, return its square as a BoundedRational.  Otherwise return null.
     * We make this determination by simple table lookup, so spurious null returns are
     * entirely possible, or even likely.
     */
    private static BoundedRational getSquare(CR cr) {
        for (int i = 0; i < sSqrts.length; ++i) {
             if (sSqrts[i] == cr) {
                return new BoundedRational(i);
             }
        }
        return null;
    }

    /**
     * Given a constructive real cr, try to determine whether cr is the logarithm of a small
     * integer.  If so, return exp(cr) as a BoundedRational.  Otherwise return null.
     * We make this determination by simple table lookup, so spurious null returns are
     * entirely possible, or even likely.
     */
    private BoundedRational getExp(CR cr) {
        for (int i = 0; i < sLogs.length; ++i) {
             if (sLogs[i] == cr) {
                return new BoundedRational(i);
             }
        }
        return null;
    }

    /**
     * If the argument is a well-known constructive real, return its name.
     * The name of "CR_ONE" is the empty string.
     * No named constructive reals are rational multiples of each other.
     * Thus two UnifiedReals with different named mCrFactors can be equal only if both
     * mRatFactors are zero or possibly if one is CR_PI and the other is CR_E.
     * (The latter is apparently an open problem.)
     */
    private static String crName(CR cr) {
        if (cr == CR_ONE) {
            return "";
        }
        if (cr == CR_PI) {
            return "\u03C0";   // GREEK SMALL LETTER PI
        }
        if (cr == CR_E) {
            return "e";
        }
        for (int i = 0; i < sSqrts.length; ++i) {
            if (cr == sSqrts[i]) {
                return "\u221A" /* SQUARE ROOT */ + i;
            }
        }
        for (int i = 0; i < sLogs.length; ++i) {
            if (cr == sLogs[i]) {
                return "ln(" + i + ")";
            }
        }
        return null;
    }

    /**
     * Would crName() return non-Null?
     */
    private static boolean isNamed(CR cr) {
        if (cr == CR_ONE || cr == CR_PI || cr == CR_E) {
            return true;
        }
        for (CR r: sSqrts) {
            if (cr == r) {
                return true;
            }
        }
        for (CR r: sLogs) {
            if (cr == r) {
                return true;
            }
        }
        return false;
    }

    /**
     * Is cr known to be algebraic (as opposed to transcendental)?
     * Currently only produces meaningful results for the above known special
     * constructive reals.
     */
    private static boolean definitelyAlgebraic(CR cr) {
        return cr == CR_ONE || getSquare(cr) != null;
    }

    /**
     * Is this number known to be rational?
     */
    public boolean definitelyRational() {
        return mCrFactor == CR_ONE || mRatFactor.signum() == 0;
    }

    /**
     * Is this number known to be irrational?
     * TODO: We could track the fact that something is irrational with an explicit flag, which
     * could cover many more cases.  Whether that matters in practice is TBD.
     */
    public boolean definitelyIrrational() {
        return !definitelyRational() && isNamed(mCrFactor);
    }

    /**
     * Is this number known to be algebraic?
     */
    public boolean definitelyAlgebraic() {
        return definitelyAlgebraic(mCrFactor) || mRatFactor.signum() == 0;
    }

    /**
     * Is this number known to be transcendental?
     */
    public boolean definitelyTranscendental() {
        return !definitelyAlgebraic() && isNamed(mCrFactor);
    }


    /**
     * Is it known that the two constructive reals differ by something other than a
     * a rational factor, i.e. is it known that two UnifiedReals
     * with those mCrFactors will compare unequal unless both mRatFactors are zero?
     * If this returns true, then a comparison of two UnifiedReals using those two
     * mCrFactors cannot diverge, though we don't know of a good runtime bound.
     */
    private static boolean definitelyIndependent(CR r1, CR r2) {
        // The question here is whether r1 = x*r2, where x is rational, where r1 and r2
        // are in our set of special known CRs, can have a solution.
        // This cannot happen if one is CR_ONE and the other is not.
        // (Since all others are irrational.)
        // This cannot happen for two named square roots, which have no repeated factors.
        // (To see this, square both sides of the equation and factor.  Each prime
        // factor in the numerator and denominator occurs twice.)
        // This cannot happen for e or pi on one side, and a square root on the other.
        // (One is transcendental, the other is algebraic.)
        // This cannot happen for two of our special natural logs.
        // (Otherwise ln(m) = (a/b)ln(n) ==> m = n^(a/b) ==> m^b = n^a, which is impossible
        // because either m or n includes a prime factor not shared by the other.)
        // This cannot happen for a log and a square root.
        // (The Lindemann-Weierstrass theorem tells us, among other things, that if
        // a is algebraic, then exp(a) is transcendental.  Thus if l in our finite
        // set of logs where algebraic, expl(l), must be transacendental.
        // But exp(l) is an integer.  Thus the logs are transcendental.  But of course the
        // square roots are algebraic.  Thus they can't be rational multiples.)
        // Unfortunately, we do not know whether e/pi is rational.
        if (r1 == r2) {
            return false;
        }
        CR other;
        if (r1 == CR_E || r1 == CR_PI) {
            return definitelyAlgebraic(r2);
        }
        if (r2 == CR_E || r2 == CR_PI) {
            return definitelyAlgebraic(r1);
        }
        return isNamed(r1) && isNamed(r2);
    }

    /**
     * Convert to String reflecting raw representation.
     * Debug or log messages only, not pretty.
     */
    public String toString() {
        return mRatFactor.toString() + "*" + mCrFactor.toString();
    }

    /**
     * Convert to readable String.
     * Intended for user output.  Produces exact expression when possible.
     */
    public String toNiceString() {
        if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
            return mRatFactor.toNiceString();
        }
        String name = crName(mCrFactor);
        if (name != null) {
            BigInteger bi = BoundedRational.asBigInteger(mRatFactor);
            if (bi != null) {
                if (bi.equals(BigInteger.ONE)) {
                    return name;
                }
                return mRatFactor.toNiceString() + name;
            }
            return "(" + mRatFactor.toNiceString() + ")" + name;
        }
        if (mRatFactor.equals(BoundedRational.ONE)) {
            return mCrFactor.toString();
        }
        return crValue().toString();
    }

    /**
     * Will toNiceString() produce an exact representation?
     */
    public boolean exactlyDisplayable() {
        return crName(mCrFactor) != null;
    }

    // Number of extra bits used in evaluation below to prefer truncation to rounding.
    // Must be <= 30.
    private final static int EXTRA_PREC = 10;

    /*
     * Returns a truncated representation of the result.
     * If exactlyTruncatable(), we round correctly towards zero. Otherwise the resulting digit
     * string may occasionally be rounded up instead.
     * Always includes a decimal point in the result.
     * The result includes n digits to the right of the decimal point.
     * @param n result precision, >= 0
     */
    public String toStringTruncated(int n) {
        if (mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO) {
            return mRatFactor.toStringTruncated(n);
        }
        final CR scaled = CR.valueOf(BigInteger.TEN.pow(n)).multiply(crValue());
        boolean negative = false;
        BigInteger intScaled;
        if (exactlyTruncatable()) {
            intScaled = scaled.get_appr(0);
            if (intScaled.signum() < 0) {
                negative = true;
                intScaled = intScaled.negate();
            }
            if (CR.valueOf(intScaled).compareTo(scaled.abs()) > 0) {
                intScaled = intScaled.subtract(BigInteger.ONE);
            }
            check(CR.valueOf(intScaled).compareTo(scaled.abs()) < 0);
        } else {
            // Approximate case.  Exact comparisons are impossible.
            intScaled = scaled.get_appr(-EXTRA_PREC);
            if (intScaled.signum() < 0) {
                negative = true;
                intScaled = intScaled.negate();
            }
            intScaled = intScaled.shiftRight(EXTRA_PREC);
        }
        String digits = intScaled.toString();
        int len = digits.length();
        if (len < n + 1) {
            digits = StringUtils.repeat('0', n + 1 - len) + digits;
            len = n + 1;
        }
        return (negative ? "-" : "") + digits.substring(0, len - n) + "."
                + digits.substring(len - n);
    }

    /*
     * Can we compute correctly truncated approximations of this number?
     */
    public boolean exactlyTruncatable() {
        // If the value is known rational, we can do exact comparisons.
        // If the value is known irrational, then we can safely compare to rational approximations;
        // equality is impossible; hence the comparison must converge.
        // The only problem cases are the ones in which we don't know.
        return mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO || definitelyIrrational();
    }

    /**
     * Return a double approximation.
     * Rational arguments are currently rounded to nearest, with ties away from zero.
     * TODO: Improve rounding.
     */
    public double doubleValue() {
        if (mCrFactor == CR_ONE) {
            return mRatFactor.doubleValue(); // Hopefully correctly rounded
        } else {
            return crValue().doubleValue(); // Approximately correctly rounded
        }
    }

    public CR crValue() {
        return mRatFactor.crValue().multiply(mCrFactor);
    }

    /**
     * Are this and r exactly comparable?
     */
    public boolean isComparable(UnifiedReal u) {
        // We check for ONE only to speed up the common case.
        // The use of a tolerance here means we can spuriously return false, not true.
        return mCrFactor == u.mCrFactor
                && (isNamed(mCrFactor) || mCrFactor.signum(DEFAULT_COMPARE_TOLERANCE) != 0)
                || mRatFactor.signum() == 0 && u.mRatFactor.signum() == 0
                || definitelyIndependent(mCrFactor, u.mCrFactor)
                || crValue().compareTo(u.crValue(), DEFAULT_COMPARE_TOLERANCE) != 0;
    }

    /**
     * Return +1 if this is greater than r, -1 if this is less than r, or 0 of the two are
     * known to be equal.
     * May diverge if the two are equal and !isComparable(r).
     */
    public int compareTo(UnifiedReal u) {
        if (definitelyZero() && u.definitelyZero()) return 0;
        if (mCrFactor == u.mCrFactor) {
            int signum = mCrFactor.signum();  // Can diverge if mCRFactor == 0.
            return signum * mRatFactor.compareTo(u.mRatFactor);
        }
        return crValue().compareTo(u.crValue());  // Can also diverge.
    }

    /**
     * Return +1 if this is greater than r, -1 if this is less than r, or possibly 0 of the two are
     * within 2^a of each other.
     */
    public int compareTo(UnifiedReal u, int a) {
        if (isComparable(u)) {
            return compareTo(u);
        } else {
            return crValue().compareTo(u.crValue(), a);
        }
    }

    /**
     * Return compareTo(ZERO, a).
     */
    public int signum(int a) {
        return compareTo(ZERO, a);
    }

    /**
     * Return compareTo(ZERO).
     * May diverge for ZERO argument if !isComparable(ZERO).
     */
    public int signum() {
        return compareTo(ZERO);
    }

    /**
     * Equality comparison.  May erroneously return true if values differ by less than 2^a,
     * and !isComparable(u).
     */
    public boolean approxEquals(UnifiedReal u, int a) {
        if (isComparable(u)) {
            if (definitelyIndependent(mCrFactor, u.mCrFactor)
                    && (mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0)) {
                // No need to actually evaluate, though we don't know which is larger.
                return false;
            } else {
                return compareTo(u) == 0;
            }
        }
        return crValue().compareTo(u.crValue(), a) == 0;
    }

    /**
     * Returns true if values are definitely known to be equal, false in all other cases.
     * This does not satisfy the contract for Object.equals().
     */
    public boolean definitelyEquals(UnifiedReal u) {
        return isComparable(u) && compareTo(u) == 0;
    }

    @Override
    public int hashCode() {
        // Better useless than wrong. Probably.
        return 0;
    }

    @Override
    public boolean equals(Object r) {
        if (r == null || !(r instanceof UnifiedReal)) {
            return false;
        }
        // This is almost certainly a programming error. Don't even try.
        throw new AssertionError("Can't compare UnifiedReals for exact equality");
    }

    /**
     * Returns true if values are definitely known not to be equal, false in all other cases.
     * Performs no approximate evaluation.
     */
    public boolean definitelyNotEquals(UnifiedReal u) {
        boolean isNamed = isNamed(mCrFactor);
        boolean uIsNamed = isNamed(u.mCrFactor);
        if (isNamed && uIsNamed) {
            if (definitelyIndependent(mCrFactor, u.mCrFactor)) {
                return mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0;
            } else if (mCrFactor == u.mCrFactor) {
                return !mRatFactor.equals(u.mRatFactor);
            }
            return !mRatFactor.equals(u.mRatFactor);
        }
        if (mRatFactor.signum() == 0) {
            return uIsNamed && u.mRatFactor.signum() != 0;
        }
        if (u.mRatFactor.signum() == 0) {
            return isNamed && mRatFactor.signum() != 0;
        }
        return false;
    }

    // And some slightly faster convenience functions for special cases:

    public boolean definitelyZero() {
        return mRatFactor.signum() == 0;
    }

    /**
     * Can this number be determined to be definitely nonzero without performing approximate
     * evaluation?
     */
    public boolean definitelyNonZero() {
        return isNamed(mCrFactor) && mRatFactor.signum() != 0;
    }

    public boolean definitelyOne() {
        return mCrFactor == CR_ONE && mRatFactor.equals(BoundedRational.ONE);
    }

    /**
     * Return equivalent BoundedRational, if known to exist, null otherwise
     */
    public BoundedRational boundedRationalValue() {
        if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
            return mRatFactor;
        }
        return null;
    }

    /**
     * Returns equivalent BigInteger result if it exists, null if not.
     */
    public BigInteger bigIntegerValue() {
        final BoundedRational r = boundedRationalValue();
        return BoundedRational.asBigInteger(r);
    }

    public UnifiedReal add(UnifiedReal u) {
        if (mCrFactor == u.mCrFactor) {
            BoundedRational nRatFactor = BoundedRational.add(mRatFactor, u.mRatFactor);
            if (nRatFactor != null) {
                return new UnifiedReal(nRatFactor, mCrFactor);
            }
        }
        if (definitelyZero()) {
            // Avoid creating new mCrFactor, even if they don't currently match.
            return u;
        }
        if (u.definitelyZero()) {
            return this;
        }
        return new UnifiedReal(crValue().add(u.crValue()));
    }

    public UnifiedReal negate() {
        return new UnifiedReal(BoundedRational.negate(mRatFactor), mCrFactor);
    }

    public UnifiedReal subtract(UnifiedReal u) {
        return add(u.negate());
    }

    public UnifiedReal multiply(UnifiedReal u) {
        // Preserve a preexisting mCrFactor when we can.
        if (mCrFactor == CR_ONE) {
            BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
            if (nRatFactor != null) {
                return new UnifiedReal(nRatFactor, u.mCrFactor);
            }
        }
        if (u.mCrFactor == CR_ONE) {
            BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
            if (nRatFactor != null) {
                return new UnifiedReal(nRatFactor, mCrFactor);
            }
        }
        if (definitelyZero() || u.definitelyZero()) {
            return ZERO;
        }
        if (mCrFactor == u.mCrFactor) {
            BoundedRational square = getSquare(mCrFactor);
            if (square != null) {
                BoundedRational nRatFactor = BoundedRational.multiply(
                        BoundedRational.multiply(square, mRatFactor), u.mRatFactor);
                if (nRatFactor != null) {
                    return new UnifiedReal(nRatFactor);
                }
            }
        }
        // Probably a bit cheaper to multiply component-wise.
        BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
        if (nRatFactor != null) {
            return new UnifiedReal(nRatFactor, mCrFactor.multiply(u.mCrFactor));
        }
        return new UnifiedReal(crValue().multiply(u.crValue()));
    }

    public static class ZeroDivisionException extends ArithmeticException {
        public ZeroDivisionException() {
            super("Division by zero");
        }
    }

    /**
     * Return the reciprocal.
     */
    public UnifiedReal inverse() {
        if (definitelyZero()) {
            throw new ZeroDivisionException();
        }
        BoundedRational square = getSquare(mCrFactor);
        if (square != null) {
            // 1/sqrt(n) = sqrt(n)/n
            BoundedRational nRatFactor = BoundedRational.inverse(
                    BoundedRational.multiply(mRatFactor, square));
            if (nRatFactor != null) {
                return new UnifiedReal(nRatFactor, mCrFactor);
            }
        }
        return new UnifiedReal(BoundedRational.inverse(mRatFactor), mCrFactor.inverse());
    }

    public UnifiedReal divide(UnifiedReal u) {
        if (mCrFactor == u.mCrFactor) {
            if (u.definitelyZero()) {
                throw new ZeroDivisionException();
            }
            BoundedRational nRatFactor = BoundedRational.divide(mRatFactor, u.mRatFactor);
            if (nRatFactor != null) {
                return new UnifiedReal(nRatFactor, CR_ONE);
            }
        }
        return multiply(u.inverse());
    }

    /**
     * Return the square root.
     * This may fail to return a known rational value, even when the result is rational.
     */
    public UnifiedReal sqrt() {
        if (definitelyZero()) {
            return ZERO;
        }
        if (mCrFactor == CR_ONE) {
            BoundedRational ratSqrt;
            // Check for all arguments of the form <perfect rational square> * small_int,
            // where small_int has a known sqrt.  This includes the small_int = 1 case.
            for (int divisor = 1; divisor < sSqrts.length; ++divisor) {
                if (sSqrts[divisor] != null) {
                    ratSqrt = BoundedRational.sqrt(
                            BoundedRational.divide(mRatFactor, new BoundedRational(divisor)));
                    if (ratSqrt != null) {
                        return new UnifiedReal(ratSqrt, sSqrts[divisor]);
                    }
                }
            }
        }
        return new UnifiedReal(crValue().sqrt());
    }

    /**
     * Return (this mod 2pi)/(pi/6) as a BigInteger, or null if that isn't easily possible.
     */
    private BigInteger getPiTwelfths() {
        if (definitelyZero()) return BigInteger.ZERO;
        if (mCrFactor == CR_PI) {
            BigInteger quotient = BoundedRational.asBigInteger(
                    BoundedRational.multiply(mRatFactor, BoundedRational.TWELVE));
            if (quotient == null) {
                return null;
            }
            return quotient.mod(BIG_24);
        }
        return null;
    }

    /**
     * Computer the sin() for an integer multiple n of pi/12, if easily representable.
     * @param n value between 0 and 23 inclusive.
     */
    private static UnifiedReal sinPiTwelfths(int n) {
        if (n >= 12) {
            UnifiedReal negResult = sinPiTwelfths(n - 12);
            return negResult == null ? null : negResult.negate();
        }
        switch (n) {
        case 0:
            return ZERO;
        case 2: // 30 degrees
            return HALF;
        case 3: // 45 degrees
            return HALF_SQRT2;
        case 4: // 60 degrees
            return HALF_SQRT3;
        case 6:
            return ONE;
        case 8:
            return HALF_SQRT3;
        case 9:
            return HALF_SQRT2;
        case 10:
            return HALF;
        default:
            return null;
        }
    }

    public UnifiedReal sin() {
        BigInteger piTwelfths = getPiTwelfths();
        if (piTwelfths != null) {
            UnifiedReal result = sinPiTwelfths(piTwelfths.intValue());
            if (result != null) {
                return result;
            }
        }
        return new UnifiedReal(crValue().sin());
    }

    private static UnifiedReal cosPiTwelfths(int n) {
        int sinArg = n + 6;
        if (sinArg >= 24) {
            sinArg -= 24;
        }
        return sinPiTwelfths(sinArg);
    }

    public UnifiedReal cos() {
        BigInteger piTwelfths = getPiTwelfths();
        if (piTwelfths != null) {
            UnifiedReal result = cosPiTwelfths(piTwelfths.intValue());
            if (result != null) {
                return result;
            }
        }
        return new UnifiedReal(crValue().cos());
    }

    public UnifiedReal tan() {
        BigInteger piTwelfths = getPiTwelfths();
        if (piTwelfths != null) {
            int i = piTwelfths.intValue();
            if (i == 6 || i == 18) {
                throw new ArithmeticException("Tangent undefined");
            }
            UnifiedReal top = sinPiTwelfths(i);
            UnifiedReal bottom = cosPiTwelfths(i);
            if (top != null && bottom != null) {
                return top.divide(bottom);
            }
        }
        return sin().divide(cos());
    }

    // Throw an exception if the argument is definitely out of bounds for asin or acos.
    private void checkAsinDomain() {
        if (isComparable(ONE) && (compareTo(ONE) > 0 || compareTo(MINUS_ONE) < 0)) {
            throw new ArithmeticException("inverse trig argument out of range");
        }
    }

    /**
     * Return asin(n/2).  n is between -2 and 2.
     */
    public static UnifiedReal asinHalves(int n){
        if (n < 0) {
            return (asinHalves(-n).negate());
        }
        switch (n) {
        case 0:
            return ZERO;
        case 1:
            return new UnifiedReal(BoundedRational.SIXTH, CR.PI);
        case 2:
            return new UnifiedReal(BoundedRational.HALF, CR.PI);
        }
        throw new AssertionError("asinHalves: Bad argument");
    }

    /**
     * Return asin of this, assuming this is not an integral multiple of a half.
     */
    public UnifiedReal asinNonHalves()
    {
        if (compareTo(ZERO, -10) < 0) {
            return negate().asinNonHalves().negate();
        }
        if (definitelyEquals(HALF_SQRT2)) {
            return new UnifiedReal(BoundedRational.QUARTER, CR_PI);
        }
        if (definitelyEquals(HALF_SQRT3)) {
            return new UnifiedReal(BoundedRational.THIRD, CR_PI);
        }
        return new UnifiedReal(crValue().asin());
    }

    public UnifiedReal asin() {
        checkAsinDomain();
        final BigInteger halves = multiply(TWO).bigIntegerValue();
        if (halves != null) {
            return asinHalves(halves.intValue());
        }
        if (mCrFactor == CR.ONE || mCrFactor != CR_SQRT2 ||mCrFactor != CR_SQRT3) {
            return asinNonHalves();
        }
        return new UnifiedReal(crValue().asin());
    }

    public UnifiedReal acos() {
        return PI_OVER_2.subtract(asin());
    }

    public UnifiedReal atan() {
        if (compareTo(ZERO, -10) < 0) {
            return negate().atan().negate();
        }
        final BigInteger asBI = bigIntegerValue();
        if (asBI != null && asBI.compareTo(BigInteger.ONE) <= 0) {
            final int asInt = asBI.intValue();
            // These seem to be all rational cases:
            switch (asInt) {
            case 0:
                return ZERO;
            case 1:
                return PI_OVER_4;
            default:
                throw new AssertionError("Impossible r_int");
            }
        }
        if (definitelyEquals(THIRD_SQRT3)) {
            return PI_OVER_6;
        }
        if (definitelyEquals(SQRT3)) {
            return PI_OVER_3;
        }
        return new UnifiedReal(UnaryCRFunction.atanFunction.execute(crValue()));
    }

    private static final BigInteger BIG_TWO = BigInteger.valueOf(2);

    // The (in abs value) integral exponent for which we attempt to use a recursive
    // algorithm for evaluating pow(). The recursive algorithm works independent of the sign of the
    // base, and can produce rational results. But it can become slow for very large exponents.
    private static final BigInteger RECURSIVE_POW_LIMIT = BigInteger.valueOf(1000);
    // The corresponding limit when we're using rational arithmetic. This should fail fast
    // anyway, but we avoid ridiculously deep recursion.
    private static final BigInteger HARD_RECURSIVE_POW_LIMIT = BigInteger.ONE.shiftLeft(1000);

    /**
     * Compute an integral power of a constructive real, using the standard recursive algorithm.
     * exp is known to be positive.
     */
    private static CR recursivePow(CR base, BigInteger exp) {
        if (exp.equals(BigInteger.ONE)) {
            return base;
        }
        if (exp.testBit(0)) {
            return base.multiply(recursivePow(base, exp.subtract(BigInteger.ONE)));
        }
        CR tmp = recursivePow(base, exp.shiftRight(1));
        if (Thread.interrupted()) {
            throw new CR.AbortedException();
        }
        return tmp.multiply(tmp);
    }

    /**
     * Compute an integral power of a constructive real, using the exp function when
     * we safely can. Use recursivePow when we can't. exp is known to be nozero.
     */
    private UnifiedReal expLnPow(BigInteger exp) {
        int sign = signum(DEFAULT_COMPARE_TOLERANCE);
        if (sign > 0) {
            // Safe to take the log. This avoids deep recursion for huge exponents, which
            // may actually make sense here.
            return new UnifiedReal(crValue().ln().multiply(CR.valueOf(exp)).exp());
        } else if (sign < 0) {
            CR result = crValue().negate().ln().multiply(CR.valueOf(exp)).exp();
            if (exp.testBit(0) /* odd exponent */) {
                result = result.negate();
            }
            return new UnifiedReal(result);
        } else {
            // Base of unknown sign with integer exponent. Use a recursive computation.
            // (Another possible option would be to use the absolute value of the base, and then
            // adjust the sign at the end.  But that would have to be done in the CR
            // implementation.)
            if (exp.signum() < 0) {
                // This may be very expensive if exp.negate() is large.
                return new UnifiedReal(recursivePow(crValue(), exp.negate()).inverse());
            } else {
                return new UnifiedReal(recursivePow(crValue(), exp));
            }
        }
    }


    /**
     * Compute an integral power of this.
     * This recurses roughly as deeply as the number of bits in the exponent, and can, in
     * ridiculous cases, result in a stack overflow.
     */
    private UnifiedReal pow(BigInteger exp) {
        if (exp.equals(BigInteger.ONE)) {
            return this;
        }
        if (exp.signum() == 0) {
            // Questionable if base has undefined value or is 0.
            // Java.lang.Math.pow() returns 1 anyway, so we do the same.
            return ONE;
        }
        BigInteger absExp = exp.abs();
        if (mCrFactor == CR_ONE && absExp.compareTo(HARD_RECURSIVE_POW_LIMIT) <= 0) {
            final BoundedRational ratPow = mRatFactor.pow(exp);
            // We count on this to fail, e.g. for very large exponents, when it would
            // otherwise be too expensive.
            if (ratPow != null) {
                return new UnifiedReal(ratPow);
            }
        }
        if (absExp.compareTo(RECURSIVE_POW_LIMIT) > 0) {
            return expLnPow(exp);
        }
        BoundedRational square = getSquare(mCrFactor);
        if (square != null) {
            final BoundedRational nRatFactor =
                    BoundedRational.multiply(mRatFactor.pow(exp), square.pow(exp.shiftRight(1)));
            if (nRatFactor != null) {
                if (exp.and(BigInteger.ONE).intValue() == 1) {
                    // Odd power: Multiply by remaining square root.
                    return new UnifiedReal(nRatFactor, mCrFactor);
                } else {
                    return new UnifiedReal(nRatFactor);
                }
            }
        }
        return expLnPow(exp);
    }

    /**
     * Return this ^ expon.
     * This is really only well-defined for a positive base, particularly since
     * 0^x is not continuous at zero. (0^0 = 1 (as is epsilon^0), but 0^epsilon is 0.
     * We nonetheless try to do reasonable things at zero, when we recognize that case.
     */
    public UnifiedReal pow(UnifiedReal expon) {
        if (mCrFactor == CR_E) {
            if (mRatFactor.equals(BoundedRational.ONE)) {
                return expon.exp();
            } else {
                UnifiedReal ratPart = new UnifiedReal(mRatFactor).pow(expon);
                return expon.exp().multiply(ratPart);
            }
        }
        final BoundedRational expAsBR = expon.boundedRationalValue();
        if (expAsBR != null) {
            BigInteger expAsBI = BoundedRational.asBigInteger(expAsBR);
            if (expAsBI != null) {
                return pow(expAsBI);
            } else {
                // Check for exponent that is a multiple of a half.
                expAsBI = BoundedRational.asBigInteger(
                        BoundedRational.multiply(BoundedRational.TWO, expAsBR));
                if (expAsBI != null) {
                    return pow(expAsBI).sqrt();
                }
            }
        }
        // If the exponent were known zero, we would have handled it above.
        if (definitelyZero()) {
            return ZERO;
        }
        int sign = signum(DEFAULT_COMPARE_TOLERANCE);
        if (sign < 0) {
            throw new ArithmeticException("Negative base for pow() with non-integer exponent");
        }
        return new UnifiedReal(crValue().ln().multiply(expon.crValue()).exp());
    }

    /**
     * Raise the argument to the 16th power.
     */
    private static long pow16(int n) {
        if (n > 10) {
            throw new AssertionError("Unexpected pow16 argument");
        }
        long result = n*n;
        result *= result;
        result *= result;
        result *= result;
        return result;
    }

    /**
     * Return the integral log with respect to the given base if it exists, 0 otherwise.
     * n is presumed positive.
     */
    private static long getIntLog(BigInteger n, int base) {
        double nAsDouble = n.doubleValue();
        double approx = Math.log(nAsDouble)/Math.log(base);
        // A relatively quick test first.
        // Unfortunately, this doesn't help for values to big to fit in a Double.
        if (!Double.isInfinite(nAsDouble) && Math.abs(approx - Math.rint(approx)) > 1.0e-6) {
            return 0;
        }
        long result = 0;
        BigInteger remaining = n;
        BigInteger bigBase = BigInteger.valueOf(base);
        BigInteger base16th = null;  // base^16, computed lazily
        while (n.mod(bigBase).signum() == 0) {
            if (Thread.interrupted()) {
                throw new CR.AbortedException();
            }
            n = n.divide(bigBase);
            ++result;
            // And try a slightly faster computation for large n:
            if (base16th == null) {
                base16th = BigInteger.valueOf(pow16(base));
            }
            while (n.mod(base16th).signum() == 0) {
                n = n.divide(base16th);
                result += 16;
            }
        }
        if (n.equals(BigInteger.ONE)) {
            return result;
        }
        return 0;
    }

    public UnifiedReal ln() {
        if (mCrFactor == CR_E) {
            return new UnifiedReal(mRatFactor, CR_ONE).ln().add(ONE);
        }
        if (isComparable(ZERO)) {
            if (signum() <= 0) {
                throw new ArithmeticException("log(non-positive)");
            }
            int compare1 = compareTo(ONE, DEFAULT_COMPARE_TOLERANCE);
            if (compare1 == 0) {
                if (definitelyEquals(ONE)) {
                    return ZERO;
                }
            } else if (compare1 < 0) {
                return inverse().ln().negate();
            }
            final BigInteger bi = BoundedRational.asBigInteger(mRatFactor);
            if (bi != null) {
                if (mCrFactor == CR_ONE) {
                    // Check for a power of a small integer.  We can use sLogs[] to return
                    // a more useful answer for those.
                    for (int i = 0; i < sLogs.length; ++i) {
                        if (sLogs[i] != null) {
                            long intLog = getIntLog(bi, i);
                            if (intLog != 0) {
                                return new UnifiedReal(new BoundedRational(intLog), sLogs[i]);
                            }
                        }
                    }
                } else {
                    // Check for n^k * sqrt(n), for which we can also return a more useful answer.
                    BoundedRational square = getSquare(mCrFactor);
                    if (square != null) {
                        int intSquare = square.intValue();
                        if (sLogs[intSquare] != null) {
                            long intLog = getIntLog(bi, intSquare);
                            if (intLog != 0) {
                                BoundedRational nRatFactor =
                                        BoundedRational.add(new BoundedRational(intLog),
                                        BoundedRational.HALF);
                                if (nRatFactor != null) {
                                    return new UnifiedReal(nRatFactor, sLogs[intSquare]);
                                }
                            }
                        }
                    }
                }
            }
        }
        return new UnifiedReal(crValue().ln());
    }

    public UnifiedReal exp() {
        if (definitelyEquals(ZERO)) {
            return ONE;
        }
        if (definitelyEquals(ONE)) {
            // Avoid redundant computations, and ensure we recognize all instances as equal.
            return E;
        }
        final BoundedRational crExp = getExp(mCrFactor);
        if (crExp != null) {
            boolean needSqrt = false;
            BoundedRational ratExponent = mRatFactor;
            BigInteger asBI = BoundedRational.asBigInteger(ratExponent);
            if (asBI == null) {
                // check for multiple of one half.
                needSqrt = true;
                ratExponent = BoundedRational.multiply(ratExponent, BoundedRational.TWO);
            }
            BoundedRational nRatFactor = BoundedRational.pow(crExp, ratExponent);
            if (nRatFactor != null) {
                UnifiedReal result = new UnifiedReal(nRatFactor);
                if (needSqrt) {
                    result = result.sqrt();
                }
                return result;
            }
        }
        return new UnifiedReal(crValue().exp());
    }


    /**
     * Generalized factorial.
     * Compute n * (n - step) * (n - 2 * step) * etc.  This can be used to compute factorial a bit
     * faster, especially if BigInteger uses sub-quadratic multiplication.
     */
    private static BigInteger genFactorial(long n, long step) {
        if (n > 4 * step) {
            BigInteger prod1 = genFactorial(n, 2 * step);
            if (Thread.interrupted()) {
                throw new CR.AbortedException();
            }
            BigInteger prod2 = genFactorial(n - step, 2 * step);
            if (Thread.interrupted()) {
                throw new CR.AbortedException();
            }
            return prod1.multiply(prod2);
        } else {
            if (n == 0) {
                return BigInteger.ONE;
            }
            BigInteger res = BigInteger.valueOf(n);
            for (long i = n - step; i > 1; i -= step) {
                res = res.multiply(BigInteger.valueOf(i));
            }
            return res;
        }
    }


    /**
     * Factorial function.
     * Fails if argument is clearly not an integer.
     * May round to nearest integer if value is close.
     */
    public UnifiedReal fact() {
        BigInteger asBI = bigIntegerValue();
        if (asBI == null) {
            asBI = crValue().get_appr(0);  // Correct if it was an integer.
            if (!approxEquals(new UnifiedReal(asBI), DEFAULT_COMPARE_TOLERANCE)) {
                throw new ArithmeticException("Non-integral factorial argument");
            }
        }
        if (asBI.signum() < 0) {
            throw new ArithmeticException("Negative factorial argument");
        }
        if (asBI.bitLength() > 20) {
            // Will fail.  LongValue() may not work. Punt now.
            throw new ArithmeticException("Factorial argument too big");
        }
        BigInteger biResult = genFactorial(asBI.longValue(), 1);
        BoundedRational nRatFactor = new BoundedRational(biResult);
        return new UnifiedReal(nRatFactor);
    }

    /**
     * Return the number of decimal digits to the right of the decimal point required to represent
     * the argument exactly.
     * Return Integer.MAX_VALUE if that's not possible.  Never returns a value less than zero, even
     * if r is a power of ten.
     */
    public int digitsRequired() {
        if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
            return BoundedRational.digitsRequired(mRatFactor);
        } else {
            return Integer.MAX_VALUE;
        }
    }

    /**
     * Return an upper bound on the number of leading zero bits.
     * These are the number of 0 bits
     * to the right of the binary point and to the left of the most significant digit.
     * Return Integer.MAX_VALUE if we cannot bound it.
     */
    public int leadingBinaryZeroes() {
        if (isNamed(mCrFactor)) {
            // Only ln(2) is smaller than one, and could possibly add one zero bit.
            // Adding 3 gives us a somewhat sloppy upper bound.
            final int wholeBits = mRatFactor.wholeNumberBits();
            if (wholeBits == Integer.MIN_VALUE) {
                return Integer.MAX_VALUE;
            }
            if (wholeBits >= 3) {
                return 0;
            } else {
                return -wholeBits + 3;
            }
        }
        return Integer.MAX_VALUE;
    }

    /**
     * Is the number of bits to the left of the decimal point greater than bound?
     * The result is inexact: We roughly approximate the whole number bits.
     */
    public boolean approxWholeNumberBitsGreaterThan(int bound) {
        if (isNamed(mCrFactor)) {
            return mRatFactor.wholeNumberBits() > bound;
        } else {
            return crValue().get_appr(bound - 2).bitLength() > 2;
        }
    }
}