blob: 351797c2ddc341f0a9dedf35ad7692786582e61b [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.
369 * The result includes n digits to the right of the decimal point.
370 * @param n result precision, >= 0
371 */
372 public String toStringTruncated(int n) {
373 if (mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO) {
374 return mRatFactor.toStringTruncated(n);
375 }
376 final CR scaled = CR.valueOf(BigInteger.TEN.pow(n)).multiply(crValue());
377 boolean negative = false;
378 BigInteger intScaled;
379 if (exactlyTruncatable()) {
380 intScaled = scaled.get_appr(0);
381 if (intScaled.signum() < 0) {
382 negative = true;
383 intScaled = intScaled.negate();
384 }
385 if (CR.valueOf(intScaled).compareTo(scaled.abs()) > 0) {
386 intScaled = intScaled.subtract(BigInteger.ONE);
387 }
388 check(CR.valueOf(intScaled).compareTo(scaled.abs()) < 0);
389 } else {
390 // Approximate case. Exact comparisons are impossible.
391 intScaled = scaled.get_appr(-EXTRA_PREC);
392 if (intScaled.signum() < 0) {
393 negative = true;
394 intScaled = intScaled.negate();
395 }
396 intScaled = intScaled.shiftRight(EXTRA_PREC);
397 }
398 String digits = intScaled.toString();
399 int len = digits.length();
400 if (len < n + 1) {
401 digits = BoundedRational.repeat('0', n + 1 - len) + digits;
402 len = n + 1;
403 }
404 return (negative ? "-" : "") + digits.substring(0, len - n) + "."
405 + digits.substring(len - n);
406 }
407
408 /*
409 * Can we compute correctly truncated approximations of this number?
410 */
411 public boolean exactlyTruncatable() {
412 // If the value is known rational, we can do exact comparisons.
413 // If the value is known irrational, then we can safely compare to rational approximations;
414 // equality is impossible; hence the comparison must converge.
415 // The only problem cases are the ones in which we don't know.
416 return mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO || definitelyIrrational();
417 }
418
419 /**
420 * Return a double approximation.
421 * TODO: Result is correctly rounded if known to be rational.
422 */
423 public double doubleValue() {
424 if (mCrFactor == CR_ONE) {
425 return mRatFactor.doubleValue(); // Hopefully correctly rounded
426 } else {
427 return crValue().doubleValue(); // Approximately correctly rounded
428 }
429 }
430
431 public CR crValue() {
432 return mRatFactor.crValue().multiply(mCrFactor);
433 }
434
435 /**
436 * Are this and r exactly comparable?
437 */
438 public boolean isComparable(UnifiedReal u) {
439 // We check for ONE only to speed up the common case.
440 // The use of a tolerance here means we can spuriously return false, not true.
441 return mCrFactor == u.mCrFactor
442 && (isNamed(mCrFactor) || mCrFactor.signum(DEFAULT_COMPARE_TOLERANCE) != 0)
443 || mRatFactor.signum() == 0 && u.mRatFactor.signum() == 0
444 || definitelyIndependent(mCrFactor, u.mCrFactor)
445 || crValue().compareTo(u.crValue(), DEFAULT_COMPARE_TOLERANCE) != 0;
446 }
447
448 /**
449 * Return +1 if this is greater than r, -1 if this is less than r, or 0 of the two are
450 * known to be equal.
451 * May diverge if the two are equal and !isComparable(r).
452 */
453 public int compareTo(UnifiedReal u) {
454 if (definitelyZero() && u.definitelyZero()) return 0;
455 if (mCrFactor == u.mCrFactor) {
456 int signum = mCrFactor.signum(); // Can diverge if mCRFactor == 0.
457 return signum * mRatFactor.compareTo(u.mRatFactor);
458 }
459 return crValue().compareTo(u.crValue()); // Can also diverge.
460 }
461
462 /**
463 * Return +1 if this is greater than r, -1 if this is less than r, or possibly 0 of the two are
464 * within 2^a of each other.
465 */
466 public int compareTo(UnifiedReal u, int a) {
467 if (isComparable(u)) {
468 return compareTo(u);
469 } else {
470 return crValue().compareTo(u.crValue(), a);
471 }
472 }
473
474 /**
475 * Return compareTo(ZERO, a).
476 */
477 public int signum(int a) {
478 return compareTo(ZERO, a);
479 }
480
481 /**
482 * Return compareTo(ZERO).
483 * May diverge for ZERO argument if !isComparable(ZERO).
484 */
485 public int signum() {
486 return compareTo(ZERO);
487 }
488
489 /**
490 * Equality comparison. May erroneously return true if values differ by less than 2^a,
491 * and !isComparable(u).
492 */
493 public boolean approxEquals(UnifiedReal u, int a) {
494 if (isComparable(u)) {
495 if (definitelyIndependent(mCrFactor, u.mCrFactor)
496 && (mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0)) {
497 // No need to actually evaluate, though we don't know which is larger.
498 return false;
499 } else {
500 return compareTo(u) == 0;
501 }
502 }
503 return crValue().compareTo(u.crValue(), a) == 0;
504 }
505
506 /**
507 * Returns true if values are definitely known to be equal, false in all other cases.
508 */
509 public boolean definitelyEquals(UnifiedReal u) {
510 return isComparable(u) && compareTo(u) == 0;
511 }
512
513 /**
514 * Returns true if values are definitely known not to be equal, false in all other cases.
515 */
516 public boolean definitelyNotEquals(UnifiedReal u) {
517 boolean isNamed = isNamed(mCrFactor);
518 boolean uIsNamed = isNamed(u.mCrFactor);
519 if (isNamed && uIsNamed) {
520 if (definitelyIndependent(mCrFactor, u.mCrFactor)) {
521 return mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0;
522 } else if (mCrFactor == u.mCrFactor) {
523 return !mRatFactor.equals(u.mRatFactor);
524 }
525 return !mRatFactor.equals(u.mRatFactor);
526 }
527 if (mRatFactor.signum() == 0) {
528 return uIsNamed && u.mRatFactor.signum() != 0;
529 }
530 if (u.mRatFactor.signum() == 0) {
531 return isNamed && mRatFactor.signum() != 0;
532 }
533 return false;
534 }
535
536 // And some slightly faster convenience functions for special cases:
537
538 public boolean definitelyZero() {
539 return mRatFactor.signum() == 0;
540 }
541
542 public boolean definitelyNonZero() {
543 return isNamed(mCrFactor) && mRatFactor.signum() != 0;
544 }
545
546 public boolean definitelyOne() {
547 return mCrFactor == CR_ONE && mRatFactor.equals(BoundedRational.ONE);
548 }
549
550 /**
551 * Return equivalent BoundedRational, if known to exist, null otherwise
552 */
553 public BoundedRational boundedRationalValue() {
554 if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
555 return mRatFactor;
556 }
557 return null;
558 }
559
560 /**
561 * Returns equivalent BigInteger result if it exists, null if not.
562 */
563 public BigInteger bigIntegerValue() {
564 final BoundedRational r = boundedRationalValue();
565 return BoundedRational.asBigInteger(r);
566 }
567
568 public UnifiedReal add(UnifiedReal u) {
569 if (mCrFactor == u.mCrFactor) {
570 BoundedRational nRatFactor = BoundedRational.add(mRatFactor, u.mRatFactor);
571 if (nRatFactor != null) {
572 return new UnifiedReal(nRatFactor, mCrFactor);
573 }
574 }
575 if (definitelyZero()) {
576 // Avoid creating new mCrFactor, even if they don't currently match.
577 return u;
578 }
579 if (u.definitelyZero()) {
580 return this;
581 }
582 return new UnifiedReal(crValue().add(u.crValue()));
583 }
584
585 public UnifiedReal negate() {
586 return new UnifiedReal(BoundedRational.negate(mRatFactor), mCrFactor);
587 }
588
589 public UnifiedReal subtract(UnifiedReal u) {
590 return add(u.negate());
591 }
592
593 public UnifiedReal multiply(UnifiedReal u) {
594 // Preserve a preexisting mCrFactor when we can.
595 if (mCrFactor == CR_ONE) {
596 BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
597 if (nRatFactor != null) {
598 return new UnifiedReal(nRatFactor, u.mCrFactor);
599 }
600 }
601 if (u.mCrFactor == CR_ONE) {
602 BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
603 if (nRatFactor != null) {
604 return new UnifiedReal(nRatFactor, mCrFactor);
605 }
606 }
607 if (definitelyZero() || u.definitelyZero()) {
608 return ZERO;
609 }
610 if (mCrFactor == u.mCrFactor) {
611 BoundedRational square = getSquare(mCrFactor);
612 if (square != null) {
613 BoundedRational nRatFactor = BoundedRational.multiply(
614 BoundedRational.multiply(square, mRatFactor), u.mRatFactor);
615 if (nRatFactor != null) {
616 return new UnifiedReal(nRatFactor);
617 }
618 }
619 }
620 // Probably a bit cheaper to multiply component-wise.
621 BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor);
622 if (nRatFactor != null) {
623 return new UnifiedReal(nRatFactor, mCrFactor.multiply(u.mCrFactor));
624 }
625 return new UnifiedReal(crValue().multiply(u.crValue()));
626 }
627
628 public static class ZeroDivisionException extends ArithmeticException {
629 public ZeroDivisionException() {
630 super("Division by zero");
631 }
632 }
633
634 /**
635 * Return the reciprocal.
636 */
637 public UnifiedReal inverse() {
638 if (definitelyZero()) {
639 throw new ZeroDivisionException();
640 }
641 BoundedRational square = getSquare(mCrFactor);
642 if (square != null) {
643 // 1/sqrt(n) = sqrt(n)/n
644 BoundedRational nRatFactor = BoundedRational.inverse(
645 BoundedRational.multiply(mRatFactor, square));
646 if (nRatFactor != null) {
647 return new UnifiedReal(nRatFactor, mCrFactor);
648 }
649 }
650 return new UnifiedReal(BoundedRational.inverse(mRatFactor), mCrFactor.inverse());
651 }
652
653 public UnifiedReal divide(UnifiedReal u) {
654 if (mCrFactor == u.mCrFactor) {
655 if (u.definitelyZero()) {
656 throw new ZeroDivisionException();
657 }
658 BoundedRational nRatFactor = BoundedRational.divide(mRatFactor, u.mRatFactor);
659 if (nRatFactor != null) {
660 return new UnifiedReal(nRatFactor, CR_ONE);
661 }
662 }
663 return multiply(u.inverse());
664 }
665
666 public UnifiedReal sqrt() {
667 if (mCrFactor == CR_ONE) {
668 BoundedRational ratSqrt;
669 // Check for all arguments of the form <perfect rational square> * small_int,
670 // where small_int has a known sqrt. This includes the small_int = 1 case.
671 for (int divisor = 1; divisor < sSqrts.length; ++divisor) {
672 if (sSqrts[divisor] != null) {
673 ratSqrt = BoundedRational.sqrt(
674 BoundedRational.divide(mRatFactor, new BoundedRational(divisor)));
675 if (ratSqrt != null) {
676 return new UnifiedReal(ratSqrt, sSqrts[divisor]);
677 }
678 }
679 }
680 }
681 return new UnifiedReal(crValue().sqrt());
682 }
683
684 /**
685 * Return (this mod 2pi)/(pi/6) as a BigInteger, or null if that isn't easily possible.
686 */
687 private BigInteger getPiTwelfths() {
688 if (definitelyZero()) return BigInteger.ZERO;
689 if (mCrFactor == CR_PI) {
690 BigInteger quotient = BoundedRational.asBigInteger(
691 BoundedRational.multiply(mRatFactor, BoundedRational.TWELVE));
692 if (quotient == null) {
693 return null;
694 }
695 return quotient.mod(BIG_24);
696 }
697 return null;
698 }
699
700 /**
701 * Computer the sin() for an integer multiple n of pi/12, if easily representable.
702 * @param n value between 0 and 23 inclusive.
703 */
704 private static UnifiedReal sinPiTwelfths(int n) {
705 if (n >= 12) {
706 UnifiedReal negResult = sinPiTwelfths(n - 12);
707 return negResult == null ? null : negResult.negate();
708 }
709 switch (n) {
710 case 0:
711 return ZERO;
712 case 2: // 30 degrees
713 return HALF;
714 case 3: // 45 degrees
715 return HALF_SQRT2;
716 case 4: // 60 degrees
717 return HALF_SQRT3;
718 case 6:
719 return ONE;
720 case 8:
721 return HALF_SQRT3;
722 case 9:
723 return HALF_SQRT2;
724 case 10:
725 return HALF;
726 default:
727 return null;
728 }
729 }
730
731 public UnifiedReal sin() {
732 BigInteger piTwelfths = getPiTwelfths();
733 if (piTwelfths != null) {
734 UnifiedReal result = sinPiTwelfths(piTwelfths.intValue());
735 if (result != null) {
736 return result;
737 }
738 }
739 return new UnifiedReal(crValue().sin());
740 }
741
742 private static UnifiedReal cosPiTwelfths(int n) {
743 int sinArg = n + 6;
744 if (sinArg >= 24) {
745 sinArg -= 24;
746 }
747 return sinPiTwelfths(sinArg);
748 }
749
750 public UnifiedReal cos() {
751 BigInteger piTwelfths = getPiTwelfths();
752 if (piTwelfths != null) {
753 UnifiedReal result = cosPiTwelfths(piTwelfths.intValue());
754 if (result != null) {
755 return result;
756 }
757 }
758 return new UnifiedReal(crValue().cos());
759 }
760
761 public UnifiedReal tan() {
762 BigInteger piTwelfths = getPiTwelfths();
763 if (piTwelfths != null) {
764 int i = piTwelfths.intValue();
765 if (i == 6 || i == 18) {
766 throw new ArithmeticException("Tangent undefined");
767 }
768 UnifiedReal top = sinPiTwelfths(i);
769 UnifiedReal bottom = cosPiTwelfths(i);
770 if (top != null && bottom != null) {
771 return top.divide(bottom);
772 }
773 }
774 return sin().divide(cos());
775 }
776
777 // Throw an exception if the argument is definitely out of bounds for asin or acos.
778 private void checkAsinDomain() {
779 if (isComparable(ONE) && (compareTo(ONE) > 0 || compareTo(MINUS_ONE) < 0)) {
780 throw new ArithmeticException("inverse trig argument out of range");
781 }
782 }
783
784 /**
785 * Return asin(n/2). n is between -2 and 2.
786 */
787 public static UnifiedReal asinHalves(int n){
788 if (n < 0) {
789 return (asinHalves(-n).negate());
790 }
791 switch (n) {
792 case 0:
793 return ZERO;
794 case 1:
795 return new UnifiedReal(BoundedRational.SIXTH, CR.PI);
796 case 2:
797 return new UnifiedReal(BoundedRational.HALF, CR.PI);
798 }
799 throw new AssertionError("asinHalves: Bad argument");
800 }
801
802 /**
803 * Return asin of this, assuming this is not an integral multiple of a half.
804 */
805 public UnifiedReal asinNonHalves()
806 {
807 if (compareTo(ZERO, -10) < 0) {
808 return negate().asinNonHalves().negate();
809 }
810 if (definitelyEquals(HALF_SQRT2)) {
811 return new UnifiedReal(BoundedRational.QUARTER, CR_PI);
812 }
813 if (definitelyEquals(HALF_SQRT3)) {
814 return new UnifiedReal(BoundedRational.THIRD, CR_PI);
815 }
816 return new UnifiedReal(crValue().asin());
817 }
818
819 public UnifiedReal asin() {
820 checkAsinDomain();
821 final BigInteger halves = multiply(TWO).bigIntegerValue();
822 if (halves != null) {
823 return asinHalves(halves.intValue());
824 }
825 if (mCrFactor == CR.ONE || mCrFactor != CR_SQRT2 ||mCrFactor != CR_SQRT3) {
826 return asinNonHalves();
827 }
828 return new UnifiedReal(crValue().asin());
829 }
830
831 public UnifiedReal acos() {
832 return PI_OVER_2.subtract(asin());
833 }
834
835 public UnifiedReal atan() {
836 if (compareTo(ZERO, -10) < 0) {
837 return negate().atan().negate();
838 }
839 final BigInteger asBI = bigIntegerValue();
840 if (asBI != null && asBI.compareTo(BigInteger.ONE) <= 0) {
841 final int asInt = asBI.intValue();
842 // These seem to be all rational cases:
843 switch (asInt) {
844 case 0:
845 return ZERO;
846 case 1:
847 return PI_OVER_4;
848 default:
849 throw new AssertionError("Impossible r_int");
850 }
851 }
852 if (definitelyEquals(THIRD_SQRT3)) {
853 return PI_OVER_6;
854 }
855 if (definitelyEquals(SQRT3)) {
856 return PI_OVER_3;
857 }
858 return new UnifiedReal(UnaryCRFunction.atanFunction.execute(crValue()));
859 }
860
861 private static final BigInteger BIG_TWO = BigInteger.valueOf(2);
862
863 /**
864 * Compute an integral power of this.
865 */
866 private UnifiedReal pow(BigInteger exp) {
867 if (exp.signum() < 0) {
868 return pow(exp.negate()).inverse();
869 }
870 if (exp.equals(BigInteger.ONE)) {
871 return this;
872 }
873 if (exp.signum() == 0) {
874 // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway,
875 // so we do the same.
876 return ONE;
877 }
878 if (mCrFactor == CR_ONE) {
879 final BoundedRational ratPow = mRatFactor.pow(exp);
880 if (ratPow != null) {
881 return new UnifiedReal(mRatFactor.pow(exp));
882 }
883 }
884 BoundedRational square = getSquare(mCrFactor);
885 if (square != null) {
886 final BoundedRational nRatFactor =
887 BoundedRational.multiply(mRatFactor.pow(exp), square.pow(exp.shiftRight(1)));
888 if (nRatFactor != null) {
889 if (exp.and(BigInteger.ONE).intValue() == 1) {
890 // Odd power: Multiply by remaining square root.
891 return new UnifiedReal(nRatFactor, mCrFactor);
892 } else {
893 return new UnifiedReal(nRatFactor);
894 }
895 }
896 }
897 return new UnifiedReal(crValue().ln().multiply(CR.valueOf(exp)).exp());
898 }
899
900 public UnifiedReal pow(UnifiedReal expon) {
901 if (mCrFactor == CR_E) {
902 if (mRatFactor.equals(BoundedRational.ONE)) {
903 return expon.exp();
904 } else {
905 UnifiedReal ratPart = new UnifiedReal(mRatFactor).pow(expon);
906 return expon.exp().multiply(ratPart);
907 }
908 }
909 final BoundedRational expAsBR = expon.boundedRationalValue();
910 if (expAsBR != null) {
911 BigInteger expAsBI = BoundedRational.asBigInteger(expAsBR);
912 if (expAsBI != null) {
913 return pow(expAsBI);
914 } else {
915 // Check for exponent that is a multiple of a half.
916 expAsBI = BoundedRational.asBigInteger(
917 BoundedRational.multiply(BoundedRational.TWO, expAsBR));
918 if (expAsBI != null) {
919 return pow(expAsBI).sqrt();
920 }
921 }
922 }
923 return new UnifiedReal(crValue().ln().multiply(expon.crValue()).exp());
924 }
925
926 /**
927 * Raise the argument to the 16th power.
928 */
929 private static long pow16(int n) {
930 if (n > 10) {
931 throw new AssertionError("Unexpexted pow16 argument");
932 }
933 long result = n*n;
934 result *= result;
935 result *= result;
936 result *= result;
937 return result;
938 }
939
940 /**
941 * Return the integral log with respect to the given base if it exists, 0 otherwise.
942 * n is presumed positive.
943 */
944 private static long getIntLog(BigInteger n, int base) {
945 double nAsDouble = n.doubleValue();
946 double approx = Math.log(nAsDouble)/Math.log(base);
947 // A relatively quick test first.
948 // Unfortunately, this doesn't help for values to big to fit in a Double.
949 if (!Double.isInfinite(nAsDouble) && Math.abs(approx - Math.rint(approx)) > 1.0e-6) {
950 return 0;
951 }
952 long result = 0;
953 BigInteger remaining = n;
954 BigInteger bigBase = BigInteger.valueOf(base);
955 BigInteger base16th = null; // base^16, computed lazily
956 while (n.mod(bigBase).signum() == 0) {
957 if (Thread.interrupted()) {
958 throw new CR.AbortedException();
959 }
960 n = n.divide(bigBase);
961 ++result;
962 // And try a slightly faster computation for large n:
963 if (base16th == null) {
964 base16th = BigInteger.valueOf(pow16(base));
965 }
966 while (n.mod(base16th).signum() == 0) {
967 n = n.divide(base16th);
968 result += 16;
969 }
970 }
971 if (n.equals(BigInteger.ONE)) {
972 return result;
973 }
974 return 0;
975 }
976
977 public UnifiedReal ln() {
978 if (isComparable(ZERO)) {
979 if (signum() <= 0) {
980 throw new ArithmeticException("log(non-positive)");
981 }
982 int compare1 = compareTo(ONE, DEFAULT_COMPARE_TOLERANCE);
983 if (compare1 == 0) {
984 if (definitelyEquals(ONE)) {
985 return ZERO;
986 }
987 } else if (compare1 < 0) {
988 return inverse().ln().negate();
989 }
990 final BigInteger bi = BoundedRational.asBigInteger(mRatFactor);
991 if (bi != null) {
992 if (mCrFactor == CR_ONE) {
993 // Check for a power of a small integer. We can use sLogs[] to return
994 // a more useful answer for those.
995 for (int i = 0; i < sLogs.length; ++i) {
996 if (sLogs[i] != null) {
997 long intLog = getIntLog(bi, i);
998 if (intLog != 0) {
999 return new UnifiedReal(new BoundedRational(intLog), sLogs[i]);
1000 }
1001 }
1002 }
1003 } else {
1004 // Check for n^k * sqrt(n), for which we can also return a more useful answer.
1005 BoundedRational square = getSquare(mCrFactor);
1006 if (square != null) {
1007 int intSquare = square.intValue();
1008 if (sLogs[intSquare] != null) {
1009 long intLog = getIntLog(bi, intSquare);
1010 if (intLog != 0) {
1011 BoundedRational nRatFactor =
1012 BoundedRational.add(new BoundedRational(intLog),
1013 BoundedRational.HALF);
1014 if (nRatFactor != null) {
1015 return new UnifiedReal(nRatFactor, sLogs[intSquare]);
1016 }
1017 }
1018 }
1019 }
1020 }
1021 }
1022 }
1023 return new UnifiedReal(crValue().ln());
1024 }
1025
1026 public UnifiedReal exp() {
1027 if (definitelyEquals(ZERO)) {
1028 return ONE;
1029 }
1030 final BoundedRational crExp = getExp(mCrFactor);
1031 if (crExp != null) {
1032 if (mRatFactor.signum() < 0) {
1033 return negate().exp().inverse();
1034 }
1035 boolean needSqrt = false;
1036 BoundedRational ratExponent = mRatFactor;
1037 BigInteger asBI = BoundedRational.asBigInteger(ratExponent);
1038 if (asBI == null) {
1039 // check for multiple of one half.
1040 needSqrt = true;
1041 ratExponent = BoundedRational.multiply(ratExponent, BoundedRational.TWO);
1042 }
1043 BoundedRational nRatFactor = BoundedRational.pow(crExp, ratExponent);
1044 if (nRatFactor != null) {
1045 UnifiedReal result = new UnifiedReal(nRatFactor);
1046 if (needSqrt) {
1047 result = result.sqrt();
1048 }
1049 return result;
1050 }
1051 }
1052 return new UnifiedReal(crValue().exp());
1053 }
1054
1055
1056 /**
1057 * Generalized factorial.
1058 * Compute n * (n - step) * (n - 2 * step) * etc. This can be used to compute factorial a bit
1059 * faster, especially if BigInteger uses sub-quadratic multiplication.
1060 */
1061 private static BigInteger genFactorial(long n, long step) {
1062 if (n > 4 * step) {
1063 BigInteger prod1 = genFactorial(n, 2 * step);
1064 if (Thread.interrupted()) {
1065 throw new CR.AbortedException();
1066 }
1067 BigInteger prod2 = genFactorial(n - step, 2 * step);
1068 if (Thread.interrupted()) {
1069 throw new CR.AbortedException();
1070 }
1071 return prod1.multiply(prod2);
1072 } else {
1073 if (n == 0) {
1074 return BigInteger.ONE;
1075 }
1076 BigInteger res = BigInteger.valueOf(n);
1077 for (long i = n - step; i > 1; i -= step) {
1078 res = res.multiply(BigInteger.valueOf(i));
1079 }
1080 return res;
1081 }
1082 }
1083
1084
1085 /**
1086 * Factorial function.
1087 * Fails if argument is clearly not an integer.
1088 * May round to nearest integer if value is close.
1089 */
1090 public UnifiedReal fact() {
1091 BigInteger asBI = bigIntegerValue();
1092 if (asBI == null) {
1093 asBI = crValue().get_appr(0); // Correct if it was an integer.
1094 if (!approxEquals(new UnifiedReal(asBI), DEFAULT_COMPARE_TOLERANCE)) {
1095 throw new ArithmeticException("Non-integral factorial argument");
1096 }
1097 }
1098 if (asBI.signum() < 0) {
1099 throw new ArithmeticException("Negative factorial argument");
1100 }
1101 if (asBI.bitLength() > 20) {
1102 // Will fail. LongValue() may not work. Punt now.
1103 throw new ArithmeticException("Factorial argument too big");
1104 }
1105 BigInteger biResult = genFactorial(asBI.longValue(), 1);
1106 BoundedRational nRatFactor = new BoundedRational(biResult);
1107 return new UnifiedReal(nRatFactor);
1108 }
1109
1110 /**
1111 * Return the number of decimal digits to the right of the decimal point required to represent
1112 * the argument exactly.
1113 * Return Integer.MAX_VALUE if that's not possible. Never returns a value less than zero, even
1114 * if r is a power of ten.
1115 */
1116 public int digitsRequired() {
1117 if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) {
1118 return BoundedRational.digitsRequired(mRatFactor);
1119 } else {
1120 return Integer.MAX_VALUE;
1121 }
1122 }
1123
1124 /**
1125 * Return an upper bound on the number of leading zero bits.
1126 * These are the number of 0 bits
1127 * to the right of the binary point and to the left of the most significant digit.
1128 * Return Integer.MAX_VALUE if we cannot bound it.
1129 */
1130 public int leadingBinaryZeroes() {
1131 if (isNamed(mCrFactor)) {
1132 // Only ln(2) is smaller than one, and could possibly add one zero bit.
1133 // Adding 3 gives us a somewhat sloppy upper bound.
1134 final int wholeBits = mRatFactor.wholeNumberBits();
1135 if (wholeBits == Integer.MIN_VALUE) {
1136 return Integer.MAX_VALUE;
1137 }
1138 if (wholeBits >= 3) {
1139 return 0;
1140 } else {
1141 return -wholeBits + 3;
1142 }
1143 }
1144 return Integer.MAX_VALUE;
1145 }
1146
1147 /**
1148 * Is the number of bits to the left of the decimal point greater than bound?
1149 * The result is inexact: We roughly approximate the whole number bits.
1150 */
1151 public boolean approxWholeNumberBitsGreaterThan(int bound) {
1152 if (isNamed(mCrFactor)) {
1153 return mRatFactor.wholeNumberBits() > bound;
1154 } else {
1155 return crValue().get_appr(bound - 2).bitLength() > 2;
1156 }
1157 }
1158}