UnifiedReal.java revision fcb1e611dcc6ffd98be6b6847cbeb56966b70721
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
17package com.android.calculator2;
18
19import java.math.BigInteger;
20import com.hp.creals.CR;
21import 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 */
42public 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