Loading src/com/android/calculator2/BoundedRational.java +148 −16 Original line number Diff line number Diff line Loading @@ -19,6 +19,7 @@ package com.android.calculator2; import com.hp.creals.CR; import java.math.BigInteger; import java.util.Objects; import java.util.Random; /** Loading Loading @@ -61,6 +62,60 @@ public class BoundedRational { mDen = BigInteger.valueOf(1); } /** * Produce BoundedRational equal to the given double. */ public static BoundedRational valueOf(double x) { final long l = Math.round(x); if ((double) l == x && Math.abs(l) <= 1000) { return valueOf(l); } final long allBits = Double.doubleToRawLongBits(Math.abs(x)); long mantissa = (allBits & ((1L << 52) - 1)); final int biased_exp = (int)(allBits >>> 52); if ((biased_exp & 0x7ff) == 0x7ff) { throw new ArithmeticException("Infinity or NaN not convertible to BoundedRational"); } final long sign = x < 0.0 ? -1 : 1; int exp = biased_exp - 1075; // 1023 + 52; we treat mantissa as integer. if (biased_exp == 0) { exp += 1; // Denormal exponent is 1 greater. } else { mantissa += (1L << 52); // Implied leading one. } BigInteger num = BigInteger.valueOf(sign * mantissa); BigInteger den = BigInteger.ONE; if (exp >= 0) { num = num.shiftLeft(exp); } else { den = den.shiftLeft(-exp); } return new BoundedRational(num, den); } /** * Produce BoundedRational equal to the given long. */ public static BoundedRational valueOf(long x) { if (x >= -2 && x <= 10) { switch((int) x) { case -2: return MINUS_TWO; case -1: return MINUS_ONE; case 0: return ZERO; case 1: return ONE; case 2: return TWO; case 10: return TEN; } } return new BoundedRational(x); } /** * Convert to String reflecting raw representation. * Debug or log messages only, not pretty. Loading Loading @@ -108,12 +163,50 @@ public class BoundedRational { /** * Return a double approximation. * The result is correctly rounded if numerator and denominator are * exactly representable as a double. * TODO: This should always be correctly rounded. * The result is correctly rounded to nearest, with ties rounded away from zero. * TODO: Should round ties to even. */ public double doubleValue() { return mNum.doubleValue() / mDen.doubleValue(); final int sign = signum(); if (sign < 0) { return -BoundedRational.negate(this).doubleValue(); } // We get the mantissa by dividing the numerator by denominator, after // suitably prescaling them so that the integral part of the result contains // enough bits. We do the prescaling to avoid any precision loss, so the division result // is correctly truncated towards zero. final int apprExp = mNum.bitLength() - mDen.bitLength(); if (apprExp < -1100 || sign == 0) { // Bail fast for clearly zero result. return 0.0; } final int neededPrec = apprExp - 80; final BigInteger dividend = neededPrec < 0 ? mNum.shiftLeft(-neededPrec) : mNum; final BigInteger divisor = neededPrec > 0 ? mDen.shiftLeft(neededPrec) : mDen; final BigInteger quotient = dividend.divide(divisor); final int qLength = quotient.bitLength(); int extraBits = qLength - 53; int exponent = neededPrec + qLength; // Exponent assuming leading binary point. if (exponent >= -1021) { // Binary point is actually to right of leading bit. --exponent; } else { // We're in the gradual underflow range. Drop more bits. extraBits += (-1022 - exponent) + 1; exponent = -1023; } final BigInteger bigMantissa = quotient.add(BigInteger.ONE.shiftLeft(extraBits - 1)).shiftRight(extraBits); if (exponent > 1024) { return Double.POSITIVE_INFINITY; } if (exponent > -1023 && bigMantissa.bitLength() != 53 || exponent <= -1023 && bigMantissa.bitLength() >= 53) { throw new AssertionError("doubleValue internal error"); } final long mantissa = bigMantissa.longValue(); final long bits = (mantissa & ((1l << 52) - 1)) | (((long) exponent + 1023) << 52); return Double.longBitsToDouble(bits); } public CR crValue() { Loading @@ -138,6 +231,10 @@ public class BoundedRational { } } /** * Is this number too big for us to continue with rational arithmetic? * We return fals for integers on the assumption that we have no better fallback. */ private boolean tooBig() { if (mDen.equals(BigInteger.ONE)) { return false; Loading Loading @@ -198,8 +295,16 @@ public class BoundedRational { return mNum.signum() * mDen.signum(); } public boolean equals(BoundedRational r) { return compareTo(r) == 0; @Override public int hashCode() { // Note that this may be too expensive to be useful. BoundedRational reduced = reduce().positiveDen(); return Objects.hash(reduced.mNum, reduced.mDen); } @Override public boolean equals(Object r) { return r != null && r instanceof BoundedRational && compareTo((BoundedRational) r) == 0; } // We use static methods for arithmetic, so that we can easily handle the null case. We try Loading Loading @@ -332,6 +437,7 @@ public class BoundedRational { public final static BoundedRational MINUS_NINETY = new BoundedRational(-90); private static final BigInteger BIG_TWO = BigInteger.valueOf(2); private static final BigInteger BIG_MINUS_ONE = BigInteger.valueOf(-1); /** * Compute integral power of this, assuming this has been reduced and exp is >= 0. Loading @@ -350,32 +456,59 @@ public class BoundedRational { if (Thread.interrupted()) { throw new CR.AbortedException(); } return rawMultiply(tmp, tmp); BoundedRational result = rawMultiply(tmp, tmp); if (result == null || result.tooBig()) { return null; } return result; } /** * Compute an integral power of this. */ public BoundedRational pow(BigInteger exp) { if (exp.signum() < 0) { return inverse(pow(exp.negate())); int expSign = exp.signum(); if (expSign == 0) { // Questionable if base has undefined or zero value. // java.lang.Math.pow() returns 1 anyway, so we do the same. return BoundedRational.ONE; } if (exp.equals(BigInteger.ONE)) { return this; } // Reducing once at the beginning means there's no point in reducing later. return reduce().rawPow(exp); BoundedRational reduced = reduce().positiveDen(); // First handle cases in which huge exponents could give compact results. if (reduced.mDen.equals(BigInteger.ONE)) { if (reduced.mNum.equals(BigInteger.ZERO)) { return ZERO; } if (reduced.mNum.equals(BigInteger.ONE)) { return ONE; } if (reduced.mNum.equals(BIG_MINUS_ONE)) { if (exp.testBit(0)) { return MINUS_ONE; } else { return ONE; } } } if (exp.bitLength() > 1000) { // Stack overflow is likely; a useful rational result is not. return null; } if (expSign < 0) { return inverse(reduced).rawPow(exp.negate()); } else { return reduced.rawPow(exp); } } public static BoundedRational pow(BoundedRational base, BoundedRational exp) { if (exp == null) { return null; } if (exp.mNum.signum() == 0) { // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway, // so we do the same. return new BoundedRational(1); } if (base == null) { return null; } Loading @@ -388,7 +521,6 @@ public class BoundedRational { private static final BigInteger BIG_FIVE = BigInteger.valueOf(5); private static final BigInteger BIG_MINUS_ONE = BigInteger.valueOf(-1); /** * Return the number of decimal digits to the right of the decimal point required to represent Loading src/com/android/calculator2/UnifiedReal.java +111 −27 Original line number Diff line number Diff line Loading @@ -84,6 +84,23 @@ public class UnifiedReal { 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; Loading Loading @@ -173,8 +190,8 @@ public class UnifiedReal { } /** * 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. * 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. */ Loading Loading @@ -419,7 +436,8 @@ public class UnifiedReal { /** * Return a double approximation. * TODO: Result is correctly rounded if known to be rational. * Rational arguments are currently rounded to nearest, with ties away from zero. * TODO: Improve rounding. */ public double doubleValue() { if (mCrFactor == CR_ONE) { Loading Loading @@ -506,11 +524,27 @@ public class UnifiedReal { /** * 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. Loading Loading @@ -669,7 +703,14 @@ public class UnifiedReal { 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, Loading Loading @@ -866,15 +907,23 @@ public class UnifiedReal { 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 constrive real, using the standard recursive algorithm. * 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.and(BigInteger.ONE).intValue() == 1) { if (exp.testBit(0)) { return base.multiply(recursivePow(base, exp.subtract(BigInteger.ONE))); } CR tmp = recursivePow(base, exp.shiftRight(1)); Loading @@ -884,28 +933,62 @@ public class UnifiedReal { 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.signum() < 0) { return pow(exp.negate()).inverse(); } if (exp.equals(BigInteger.ONE)) { return this; } if (exp.signum() == 0) { // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway, // so we do the same. // Questionable if base has undefined value or is 0. // Java.lang.Math.pow() returns 1 anyway, so we do the same. return ONE; } if (mCrFactor == CR_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(mRatFactor.pow(exp)); return new UnifiedReal(ratPow); } } if (absExp.compareTo(RECURSIVE_POW_LIMIT) > 0) { return expLnPow(exp); } BoundedRational square = getSquare(mCrFactor); if (square != null) { Loading @@ -920,19 +1003,15 @@ public class UnifiedReal { } } } if (signum(DEFAULT_COMPARE_TOLERANCE) > 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 { // Possibly negative base 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.) return new UnifiedReal(recursivePow(crValue(), exp)); } 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)) { Loading @@ -956,6 +1035,14 @@ public class UnifiedReal { } } } // 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()); } Loading @@ -964,7 +1051,7 @@ public class UnifiedReal { */ private static long pow16(int n) { if (n > 10) { throw new AssertionError("Unexpexted pow16 argument"); throw new AssertionError("Unexpected pow16 argument"); } long result = n*n; result *= result; Loading Loading @@ -1072,9 +1159,6 @@ public class UnifiedReal { } final BoundedRational crExp = getExp(mCrFactor); if (crExp != null) { if (mRatFactor.signum() < 0) { return negate().exp().inverse(); } boolean needSqrt = false; BoundedRational ratExponent = mRatFactor; BigInteger asBI = BoundedRational.asBigInteger(ratExponent); Loading Loading
src/com/android/calculator2/BoundedRational.java +148 −16 Original line number Diff line number Diff line Loading @@ -19,6 +19,7 @@ package com.android.calculator2; import com.hp.creals.CR; import java.math.BigInteger; import java.util.Objects; import java.util.Random; /** Loading Loading @@ -61,6 +62,60 @@ public class BoundedRational { mDen = BigInteger.valueOf(1); } /** * Produce BoundedRational equal to the given double. */ public static BoundedRational valueOf(double x) { final long l = Math.round(x); if ((double) l == x && Math.abs(l) <= 1000) { return valueOf(l); } final long allBits = Double.doubleToRawLongBits(Math.abs(x)); long mantissa = (allBits & ((1L << 52) - 1)); final int biased_exp = (int)(allBits >>> 52); if ((biased_exp & 0x7ff) == 0x7ff) { throw new ArithmeticException("Infinity or NaN not convertible to BoundedRational"); } final long sign = x < 0.0 ? -1 : 1; int exp = biased_exp - 1075; // 1023 + 52; we treat mantissa as integer. if (biased_exp == 0) { exp += 1; // Denormal exponent is 1 greater. } else { mantissa += (1L << 52); // Implied leading one. } BigInteger num = BigInteger.valueOf(sign * mantissa); BigInteger den = BigInteger.ONE; if (exp >= 0) { num = num.shiftLeft(exp); } else { den = den.shiftLeft(-exp); } return new BoundedRational(num, den); } /** * Produce BoundedRational equal to the given long. */ public static BoundedRational valueOf(long x) { if (x >= -2 && x <= 10) { switch((int) x) { case -2: return MINUS_TWO; case -1: return MINUS_ONE; case 0: return ZERO; case 1: return ONE; case 2: return TWO; case 10: return TEN; } } return new BoundedRational(x); } /** * Convert to String reflecting raw representation. * Debug or log messages only, not pretty. Loading Loading @@ -108,12 +163,50 @@ public class BoundedRational { /** * Return a double approximation. * The result is correctly rounded if numerator and denominator are * exactly representable as a double. * TODO: This should always be correctly rounded. * The result is correctly rounded to nearest, with ties rounded away from zero. * TODO: Should round ties to even. */ public double doubleValue() { return mNum.doubleValue() / mDen.doubleValue(); final int sign = signum(); if (sign < 0) { return -BoundedRational.negate(this).doubleValue(); } // We get the mantissa by dividing the numerator by denominator, after // suitably prescaling them so that the integral part of the result contains // enough bits. We do the prescaling to avoid any precision loss, so the division result // is correctly truncated towards zero. final int apprExp = mNum.bitLength() - mDen.bitLength(); if (apprExp < -1100 || sign == 0) { // Bail fast for clearly zero result. return 0.0; } final int neededPrec = apprExp - 80; final BigInteger dividend = neededPrec < 0 ? mNum.shiftLeft(-neededPrec) : mNum; final BigInteger divisor = neededPrec > 0 ? mDen.shiftLeft(neededPrec) : mDen; final BigInteger quotient = dividend.divide(divisor); final int qLength = quotient.bitLength(); int extraBits = qLength - 53; int exponent = neededPrec + qLength; // Exponent assuming leading binary point. if (exponent >= -1021) { // Binary point is actually to right of leading bit. --exponent; } else { // We're in the gradual underflow range. Drop more bits. extraBits += (-1022 - exponent) + 1; exponent = -1023; } final BigInteger bigMantissa = quotient.add(BigInteger.ONE.shiftLeft(extraBits - 1)).shiftRight(extraBits); if (exponent > 1024) { return Double.POSITIVE_INFINITY; } if (exponent > -1023 && bigMantissa.bitLength() != 53 || exponent <= -1023 && bigMantissa.bitLength() >= 53) { throw new AssertionError("doubleValue internal error"); } final long mantissa = bigMantissa.longValue(); final long bits = (mantissa & ((1l << 52) - 1)) | (((long) exponent + 1023) << 52); return Double.longBitsToDouble(bits); } public CR crValue() { Loading @@ -138,6 +231,10 @@ public class BoundedRational { } } /** * Is this number too big for us to continue with rational arithmetic? * We return fals for integers on the assumption that we have no better fallback. */ private boolean tooBig() { if (mDen.equals(BigInteger.ONE)) { return false; Loading Loading @@ -198,8 +295,16 @@ public class BoundedRational { return mNum.signum() * mDen.signum(); } public boolean equals(BoundedRational r) { return compareTo(r) == 0; @Override public int hashCode() { // Note that this may be too expensive to be useful. BoundedRational reduced = reduce().positiveDen(); return Objects.hash(reduced.mNum, reduced.mDen); } @Override public boolean equals(Object r) { return r != null && r instanceof BoundedRational && compareTo((BoundedRational) r) == 0; } // We use static methods for arithmetic, so that we can easily handle the null case. We try Loading Loading @@ -332,6 +437,7 @@ public class BoundedRational { public final static BoundedRational MINUS_NINETY = new BoundedRational(-90); private static final BigInteger BIG_TWO = BigInteger.valueOf(2); private static final BigInteger BIG_MINUS_ONE = BigInteger.valueOf(-1); /** * Compute integral power of this, assuming this has been reduced and exp is >= 0. Loading @@ -350,32 +456,59 @@ public class BoundedRational { if (Thread.interrupted()) { throw new CR.AbortedException(); } return rawMultiply(tmp, tmp); BoundedRational result = rawMultiply(tmp, tmp); if (result == null || result.tooBig()) { return null; } return result; } /** * Compute an integral power of this. */ public BoundedRational pow(BigInteger exp) { if (exp.signum() < 0) { return inverse(pow(exp.negate())); int expSign = exp.signum(); if (expSign == 0) { // Questionable if base has undefined or zero value. // java.lang.Math.pow() returns 1 anyway, so we do the same. return BoundedRational.ONE; } if (exp.equals(BigInteger.ONE)) { return this; } // Reducing once at the beginning means there's no point in reducing later. return reduce().rawPow(exp); BoundedRational reduced = reduce().positiveDen(); // First handle cases in which huge exponents could give compact results. if (reduced.mDen.equals(BigInteger.ONE)) { if (reduced.mNum.equals(BigInteger.ZERO)) { return ZERO; } if (reduced.mNum.equals(BigInteger.ONE)) { return ONE; } if (reduced.mNum.equals(BIG_MINUS_ONE)) { if (exp.testBit(0)) { return MINUS_ONE; } else { return ONE; } } } if (exp.bitLength() > 1000) { // Stack overflow is likely; a useful rational result is not. return null; } if (expSign < 0) { return inverse(reduced).rawPow(exp.negate()); } else { return reduced.rawPow(exp); } } public static BoundedRational pow(BoundedRational base, BoundedRational exp) { if (exp == null) { return null; } if (exp.mNum.signum() == 0) { // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway, // so we do the same. return new BoundedRational(1); } if (base == null) { return null; } Loading @@ -388,7 +521,6 @@ public class BoundedRational { private static final BigInteger BIG_FIVE = BigInteger.valueOf(5); private static final BigInteger BIG_MINUS_ONE = BigInteger.valueOf(-1); /** * Return the number of decimal digits to the right of the decimal point required to represent Loading
src/com/android/calculator2/UnifiedReal.java +111 −27 Original line number Diff line number Diff line Loading @@ -84,6 +84,23 @@ public class UnifiedReal { 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; Loading Loading @@ -173,8 +190,8 @@ public class UnifiedReal { } /** * 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. * 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. */ Loading Loading @@ -419,7 +436,8 @@ public class UnifiedReal { /** * Return a double approximation. * TODO: Result is correctly rounded if known to be rational. * Rational arguments are currently rounded to nearest, with ties away from zero. * TODO: Improve rounding. */ public double doubleValue() { if (mCrFactor == CR_ONE) { Loading Loading @@ -506,11 +524,27 @@ public class UnifiedReal { /** * 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. Loading Loading @@ -669,7 +703,14 @@ public class UnifiedReal { 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, Loading Loading @@ -866,15 +907,23 @@ public class UnifiedReal { 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 constrive real, using the standard recursive algorithm. * 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.and(BigInteger.ONE).intValue() == 1) { if (exp.testBit(0)) { return base.multiply(recursivePow(base, exp.subtract(BigInteger.ONE))); } CR tmp = recursivePow(base, exp.shiftRight(1)); Loading @@ -884,28 +933,62 @@ public class UnifiedReal { 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.signum() < 0) { return pow(exp.negate()).inverse(); } if (exp.equals(BigInteger.ONE)) { return this; } if (exp.signum() == 0) { // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway, // so we do the same. // Questionable if base has undefined value or is 0. // Java.lang.Math.pow() returns 1 anyway, so we do the same. return ONE; } if (mCrFactor == CR_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(mRatFactor.pow(exp)); return new UnifiedReal(ratPow); } } if (absExp.compareTo(RECURSIVE_POW_LIMIT) > 0) { return expLnPow(exp); } BoundedRational square = getSquare(mCrFactor); if (square != null) { Loading @@ -920,19 +1003,15 @@ public class UnifiedReal { } } } if (signum(DEFAULT_COMPARE_TOLERANCE) > 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 { // Possibly negative base 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.) return new UnifiedReal(recursivePow(crValue(), exp)); } 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)) { Loading @@ -956,6 +1035,14 @@ public class UnifiedReal { } } } // 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()); } Loading @@ -964,7 +1051,7 @@ public class UnifiedReal { */ private static long pow16(int n) { if (n > 10) { throw new AssertionError("Unexpexted pow16 argument"); throw new AssertionError("Unexpected pow16 argument"); } long result = n*n; result *= result; Loading Loading @@ -1072,9 +1159,6 @@ public class UnifiedReal { } final BoundedRational crExp = getExp(mCrFactor); if (crExp != null) { if (mRatFactor.signum() < 0) { return negate().exp().inverse(); } boolean needSqrt = false; BoundedRational ratExponent = mRatFactor; BigInteger asBI = BoundedRational.asBigInteger(ratExponent); Loading