Extend BoundedRational and UnifiedReal and fix pow(), exp()android-o-mr1-iot-preview-6o-mr1-iot-preview-6

Bug: 64852569

Clean up exp() and pow() implementation so that it always
terminates in reasonable time, if the result has reasonable size.
Do not negate the exponent and take the reciprocal of the
result, since that makes some near-zero results extremely
expensive to compute.

Lean more towards using the exp() and ln() implementation
of pow() when there is a choice. The recursive implementation
can be problematic with huge exponents.

Add accurate conversion function from Double.

Improve accuracy of conversion in the other direction. The old version
was prone to just saying NaN or Infinite for rationals with long
representations.

Do something more reasonable about hashCode() and equals() for
these two types, to make them safer to use.

Some minor cleanups, including some minor performance fixes.

A lot of this was driven by attempts to compare UnifiedReal results
to java.lang.Math generated double results.

Test: Ran unit tests (on host only, so far).
Change-Id: If2e47d99841b3b1fec2349acb31608136a71828e
(cherry picked from commit ff47d599ad952853119dbffbe1e4a8ac06924b4a)

2 files changed, 259 insertions, 43 deletions

diff --git a/src/com/android/calculator2/BoundedRational.java b/src/com/android/calculator2/BoundedRational.java
index f3452a2..16fa581 100644
--- a/src/com/android/calculator2/BoundedRational.java
+++ b/src/com/android/calculator2/BoundedRational.java
@@ -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;
 
 /**
@@ -62,6 +63,60 @@ public class BoundedRational {
     }
 
     /**
+     * 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.
      */
@@ -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() {
@@ -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;
@@ -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
@@ -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.
@@ -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;
         }
@@ -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
diff --git a/src/com/android/calculator2/UnifiedReal.java b/src/com/android/calculator2/UnifiedReal.java
index 61cac29..f85dd3e 100644
--- a/src/com/android/calculator2/UnifiedReal.java
+++ b/src/com/android/calculator2/UnifiedReal.java
@@ -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;
@@ -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.
      */
@@ -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) {
@@ -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.
@@ -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,
@@ -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));
@@ -885,28 +934,62 @@ public class UnifiedReal {
     }
 
     /**
+     * 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) {
             final BoundedRational nRatFactor =
@@ -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)) {
@@ -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());
     }
 
@@ -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;
@@ -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);