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