• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
check(boolean b)56     private static void check(boolean b) {
57         if (!b) {
58             throw new AssertionError();
59         }
60     }
61 
UnifiedReal(BoundedRational rat, CR cr)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 
UnifiedReal(CR cr)71     public UnifiedReal(CR cr) {
72         this(BoundedRational.ONE, cr);
73     }
74 
UnifiedReal(BoundedRational rat)75     public UnifiedReal(BoundedRational rat) {
76         this(rat, CR_ONE);
77     }
78 
UnifiedReal(BigInteger n)79     public UnifiedReal(BigInteger n) {
80         this(new BoundedRational(n));
81     }
82 
UnifiedReal(long n)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      */
getSquare(CR cr)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      */
getExp(CR cr)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      */
crName(CR cr)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      */
isNamed(CR cr)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      */
definitelyAlgebraic(CR cr)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      */
definitelyRational()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      */
definitelyIrrational()262     public boolean definitelyIrrational() {
263         return !definitelyRational() && isNamed(mCrFactor);
264     }
265 
266     /**
267      * Is this number known to be algebraic?
268      */
definitelyAlgebraic()269     public boolean definitelyAlgebraic() {
270         return definitelyAlgebraic(mCrFactor) || mRatFactor.signum() == 0;
271     }
272 
273     /**
274      * Is this number known to be transcendental?
275      */
definitelyTranscendental()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      */
definitelyIndependent(CR r1, CR r2)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      */
toString()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      */
toNiceString()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      */
exactlyDisplayable()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      */
toStringTruncated(int n)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