Home | History | Annotate | Download | only in calculator2
      1 /*
      2  * Copyright (C) 2016 The Android Open Source Project
      3  *
      4  * Licensed under the Apache License, Version 2.0 (the "License");
      5  * you may not use this file except in compliance with the License.
      6  * You may obtain a copy of the License at
      7  *
      8  *   http://www.apache.org/licenses/LICENSE-2.0
      9  *
     10  * Unless required by applicable law or agreed to in writing, software
     11  * distributed under the License is distributed on an "AS IS" BASIS,
     12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
     13  * See the License for the specific language governing permissions and
     14  * limitations under the License.
     15  */
     16 
     17 package com.android.calculator2;
     18 
     19 import java.math.BigInteger;
     20 import com.hp.creals.CR;
     21 import com.hp.creals.UnaryCRFunction;
     22 
     23 /**
     24  * Computable real numbers, represented so that we can get exact decidable comparisons
     25  * for a number of interesting special cases, including rational computations.
     26  *
     27  * A real number is represented as the product of two numbers with different representations:
     28  * A) A BoundedRational that can only represent a subset of the rationals, but supports
     29  *    exact computable comparisons.
     30  * B) A lazily evaluated "constructive real number" that provides operations to evaluate
     31  *    itself to any requested number of digits.
     32  * Whenever possible, we choose (B) to be one of a small set of known constants about which we
     33  * know more.  For example, whenever we can, we represent rationals such that (B) is 1.
     34  * This scheme allows us to do some very limited symbolic computation on numbers when both
     35  * have the same (B) value, as well as in some other situations.  We try to maximize that
     36  * possibility.
     37  *
     38  * Arithmetic operations and operations that produce finite approximations may throw unchecked
     39  * exceptions produced by the underlying CR and BoundedRational packages, including
     40  * CR.PrecisionOverflowException and CR.AbortedException.
     41  */
     42 public class UnifiedReal {
     43 
     44     private final BoundedRational mRatFactor;
     45     private final CR mCrFactor;
     46     // TODO: It would be helpful to add flags to indicate whether the result is known
     47     // irrational, etc.  This sometimes happens even if mCrFactor is not one of the known ones.
     48     // And exact comparisons between rationals and known irrationals are decidable.
     49 
     50     /**
     51      * Perform some nontrivial consistency checks.
     52      * @hide
     53      */
     54     public static boolean enableChecks = true;
     55 
     56     private static void check(boolean b) {
     57         if (!b) {
     58             throw new AssertionError();
     59         }
     60     }
     61 
     62     private UnifiedReal(BoundedRational rat, CR cr) {
     63         if (rat == null) {
     64             throw new ArithmeticException("Building UnifiedReal from null");
     65         }
     66         // We don't normally traffic in null CRs, and hence don't test explicitly.
     67         mCrFactor = cr;
     68         mRatFactor = rat;
     69     }
     70 
     71     public UnifiedReal(CR cr) {
     72         this(BoundedRational.ONE, cr);
     73     }
     74 
     75     public UnifiedReal(BoundedRational rat) {
     76         this(rat, CR_ONE);
     77     }
     78 
     79     public UnifiedReal(BigInteger n) {
     80         this(new BoundedRational(n));
     81     }
     82 
     83     public UnifiedReal(long n) {
     84         this(new BoundedRational(n));
     85     }
     86 
     87     // Various helpful constants
     88     private final static BigInteger BIG_24 = BigInteger.valueOf(24);
     89     private final static int DEFAULT_COMPARE_TOLERANCE = -1000;
     90 
     91     // Well-known CR constants we try to use in the mCrFactor position:
     92     private final static CR CR_ONE = CR.ONE;
     93     private final static CR CR_PI = CR.PI;
     94     private final static CR CR_E = CR.ONE.exp();
     95     private final static CR CR_SQRT2 = CR.valueOf(2).sqrt();
     96     private final static CR CR_SQRT3 = CR.valueOf(3).sqrt();
     97     private final static CR CR_LN2 = CR.valueOf(2).ln();
     98     private final static CR CR_LN3 = CR.valueOf(3).ln();
     99     private final static CR CR_LN5 = CR.valueOf(5).ln();
    100     private final static CR CR_LN6 = CR.valueOf(6).ln();
    101     private final static CR CR_LN7 = CR.valueOf(7).ln();
    102     private final static CR CR_LN10 = CR.valueOf(10).ln();
    103 
    104     // Square roots that we try to recognize.
    105     // We currently recognize only a small fixed collection, since the sqrt() function needs to
    106     // identify numbers of the form <SQRT[i]>*n^2, and we don't otherwise know of a good
    107     // algorithm for that.
    108     private final static CR sSqrts[] = {
    109             null,
    110             CR.ONE,
    111             CR_SQRT2,
    112             CR_SQRT3,
    113             null,
    114             CR.valueOf(5).sqrt(),
    115             CR.valueOf(6).sqrt(),
    116             CR.valueOf(7).sqrt(),
    117             null,
    118             null,
    119             CR.valueOf(10).sqrt() };
    120 
    121     // Natural logs of small integers that we try to recognize.
    122     private final static CR sLogs[] = {
    123             null,
    124             null,
    125             CR_LN2,
    126             CR_LN3,
    127             null,
    128             CR_LN5,
    129             CR_LN6,
    130             CR_LN7,
    131             null,
    132             null,
    133             CR_LN10 };
    134 
    135 
    136     // Some convenient UnifiedReal constants.
    137     public static final UnifiedReal PI = new UnifiedReal(CR_PI);
    138     public static final UnifiedReal E = new UnifiedReal(CR_E);
    139     public static final UnifiedReal ZERO = new UnifiedReal(BoundedRational.ZERO);
    140     public static final UnifiedReal ONE = new UnifiedReal(BoundedRational.ONE);
    141     public static final UnifiedReal MINUS_ONE = new UnifiedReal(BoundedRational.MINUS_ONE);
    142     public static final UnifiedReal TWO = new UnifiedReal(BoundedRational.TWO);
    143     public static final UnifiedReal MINUS_TWO = new UnifiedReal(BoundedRational.MINUS_TWO);
    144     public static final UnifiedReal HALF = new UnifiedReal(BoundedRational.HALF);
    145     public static final UnifiedReal MINUS_HALF = new UnifiedReal(BoundedRational.MINUS_HALF);
    146     public static final UnifiedReal TEN = new UnifiedReal(BoundedRational.TEN);
    147     public static final UnifiedReal RADIANS_PER_DEGREE
    148             = new UnifiedReal(new BoundedRational(1, 180), CR_PI);
    149     private static final UnifiedReal SIX = new UnifiedReal(6);
    150     private static final UnifiedReal HALF_SQRT2 = new UnifiedReal(BoundedRational.HALF, CR_SQRT2);
    151     private static final UnifiedReal SQRT3 = new UnifiedReal(CR_SQRT3);
    152     private static final UnifiedReal HALF_SQRT3 = new UnifiedReal(BoundedRational.HALF, CR_SQRT3);
    153     private static final UnifiedReal THIRD_SQRT3 = new UnifiedReal(BoundedRational.THIRD, CR_SQRT3);
    154     private static final UnifiedReal PI_OVER_2 = new UnifiedReal(BoundedRational.HALF, CR_PI);
    155     private static final UnifiedReal PI_OVER_3 = new UnifiedReal(BoundedRational.THIRD, CR_PI);
    156     private static final UnifiedReal PI_OVER_4 = new UnifiedReal(BoundedRational.QUARTER, CR_PI);
    157     private static final UnifiedReal PI_OVER_6 = new UnifiedReal(BoundedRational.SIXTH, CR_PI);
    158 
    159 
    160     /**
    161      * Given a constructive real cr, try to determine whether cr is the square root of
    162      * a small integer.  If so, return its square as a BoundedRational.  Otherwise return null.
    163      * We make this determination by simple table lookup, so spurious null returns are
    164      * entirely possible, or even likely.
    165      */
    166     private static BoundedRational getSquare(CR cr) {
    167         for (int i = 0; i < sSqrts.length; ++i) {
    168              if (sSqrts[i] == cr) {
    169                 return new BoundedRational(i);
    170              }
    171         }
    172         return null;
    173     }
    174 
    175     /**
    176      * Given a constructive real cr, try to determine whether cr is the square root of
    177      * a small integer.  If so, return its square as a BoundedRational.  Otherwise return null.
    178      * We make this determination by simple table lookup, so spurious null returns are
    179      * entirely possible, or even likely.
    180      */
    181     private BoundedRational getExp(CR cr) {
    182         for (int i = 0; i < sLogs.length; ++i) {
    183              if (sLogs[i] == cr) {
    184                 return new BoundedRational(i);
    185              }
    186         }
    187         return null;
    188     }
    189 
    190     /**
    191      * If the argument is a well-known constructive real, return its name.
    192      * The name of "CR_ONE" is the empty string.
    193      * No named constructive reals are rational multiples of each other.
    194      * Thus two UnifiedReals with different named mCrFactors can be equal only if both
    195      * mRatFactors are zero or possibly if one is CR_PI and the other is CR_E.
    196      * (The latter is apparently an open problem.)
    197      */
    198     private static String crName(CR cr) {
    199         if (cr == CR_ONE) {
    200             return "";
    201         }
    202         if (cr == CR_PI) {
    203             return "\u03C0";   // GREEK SMALL LETTER PI
    204         }
    205         if (cr == CR_E) {
    206             return "e";
    207         }
    208         for (int i = 0; i < sSqrts.length; ++i) {
    209             if (cr == sSqrts[i]) {
    210                 return "\u221A" /* SQUARE ROOT */ + i;
    211             }
    212         }
    213         for (int i = 0; i < sLogs.length; ++i) {
    214             if (cr == sLogs[i]) {
    215                 return "ln(" + i + ")";
    216             }
    217         }
    218         return null;
    219     }
    220 
    221     /**
    222      * Would crName() return non-Null?
    223      */
    224     private static boolean isNamed(CR cr) {
    225         if (cr == CR_ONE || cr == CR_PI || cr == CR_E) {
    226             return true;
    227         }
    228         for (CR r: sSqrts) {
    229             if (cr == r) {
    230                 return true;
    231             }
    232         }
    233         for (CR r: sLogs) {
    234             if (cr == r) {
    235                 return true;
    236             }
    237         }
    238         return false;
    239     }
    240 
    241     /**
    242      * Is cr known to be algebraic (as opposed to transcendental)?
    243      * Currently only produces meaningful results for the above known special
    244      * constructive reals.
    245      */
    246     private static boolean definitelyAlgebraic(CR cr) {
    247         return cr == CR_ONE || getSquare(cr) != null;
    248     }
    249 
    250     /**
    251      * Is this number known to be rational?
    252      */
    253     public boolean definitelyRational() {
    254         return mCrFactor == CR_ONE || mRatFactor.signum() == 0;
    255     }
    256 
    257     /**
    258      * Is this number known to be irrational?
    259      * TODO: We could track the fact that something is irrational with an explicit flag, which
    260      * could cover many more cases.  Whether that matters in practice is TBD.
    261      */
    262     public boolean definitelyIrrational() {
    263         return !definitelyRational() && isNamed(mCrFactor);
    264     }
    265 
    266     /**
    267      * Is this number known to be algebraic?
    268      */
    269     public boolean definitelyAlgebraic() {
    270         return definitelyAlgebraic(mCrFactor) || mRatFactor.signum() == 0;
    271     }
    272 
    273     /**
    274      * Is this number known to be transcendental?
    275      */
    276     public boolean definitelyTranscendental() {
    277         return !definitelyAlgebraic() && isNamed(mCrFactor);
    278     }
    279 
    280 
    281     /**
    282      * Is it known that the two constructive reals differ by something other than a
    283      * a rational factor, i.e. is it known that two UnifiedReals
    284      * with those mCrFactors will compare unequal unless both mRatFactors are zero?
    285      * If this returns true, then a comparison of two UnifiedReals using those two
    286      * mCrFactors cannot diverge, though we don't know of a good runtime bound.
    287      */
    288     private static boolean definitelyIndependent(CR r1, CR r2) {
    289         // The question here is whether r1 = x*r2, where x is rational, where r1 and r2
    290         // are in our set of special known CRs, can have a solution.
    291         // This cannot happen if one is CR_ONE and the other is not.
    292         // (Since all others are irrational.)
    293         // This cannot happen for two named square roots, which have no repeated factors.
    294         // (To see this, square both sides of the equation and factor.  Each prime
    295         // factor in the numerator and denominator occurs twice.)
    296         // This cannot happen for e or pi on one side, and a square root on the other.
    297         // (One is transcendental, the other is algebraic.)
    298         // This cannot happen for two of our special natural logs.
    299         // (Otherwise ln(m) = (a/b)ln(n) ==> m = n^(a/b) ==> m^b = n^a, which is impossible
    300         // because either m or n includes a prime factor not shared by the other.)
    301         // This cannot happen for a log and a square root.
    302         // (The Lindemann-Weierstrass theorem tells us, among other things, that if
    303         // a is algebraic, then exp(a) is transcendental.  Thus if l in our finite
    304         // set of logs where algebraic, expl(l), must be transacendental.
    305         // But exp(l) is an integer.  Thus the logs are transcendental.  But of course the
    306         // square roots are algebraic.  Thus they can't be rational multiples.)
    307         // Unfortunately, we do not know whether e/pi is rational.
    308         if (r1 == r2) {
    309             return false;
    310         }
    311         CR other;
    312         if (r1 == CR_E || r1 == CR_PI) {
    313             return definitelyAlgebraic(r2);
    314         }
    315         if (r2 == CR_E || r2 == CR_PI) {
    316             return definitelyAlgebraic(r1);
    317         }
    318         return isNamed(r1) && isNamed(r2);
    319     }
    320 
    321     /**
    322      * Convert to String reflecting raw representation.
    323      * Debug or log messages only, not pretty.
    324      */
    325     public String toString() {
    326         return mRatFactor.toString() + "*" + mCrFactor.toString();
    327     }
    328 
    329     /**
    330      * Convert to readable String.
    331      * Intended for user output.  Produces exact expression when possible.
    332      */
    333     public String toNiceString() {
    334         if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
    335             return mRatFactor.toNiceString();
    336         }
    337         String name = crName(mCrFactor);
    338         if (name != null) {
    339             BigInteger bi = BoundedRational.asBigInteger(mRatFactor);
    340             if (bi != null) {
    341                 if (bi.equals(BigInteger.ONE)) {
    342                     return name;
    343                 }
    344                 return mRatFactor.toNiceString() + name;
    345             }
    346             return "(" + mRatFactor.toNiceString() + ")" + name;
    347         }
    348         if (mRatFactor.equals(BoundedRational.ONE)) {
    349             return mCrFactor.toString();
    350         }
    351         return crValue().toString();
    352     }
    353 
    354     /**
    355      * Will toNiceString() produce an exact representation?
    356      */
    357     public boolean exactlyDisplayable() {
    358         return crName(mCrFactor) != null;
    359     }
    360 
    361     // Number of extra bits used in evaluation below to prefer truncation to rounding.
    362     // Must be <= 30.
    363     private final static int EXTRA_PREC = 10;
    364 
    365     /*
    366      * Returns a truncated representation of the result.
    367      * If exactlyTruncatable(), we round correctly towards zero. Otherwise the resulting digit
    368      * string may occasionally be rounded up instead.
    369      * Always includes a decimal point in the result.
    370      * The result includes n digits to the right of the decimal point.
    371      * @param n result precision, >= 0
    372      */
    373     public String toStringTruncated(int n) {
    374         if (mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO) {
    375             return mRatFactor.toStringTruncated(n);
    376         }
    377         final CR scaled = CR.valueOf(BigInteger.TEN.pow(n)).multiply(crValue());
    378         boolean negative = false;
    379         BigInteger intScaled;
    380         if (exactlyTruncatable()) {
    381             intScaled = scaled.get_appr(0);
    382             if (intScaled.signum() < 0) {
    383                 negative = true;
    384                 intScaled = intScaled.negate();
    385             }
    386             if (CR.valueOf(intScaled).compareTo(scaled.abs()) > 0) {
    387                 intScaled = intScaled.subtract(BigInteger.ONE);
    388             }
    389             check(CR.valueOf(intScaled).compareTo(scaled.abs()) < 0);
    390         } else {
    391             // Approximate case.  Exact comparisons are impossible.
    392             intScaled = scaled.get_appr(-EXTRA_PREC);
    393             if (intScaled.signum() < 0) {
    394                 negative = true;
    395                 intScaled = intScaled.negate();
    396             }
    397             intScaled = intScaled.shiftRight(EXTRA_PREC);
    398         }
    399         String digits = intScaled.toString();
    400         int len = digits.length();
    401         if (len < n + 1) {
    402             digits = StringUtils.repeat('0', n + 1 - len) + digits;
    403             len = n + 1;
    404         }
    405         return (negative ? "-" : "") + digits.substring(0, len - n) + "."
    406                 + digits.substring(len - n);
    407     }
    408 
    409     /*
    410      * Can we compute correctly truncated approximations of this number?
    411      */
    412     public boolean exactlyTruncatable() {
    413         // If the value is known rational, we can do exact comparisons.
    414         // If the value is known irrational, then we can safely compare to rational approximations;
    415         // equality is impossible; hence the comparison must converge.
    416         // The only problem cases are the ones in which we don't know.
    417         return mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO || definitelyIrrational();
    418     }
    419 
    420     /**
    421      * Return a double approximation.
    422      * TODO: Result is correctly rounded if known to be rational.
    423      */
    424     public double doubleValue() {
    425         if (mCrFactor == CR_ONE) {
    426             return mRatFactor.doubleValue(); // Hopefully correctly rounded
    427         } else {
    428             return crValue().doubleValue(); // Approximately correctly rounded
    429         }
    430     }
    431 
    432     public CR crValue() {
    433         return mRatFactor.crValue().multiply(mCrFactor);
    434     }
    435 
    436     /**
    437      * Are this and r exactly comparable?
    438      */
    439     public boolean isComparable(UnifiedReal u) {
    440         // We check for ONE only to speed up the common case.
    441         // The use of a tolerance here means we can spuriously return false, not true.
    442         return mCrFactor == u.mCrFactor
    443                 && (isNamed(mCrFactor) || mCrFactor.signum(DEFAULT_COMPARE_TOLERANCE) != 0)
    444                 || mRatFactor.signum() == 0 && u.mRatFactor.signum() == 0
    445                 || definitelyIndependent(mCrFactor, u.mCrFactor)
    446                 || crValue().compareTo(u.crValue(), DEFAULT_COMPARE_TOLERANCE) != 0;
    447     }
    448 
    449     /**
    450      * Return +1 if this is greater than r, -1 if this is less than r, or 0 of the two are
    451      * known to be equal.
    452      * May diverge if the two are equal and !isComparable(r).
    453      */
    454     public int compareTo(UnifiedReal u) {
    455         if (definitelyZero() && u.definitelyZero()) return 0;
    456         if (mCrFactor == u.mCrFactor) {
    457             int signum = mCrFactor.signum();  // Can diverge if mCRFactor == 0.
    458             return signum * mRatFactor.compareTo(u.mRatFactor);
    459         }
    460         return crValue().compareTo(u.crValue());  // Can also diverge.
    461     }
    462 
    463     /**
    464      * Return +1 if this is greater than r, -1 if this is less than r, or possibly 0 of the two are
    465      * within 2^a of each other.
    466      */
    467     public int compareTo(UnifiedReal u, int a) {
    468         if (isComparable(u)) {
    469             return compareTo(u);
    470         } else {
    471             return crValue().compareTo(u.crValue(), a);
    472         }
    473     }
    474 
    475     /**
    476      * Return compareTo(ZERO, a).
    477      */
    478     public int signum(int a) {
    479         return compareTo(ZERO, a);
    480     }
    481 
    482     /**
    483      * Return compareTo(ZERO).
    484      * May diverge for ZERO argument if !isComparable(ZERO).
    485      */
    486     public int signum() {
    487         return compareTo(ZERO);
    488     }
    489 
    490     /**
    491      * Equality comparison.  May erroneously return true if values differ by less than 2^a,
    492      * and !isComparable(u).
    493      */
    494     public boolean approxEquals(UnifiedReal u, int a) {
    495         if (isComparable(u)) {
    496             if (definitelyIndependent(mCrFactor, u.mCrFactor)
    497                     && (mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0)) {
    498                 // No need to actually evaluate, though we don't know which is larger.
    499                 return false;
    500             } else {
    501                 return compareTo(u) == 0;
    502             }
    503         }
    504         return crValue().compareTo(u.crValue(), a) == 0;
    505     }
    506 
    507     /**
    508      * Returns true if values are definitely known to be equal, false in all other cases.
    509      */
    510     public boolean definitelyEquals(UnifiedReal u) {
    511         return isComparable(u) && compareTo(u) == 0;
    512     }
    513 
    514     /**
    515      * Returns true if values are definitely known not to be equal, false in all other cases.
    516      * Performs no approximate evaluation.
    517      */
    518     public boolean definitelyNotEquals(UnifiedReal u) {
    519         boolean isNamed = isNamed(mCrFactor);
    520         boolean uIsNamed = isNamed(u.mCrFactor);
    521         if (isNamed && uIsNamed) {
    522             if (definitelyIndependent(mCrFactor, u.mCrFactor)) {
    523                 return mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0;
    524             } else if (mCrFactor == u.mCrFactor) {
    525                 return !mRatFactor.equals(u.mRatFactor);
    526             }
    527             return !mRatFactor.equals(u.mRatFactor);
    528         }
    529         if (mRatFactor.signum() == 0) {
    530             return uIsNamed && u.mRatFactor.signum() != 0;
    531         }
    532         if (u.mRatFactor.signum() == 0) {
    533             return isNamed && mRatFactor.signum() != 0;
    534         }
    535         return false;
    536     }
    537 
    538     // And some slightly faster convenience functions for special cases:
    539 
    540     public boolean definitelyZero() {
    541         return mRatFactor.signum() == 0;
    542     }
    543 
    544     /**
    545      * Can this number be determined to be definitely nonzero without performing approximate
    546      * evaluation?
    547      */
    548     public boolean definitelyNonZero() {
    549         return isNamed(mCrFactor) && mRatFactor.signum() != 0;
    550     }
    551 
    552     public boolean definitelyOne() {
    553         return mCrFactor == CR_ONE && mRatFactor.equals(BoundedRational.ONE);
    554     }
    555 
    556     /**
    557      * Return equivalent BoundedRational, if known to exist, null otherwise
    558      */
    559     public BoundedRational boundedRationalValue() {
    560         if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
    561             return mRatFactor;
    562         }
    563         return null;
    564     }
    565 
    566     /**
    567      * Returns equivalent BigInteger result if it exists, null if not.
    568      */
    569     public BigInteger bigIntegerValue() {
    570         final BoundedRational r = boundedRationalValue();
    571         return BoundedRational.asBigInteger(r);
    572     }
    573 
    574     public UnifiedReal add(UnifiedReal u) {
    575         if (mCrFactor == u.mCrFactor) {
    576             BoundedRational nRatFactor = BoundedRational.add(mRatFactor, u.mRatFactor);
    577             if (nRatFactor != null) {
    578                 return new UnifiedReal(nRatFactor, mCrFactor);
    579             }
    580         }
    581         if (definitelyZero()) {
    582             // Avoid creating new mCrFactor, even if they don't currently match.
    583             return u;
    584         }
    585         if (u.definitelyZero()) {
    586             return this;
    587         }
    588         return new UnifiedReal(crValue().add(u.crValue()));
    589     }
    590 
    591     public UnifiedReal negate() {
    592         return new UnifiedReal(BoundedRational.negate(mRatFactor), mCrFactor);
    593     }
    594 
    595     public UnifiedReal subtract(UnifiedReal u) {
    596         return add(u.negate());
    597     }
    598 
    599     public UnifiedReal multiply(UnifiedReal u) {
    600         // Preserve a preexisting mCrFactor when we can.
    601         if (mCrFactor == CR_ONE) {
    602             BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
    603             if (nRatFactor != null) {
    604                 return new UnifiedReal(nRatFactor, u.mCrFactor);
    605             }
    606         }
    607         if (u.mCrFactor == CR_ONE) {
    608             BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
    609             if (nRatFactor != null) {
    610                 return new UnifiedReal(nRatFactor, mCrFactor);
    611             }
    612         }
    613         if (definitelyZero() || u.definitelyZero()) {
    614             return ZERO;
    615         }
    616         if (mCrFactor == u.mCrFactor) {
    617             BoundedRational square = getSquare(mCrFactor);
    618             if (square != null) {
    619                 BoundedRational nRatFactor = BoundedRational.multiply(
    620                         BoundedRational.multiply(square, mRatFactor), u.mRatFactor);
    621                 if (nRatFactor != null) {
    622                     return new UnifiedReal(nRatFactor);
    623                 }
    624             }
    625         }
    626         // Probably a bit cheaper to multiply component-wise.
    627         BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
    628         if (nRatFactor != null) {
    629             return new UnifiedReal(nRatFactor, mCrFactor.multiply(u.mCrFactor));
    630         }
    631         return new UnifiedReal(crValue().multiply(u.crValue()));
    632     }
    633 
    634     public static class ZeroDivisionException extends ArithmeticException {
    635         public ZeroDivisionException() {
    636             super("Division by zero");
    637         }
    638     }
    639 
    640     /**
    641      * Return the reciprocal.
    642      */
    643     public UnifiedReal inverse() {
    644         if (definitelyZero()) {
    645             throw new ZeroDivisionException();
    646         }
    647         BoundedRational square = getSquare(mCrFactor);
    648         if (square != null) {
    649             // 1/sqrt(n) = sqrt(n)/n
    650             BoundedRational nRatFactor = BoundedRational.inverse(
    651                     BoundedRational.multiply(mRatFactor, square));
    652             if (nRatFactor != null) {
    653                 return new UnifiedReal(nRatFactor, mCrFactor);
    654             }
    655         }
    656         return new UnifiedReal(BoundedRational.inverse(mRatFactor), mCrFactor.inverse());
    657     }
    658 
    659     public UnifiedReal divide(UnifiedReal u) {
    660         if (mCrFactor == u.mCrFactor) {
    661             if (u.definitelyZero()) {
    662                 throw new ZeroDivisionException();
    663             }
    664             BoundedRational nRatFactor = BoundedRational.divide(mRatFactor, u.mRatFactor);
    665             if (nRatFactor != null) {
    666                 return new UnifiedReal(nRatFactor, CR_ONE);
    667             }
    668         }
    669         return multiply(u.inverse());
    670     }
    671 
    672     public UnifiedReal sqrt() {
    673         if (mCrFactor == CR_ONE) {
    674             BoundedRational ratSqrt;
    675             // Check for all arguments of the form <perfect rational square> * small_int,
    676             // where small_int has a known sqrt.  This includes the small_int = 1 case.
    677             for (int divisor = 1; divisor < sSqrts.length; ++divisor) {
    678                 if (sSqrts[divisor] != null) {
    679                     ratSqrt = BoundedRational.sqrt(
    680                             BoundedRational.divide(mRatFactor, new BoundedRational(divisor)));
    681                     if (ratSqrt != null) {
    682                         return new UnifiedReal(ratSqrt, sSqrts[divisor]);
    683                     }
    684                 }
    685             }
    686         }
    687         return new UnifiedReal(crValue().sqrt());
    688     }
    689 
    690     /**
    691      * Return (this mod 2pi)/(pi/6) as a BigInteger, or null if that isn't easily possible.
    692      */
    693     private BigInteger getPiTwelfths() {
    694         if (definitelyZero()) return BigInteger.ZERO;
    695         if (mCrFactor == CR_PI) {
    696             BigInteger quotient = BoundedRational.asBigInteger(
    697                     BoundedRational.multiply(mRatFactor, BoundedRational.TWELVE));
    698             if (quotient == null) {
    699                 return null;
    700             }
    701             return quotient.mod(BIG_24);
    702         }
    703         return null;
    704     }
    705 
    706     /**
    707      * Computer the sin() for an integer multiple n of pi/12, if easily representable.
    708      * @param n value between 0 and 23 inclusive.
    709      */
    710     private static UnifiedReal sinPiTwelfths(int n) {
    711         if (n >= 12) {
    712             UnifiedReal negResult = sinPiTwelfths(n - 12);
    713             return negResult == null ? null : negResult.negate();
    714         }
    715         switch (n) {
    716         case 0:
    717             return ZERO;
    718         case 2: // 30 degrees
    719             return HALF;
    720         case 3: // 45 degrees
    721             return HALF_SQRT2;
    722         case 4: // 60 degrees
    723             return HALF_SQRT3;
    724         case 6:
    725             return ONE;
    726         case 8:
    727             return HALF_SQRT3;
    728         case 9:
    729             return HALF_SQRT2;
    730         case 10:
    731             return HALF;
    732         default:
    733             return null;
    734         }
    735     }
    736 
    737     public UnifiedReal sin() {
    738         BigInteger piTwelfths = getPiTwelfths();
    739         if (piTwelfths != null) {
    740             UnifiedReal result = sinPiTwelfths(piTwelfths.intValue());
    741             if (result != null) {
    742                 return result;
    743             }
    744         }
    745         return new UnifiedReal(crValue().sin());
    746     }
    747 
    748     private static UnifiedReal cosPiTwelfths(int n) {
    749         int sinArg = n + 6;
    750         if (sinArg >= 24) {
    751             sinArg -= 24;
    752         }
    753         return sinPiTwelfths(sinArg);
    754     }
    755 
    756     public UnifiedReal cos() {
    757         BigInteger piTwelfths = getPiTwelfths();
    758         if (piTwelfths != null) {
    759             UnifiedReal result = cosPiTwelfths(piTwelfths.intValue());
    760             if (result != null) {
    761                 return result;
    762             }
    763         }
    764         return new UnifiedReal(crValue().cos());
    765     }
    766 
    767     public UnifiedReal tan() {
    768         BigInteger piTwelfths = getPiTwelfths();
    769         if (piTwelfths != null) {
    770             int i = piTwelfths.intValue();
    771             if (i == 6 || i == 18) {
    772                 throw new ArithmeticException("Tangent undefined");
    773             }
    774             UnifiedReal top = sinPiTwelfths(i);
    775             UnifiedReal bottom = cosPiTwelfths(i);
    776             if (top != null && bottom != null) {
    777                 return top.divide(bottom);
    778             }
    779         }
    780         return sin().divide(cos());
    781     }
    782 
    783     // Throw an exception if the argument is definitely out of bounds for asin or acos.
    784     private void checkAsinDomain() {
    785         if (isComparable(ONE) && (compareTo(ONE) > 0 || compareTo(MINUS_ONE) < 0)) {
    786             throw new ArithmeticException("inverse trig argument out of range");
    787         }
    788     }
    789 
    790     /**
    791      * Return asin(n/2).  n is between -2 and 2.
    792      */
    793     public static UnifiedReal asinHalves(int n){
    794         if (n < 0) {
    795             return (asinHalves(-n).negate());
    796         }
    797         switch (n) {
    798         case 0:
    799             return ZERO;
    800         case 1:
    801             return new UnifiedReal(BoundedRational.SIXTH, CR.PI);
    802         case 2:
    803             return new UnifiedReal(BoundedRational.HALF, CR.PI);
    804         }
    805         throw new AssertionError("asinHalves: Bad argument");
    806     }
    807 
    808     /**
    809      * Return asin of this, assuming this is not an integral multiple of a half.
    810      */
    811     public UnifiedReal asinNonHalves()
    812     {
    813         if (compareTo(ZERO, -10) < 0) {
    814             return negate().asinNonHalves().negate();
    815         }
    816         if (definitelyEquals(HALF_SQRT2)) {
    817             return new UnifiedReal(BoundedRational.QUARTER, CR_PI);
    818         }
    819         if (definitelyEquals(HALF_SQRT3)) {
    820             return new UnifiedReal(BoundedRational.THIRD, CR_PI);
    821         }
    822         return new UnifiedReal(crValue().asin());
    823     }
    824 
    825     public UnifiedReal asin() {
    826         checkAsinDomain();
    827         final BigInteger halves = multiply(TWO).bigIntegerValue();
    828         if (halves != null) {
    829             return asinHalves(halves.intValue());
    830         }
    831         if (mCrFactor == CR.ONE || mCrFactor != CR_SQRT2 ||mCrFactor != CR_SQRT3) {
    832             return asinNonHalves();
    833         }
    834         return new UnifiedReal(crValue().asin());
    835     }
    836 
    837     public UnifiedReal acos() {
    838         return PI_OVER_2.subtract(asin());
    839     }
    840 
    841     public UnifiedReal atan() {
    842         if (compareTo(ZERO, -10) < 0) {
    843             return negate().atan().negate();
    844         }
    845         final BigInteger asBI = bigIntegerValue();
    846         if (asBI != null && asBI.compareTo(BigInteger.ONE) <= 0) {
    847             final int asInt = asBI.intValue();
    848             // These seem to be all rational cases:
    849             switch (asInt) {
    850             case 0:
    851                 return ZERO;
    852             case 1:
    853                 return PI_OVER_4;
    854             default:
    855                 throw new AssertionError("Impossible r_int");
    856             }
    857         }
    858         if (definitelyEquals(THIRD_SQRT3)) {
    859             return PI_OVER_6;
    860         }
    861         if (definitelyEquals(SQRT3)) {
    862             return PI_OVER_3;
    863         }
    864         return new UnifiedReal(UnaryCRFunction.atanFunction.execute(crValue()));
    865     }
    866 
    867     private static final BigInteger BIG_TWO = BigInteger.valueOf(2);
    868 
    869     /**
    870      * Compute an integral power of a constrive real, using the standard recursive algorithm.
    871      * exp is known to be positive.
    872      */
    873     private static CR recursivePow(CR base, BigInteger exp) {
    874         if (exp.equals(BigInteger.ONE)) {
    875             return base;
    876         }
    877         if (exp.and(BigInteger.ONE).intValue() == 1) {
    878             return base.multiply(recursivePow(base, exp.subtract(BigInteger.ONE)));
    879         }
    880         CR tmp = recursivePow(base, exp.shiftRight(1));
    881         if (Thread.interrupted()) {
    882             throw new CR.AbortedException();
    883         }
    884         return tmp.multiply(tmp);
    885     }
    886 
    887     /**
    888      * Compute an integral power of this.
    889      * This recurses roughly as deeply as the number of bits in the exponent, and can, in
    890      * ridiculous cases, result in a stack overflow.
    891      */
    892     private UnifiedReal pow(BigInteger exp) {
    893         if (exp.signum() < 0) {
    894             return pow(exp.negate()).inverse();
    895         }
    896         if (exp.equals(BigInteger.ONE)) {
    897             return this;
    898         }
    899         if (exp.signum() == 0) {
    900             // Questionable if base has undefined value.  Java.lang.Math.pow() returns 1 anyway,
    901             // so we do the same.
    902             return ONE;
    903         }
    904         if (mCrFactor == CR_ONE) {
    905             final BoundedRational ratPow = mRatFactor.pow(exp);
    906             if (ratPow != null) {
    907                 return new UnifiedReal(mRatFactor.pow(exp));
    908             }
    909         }
    910         BoundedRational square = getSquare(mCrFactor);
    911         if (square != null) {
    912             final BoundedRational nRatFactor =
    913                     BoundedRational.multiply(mRatFactor.pow(exp), square.pow(exp.shiftRight(1)));
    914             if (nRatFactor != null) {
    915                 if (exp.and(BigInteger.ONE).intValue() == 1) {
    916                     // Odd power: Multiply by remaining square root.
    917                     return new UnifiedReal(nRatFactor, mCrFactor);
    918                 } else {
    919                     return new UnifiedReal(nRatFactor);
    920                 }
    921             }
    922         }
    923         if (signum(DEFAULT_COMPARE_TOLERANCE) > 0) {
    924             // Safe to take the log. This avoids deep recursion for huge exponents, which
    925             // may actually make sense here.
    926             return new UnifiedReal(crValue().ln().multiply(CR.valueOf(exp)).exp());
    927         } else {
    928             // Possibly negative base with integer exponent. Use a recursive computation.
    929             // (Another possible option would be to use the absolute value of the base, and then
    930             // adjust the sign at the end.  But that would have to be done in the CR
    931             // implementation.)
    932             return new UnifiedReal(recursivePow(crValue(), exp));
    933         }
    934     }
    935 
    936     public UnifiedReal pow(UnifiedReal expon) {
    937         if (mCrFactor == CR_E) {
    938             if (mRatFactor.equals(BoundedRational.ONE)) {
    939                 return expon.exp();
    940             } else {
    941                 UnifiedReal ratPart = new UnifiedReal(mRatFactor).pow(expon);
    942                 return expon.exp().multiply(ratPart);
    943             }
    944         }
    945         final BoundedRational expAsBR = expon.boundedRationalValue();
    946         if (expAsBR != null) {
    947             BigInteger expAsBI = BoundedRational.asBigInteger(expAsBR);
    948             if (expAsBI != null) {
    949                 return pow(expAsBI);
    950             } else {
    951                 // Check for exponent that is a multiple of a half.
    952                 expAsBI = BoundedRational.asBigInteger(
    953                         BoundedRational.multiply(BoundedRational.TWO, expAsBR));
    954                 if (expAsBI != null) {
    955                     return pow(expAsBI).sqrt();
    956                 }
    957             }
    958         }
    959         return new UnifiedReal(crValue().ln().multiply(expon.crValue()).exp());
    960     }
    961 
    962     /**
    963      * Raise the argument to the 16th power.
    964      */
    965     private static long pow16(int n) {
    966         if (n > 10) {
    967             throw new AssertionError("Unexpexted pow16 argument");
    968         }
    969         long result = n*n;
    970         result *= result;
    971         result *= result;
    972         result *= result;
    973         return result;
    974     }
    975 
    976     /**
    977      * Return the integral log with respect to the given base if it exists, 0 otherwise.
    978      * n is presumed positive.
    979      */
    980     private static long getIntLog(BigInteger n, int base) {
    981         double nAsDouble = n.doubleValue();
    982         double approx = Math.log(nAsDouble)/Math.log(base);
    983         // A relatively quick test first.
    984         // Unfortunately, this doesn't help for values to big to fit in a Double.
    985         if (!Double.isInfinite(nAsDouble) && Math.abs(approx - Math.rint(approx)) > 1.0e-6) {
    986             return 0;
    987         }
    988         long result = 0;
    989         BigInteger remaining = n;
    990         BigInteger bigBase = BigInteger.valueOf(base);
    991         BigInteger base16th = null;  // base^16, computed lazily
    992         while (n.mod(bigBase).signum() == 0) {
    993             if (Thread.interrupted()) {
    994                 throw new CR.AbortedException();
    995             }
    996             n = n.divide(bigBase);
    997             ++result;
    998             // And try a slightly faster computation for large n:
    999             if (base16th == null) {
   1000                 base16th = BigInteger.valueOf(pow16(base));
   1001             }
   1002             while (n.mod(base16th).signum() == 0) {
   1003                 n = n.divide(base16th);
   1004                 result += 16;
   1005             }
   1006         }
   1007         if (n.equals(BigInteger.ONE)) {
   1008             return result;
   1009         }
   1010         return 0;
   1011     }
   1012 
   1013     public UnifiedReal ln() {
   1014         if (isComparable(ZERO)) {
   1015             if (signum() <= 0) {
   1016                 throw new ArithmeticException("log(non-positive)");
   1017             }
   1018             int compare1 = compareTo(ONE, DEFAULT_COMPARE_TOLERANCE);
   1019             if (compare1 == 0) {
   1020                 if (definitelyEquals(ONE)) {
   1021                     return ZERO;
   1022                 }
   1023             } else if (compare1 < 0) {
   1024                 return inverse().ln().negate();
   1025             }
   1026             final BigInteger bi = BoundedRational.asBigInteger(mRatFactor);
   1027             if (bi != null) {
   1028                 if (mCrFactor == CR_ONE) {
   1029                     // Check for a power of a small integer.  We can use sLogs[] to return
   1030                     // a more useful answer for those.
   1031                     for (int i = 0; i < sLogs.length; ++i) {
   1032                         if (sLogs[i] != null) {
   1033                             long intLog = getIntLog(bi, i);
   1034                             if (intLog != 0) {
   1035                                 return new UnifiedReal(new BoundedRational(intLog), sLogs[i]);
   1036                             }
   1037                         }
   1038                     }
   1039                 } else {
   1040                     // Check for n^k * sqrt(n), for which we can also return a more useful answer.
   1041                     BoundedRational square = getSquare(mCrFactor);
   1042                     if (square != null) {
   1043                         int intSquare = square.intValue();
   1044                         if (sLogs[intSquare] != null) {
   1045                             long intLog = getIntLog(bi, intSquare);
   1046                             if (intLog != 0) {
   1047                                 BoundedRational nRatFactor =
   1048                                         BoundedRational.add(new BoundedRational(intLog),
   1049                                         BoundedRational.HALF);
   1050                                 if (nRatFactor != null) {
   1051                                     return new UnifiedReal(nRatFactor, sLogs[intSquare]);
   1052                                 }
   1053                             }
   1054                         }
   1055                     }
   1056                 }
   1057             }
   1058         }
   1059         return new UnifiedReal(crValue().ln());
   1060     }
   1061 
   1062     public UnifiedReal exp() {
   1063         if (definitelyEquals(ZERO)) {
   1064             return ONE;
   1065         }
   1066         if (definitelyEquals(ONE)) {
   1067             // Avoid redundant computations, and ensure we recognize all instances as equal.
   1068             return E;
   1069         }
   1070         final BoundedRational crExp = getExp(mCrFactor);
   1071         if (crExp != null) {
   1072             if (mRatFactor.signum() < 0) {
   1073                 return negate().exp().inverse();
   1074             }
   1075             boolean needSqrt = false;
   1076             BoundedRational ratExponent = mRatFactor;
   1077             BigInteger asBI = BoundedRational.asBigInteger(ratExponent);
   1078             if (asBI == null) {
   1079                 // check for multiple of one half.
   1080                 needSqrt = true;
   1081                 ratExponent = BoundedRational.multiply(ratExponent, BoundedRational.TWO);
   1082             }
   1083             BoundedRational nRatFactor = BoundedRational.pow(crExp, ratExponent);
   1084             if (nRatFactor != null) {
   1085                 UnifiedReal result = new UnifiedReal(nRatFactor);
   1086                 if (needSqrt) {
   1087                     result = result.sqrt();
   1088                 }
   1089                 return result;
   1090             }
   1091         }
   1092         return new UnifiedReal(crValue().exp());
   1093     }
   1094 
   1095 
   1096     /**
   1097      * Generalized factorial.
   1098      * Compute n * (n - step) * (n - 2 * step) * etc.  This can be used to compute factorial a bit
   1099      * faster, especially if BigInteger uses sub-quadratic multiplication.
   1100      */
   1101     private static BigInteger genFactorial(long n, long step) {
   1102         if (n > 4 * step) {
   1103             BigInteger prod1 = genFactorial(n, 2 * step);
   1104             if (Thread.interrupted()) {
   1105                 throw new CR.AbortedException();
   1106             }
   1107             BigInteger prod2 = genFactorial(n - step, 2 * step);
   1108             if (Thread.interrupted()) {
   1109                 throw new CR.AbortedException();
   1110             }
   1111             return prod1.multiply(prod2);
   1112         } else {
   1113             if (n == 0) {
   1114                 return BigInteger.ONE;
   1115             }
   1116             BigInteger res = BigInteger.valueOf(n);
   1117             for (long i = n - step; i > 1; i -= step) {
   1118                 res = res.multiply(BigInteger.valueOf(i));
   1119             }
   1120             return res;
   1121         }
   1122     }
   1123 
   1124 
   1125     /**
   1126      * Factorial function.
   1127      * Fails if argument is clearly not an integer.
   1128      * May round to nearest integer if value is close.
   1129      */
   1130     public UnifiedReal fact() {
   1131         BigInteger asBI = bigIntegerValue();
   1132         if (asBI == null) {
   1133             asBI = crValue().get_appr(0);  // Correct if it was an integer.
   1134             if (!approxEquals(new UnifiedReal(asBI), DEFAULT_COMPARE_TOLERANCE)) {
   1135                 throw new ArithmeticException("Non-integral factorial argument");
   1136             }
   1137         }
   1138         if (asBI.signum() < 0) {
   1139             throw new ArithmeticException("Negative factorial argument");
   1140         }
   1141         if (asBI.bitLength() > 20) {
   1142             // Will fail.  LongValue() may not work. Punt now.
   1143             throw new ArithmeticException("Factorial argument too big");
   1144         }
   1145         BigInteger biResult = genFactorial(asBI.longValue(), 1);
   1146         BoundedRational nRatFactor = new BoundedRational(biResult);
   1147         return new UnifiedReal(nRatFactor);
   1148     }
   1149 
   1150     /**
   1151      * Return the number of decimal digits to the right of the decimal point required to represent
   1152      * the argument exactly.
   1153      * Return Integer.MAX_VALUE if that's not possible.  Never returns a value less than zero, even
   1154      * if r is a power of ten.
   1155      */
   1156     public int digitsRequired() {
   1157         if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
   1158             return BoundedRational.digitsRequired(mRatFactor);
   1159         } else {
   1160             return Integer.MAX_VALUE;
   1161         }
   1162     }
   1163 
   1164     /**
   1165      * Return an upper bound on the number of leading zero bits.
   1166      * These are the number of 0 bits
   1167      * to the right of the binary point and to the left of the most significant digit.
   1168      * Return Integer.MAX_VALUE if we cannot bound it.
   1169      */
   1170     public int leadingBinaryZeroes() {
   1171         if (isNamed(mCrFactor)) {
   1172             // Only ln(2) is smaller than one, and could possibly add one zero bit.
   1173             // Adding 3 gives us a somewhat sloppy upper bound.
   1174             final int wholeBits = mRatFactor.wholeNumberBits();
   1175             if (wholeBits == Integer.MIN_VALUE) {
   1176                 return Integer.MAX_VALUE;
   1177             }
   1178             if (wholeBits >= 3) {
   1179                 return 0;
   1180             } else {
   1181                 return -wholeBits + 3;
   1182             }
   1183         }
   1184         return Integer.MAX_VALUE;
   1185     }
   1186 
   1187     /**
   1188      * Is the number of bits to the left of the decimal point greater than bound?
   1189      * The result is inexact: We roughly approximate the whole number bits.
   1190      */
   1191     public boolean approxWholeNumberBitsGreaterThan(int bound) {
   1192         if (isNamed(mCrFactor)) {
   1193             return mRatFactor.wholeNumberBits() > bound;
   1194         } else {
   1195             return crValue().get_appr(bound - 2).bitLength() > 2;
   1196         }
   1197     }
   1198 }
   1199