blob: 7f4bfc85ad58bb1523c97d65a34763fdab9de04b [file] [log] [blame]
Hans Boehm995e5eb2016-02-08 11:03:01 -08001/*
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.
Hans Boehm3465c782016-12-12 17:28:10 -0800369 * Always includes a decimal point in the result.
Hans Boehm995e5eb2016-02-08 11:03:01 -0800370 * 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) {
Hans Boehm24c91ed2016-06-30 18:53:44 -0700402 digits = StringUtils.repeat('0', n + 1 - len) + digits;
Hans Boehm995e5eb2016-02-08 11:03:01 -0800403 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.
Hans Boehm4452c782016-12-07 14:52:05 -0800516 * Performs no approximate evaluation.
Hans Boehm995e5eb2016-02-08 11:03:01 -0800517 */
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
Hans Boehm4452c782016-12-07 14:52:05 -0800544 /**
545 * Can this number be determined to be definitely nonzero without performing approximate
546 * evaluation?
547 */
Hans Boehm995e5eb2016-02-08 11:03:01 -0800548 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 /**
Hans Boehm4452c782016-12-07 14:52:05 -0800870 * 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 /**
Hans Boehm995e5eb2016-02-08 11:03:01 -0800888 * Compute an integral power of this.
Hans Boehm4452c782016-12-07 14:52:05 -0800889 * This recurses roughly as deeply as the number of bits in the exponent, and can, in
890 * ridiculous cases, result in a stack overflow.
Hans Boehm995e5eb2016-02-08 11:03:01 -0800891 */
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 }
Hans Boehm4452c782016-12-07 14:52:05 -0800923 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 }
Hans Boehm995e5eb2016-02-08 11:03:01 -0800934 }
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 final BoundedRational crExp = getExp(mCrFactor);
1067 if (crExp != null) {
1068 if (mRatFactor.signum() < 0) {
1069 return negate().exp().inverse();
1070 }
1071 boolean needSqrt = false;
1072 BoundedRational ratExponent = mRatFactor;
1073 BigInteger asBI = BoundedRational.asBigInteger(ratExponent);
1074 if (asBI == null) {
1075 // check for multiple of one half.
1076 needSqrt = true;
1077 ratExponent = BoundedRational.multiply(ratExponent, BoundedRational.TWO);
1078 }
1079 BoundedRational nRatFactor = BoundedRational.pow(crExp, ratExponent);
1080 if (nRatFactor != null) {
1081 UnifiedReal result = new UnifiedReal(nRatFactor);
1082 if (needSqrt) {
1083 result = result.sqrt();
1084 }
1085 return result;
1086 }
1087 }
1088 return new UnifiedReal(crValue().exp());
1089 }
1090
1091
1092 /**
1093 * Generalized factorial.
1094 * Compute n * (n - step) * (n - 2 * step) * etc. This can be used to compute factorial a bit
1095 * faster, especially if BigInteger uses sub-quadratic multiplication.
1096 */
1097 private static BigInteger genFactorial(long n, long step) {
1098 if (n > 4 * step) {
1099 BigInteger prod1 = genFactorial(n, 2 * step);
1100 if (Thread.interrupted()) {
1101 throw new CR.AbortedException();
1102 }
1103 BigInteger prod2 = genFactorial(n - step, 2 * step);
1104 if (Thread.interrupted()) {
1105 throw new CR.AbortedException();
1106 }
1107 return prod1.multiply(prod2);
1108 } else {
1109 if (n == 0) {
1110 return BigInteger.ONE;
1111 }
1112 BigInteger res = BigInteger.valueOf(n);
1113 for (long i = n - step; i > 1; i -= step) {
1114 res = res.multiply(BigInteger.valueOf(i));
1115 }
1116 return res;
1117 }
1118 }
1119
1120
1121 /**
1122 * Factorial function.
1123 * Fails if argument is clearly not an integer.
1124 * May round to nearest integer if value is close.
1125 */
1126 public UnifiedReal fact() {
1127 BigInteger asBI = bigIntegerValue();
1128 if (asBI == null) {
1129 asBI = crValue().get_appr(0); // Correct if it was an integer.
1130 if (!approxEquals(new UnifiedReal(asBI), DEFAULT_COMPARE_TOLERANCE)) {
1131 throw new ArithmeticException("Non-integral factorial argument");
1132 }
1133 }
1134 if (asBI.signum() < 0) {
1135 throw new ArithmeticException("Negative factorial argument");
1136 }
1137 if (asBI.bitLength() > 20) {
1138 // Will fail. LongValue() may not work. Punt now.
1139 throw new ArithmeticException("Factorial argument too big");
1140 }
1141 BigInteger biResult = genFactorial(asBI.longValue(), 1);
1142 BoundedRational nRatFactor = new BoundedRational(biResult);
1143 return new UnifiedReal(nRatFactor);
1144 }
1145
1146 /**
1147 * Return the number of decimal digits to the right of the decimal point required to represent
1148 * the argument exactly.
1149 * Return Integer.MAX_VALUE if that's not possible. Never returns a value less than zero, even
1150 * if r is a power of ten.
1151 */
1152 public int digitsRequired() {
1153 if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
1154 return BoundedRational.digitsRequired(mRatFactor);
1155 } else {
1156 return Integer.MAX_VALUE;
1157 }
1158 }
1159
1160 /**
1161 * Return an upper bound on the number of leading zero bits.
1162 * These are the number of 0 bits
1163 * to the right of the binary point and to the left of the most significant digit.
1164 * Return Integer.MAX_VALUE if we cannot bound it.
1165 */
1166 public int leadingBinaryZeroes() {
1167 if (isNamed(mCrFactor)) {
1168 // Only ln(2) is smaller than one, and could possibly add one zero bit.
1169 // Adding 3 gives us a somewhat sloppy upper bound.
1170 final int wholeBits = mRatFactor.wholeNumberBits();
1171 if (wholeBits == Integer.MIN_VALUE) {
1172 return Integer.MAX_VALUE;
1173 }
1174 if (wholeBits >= 3) {
1175 return 0;
1176 } else {
1177 return -wholeBits + 3;
1178 }
1179 }
1180 return Integer.MAX_VALUE;
1181 }
1182
1183 /**
1184 * Is the number of bits to the left of the decimal point greater than bound?
1185 * The result is inexact: We roughly approximate the whole number bits.
1186 */
1187 public boolean approxWholeNumberBitsGreaterThan(int bound) {
1188 if (isNamed(mCrFactor)) {
1189 return mRatFactor.wholeNumberBits() > bound;
1190 } else {
1191 return crValue().get_appr(bound - 2).bitLength() > 2;
1192 }
1193 }
1194}