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