/* * 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 *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 * 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; } } }