1 /* 2 * Copyright (C) 2015 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 /* 18 * The above license covers additions and changes by AOSP authors. 19 * The original code is licensed as follows: 20 */ 21 22 // 23 // Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED 24 // 25 // Permission is granted free of charge to copy, modify, use and distribute 26 // this software provided you include the entirety of this notice in all 27 // copies made. 28 // 29 // THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY 30 // KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, 31 // WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT 32 // FOR A PARTICULAR PURPOSE OR NON-INFRINGING. SGI ASSUMES NO RISK AS TO THE 33 // QUALITY AND PERFORMANCE OF THE SOFTWARE. SHOULD THE SOFTWARE PROVE 34 // DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY 35 // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES 36 // AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS 37 // AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. 38 // 39 // UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, 40 // WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR 41 // OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, 42 // INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE 43 // SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK 44 // STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL 45 // OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF 46 // THE POSSIBILITY OF SUCH DAMAGES. THIS LIMITATION OF LIABILITY SHALL NOT 47 // APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE 48 // LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE 49 // EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT 50 // EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. 51 // 52 // These license terms shall be governed by and construed in accordance with 53 // the laws of the United States and the State of California as applied to 54 // agreements entered into and to be performed entirely within California 55 // between California residents. Any litigation relating to these license 56 // terms shall be subject to the exclusive jurisdiction of the Federal Courts 57 // of the Northern District of California (or, absent subject matter 58 // jurisdiction in such courts, the courts of the State of California), with 59 // venue lying exclusively in Santa Clara County, California. 60 61 // Copyright (c) 2001-2004, Hewlett-Packard Development Company, L.P. 62 // 63 // Permission is granted free of charge to copy, modify, use and distribute 64 // this software provided you include the entirety of this notice in all 65 // copies made. 66 // 67 // THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY 68 // KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, 69 // WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT 70 // FOR A PARTICULAR PURPOSE OR NON-INFRINGING. HEWLETT-PACKARD ASSUMES 71 // NO RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE. 72 // SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, 73 // HEWLETT-PACKARD ASSUMES NO COST OR LIABILITY FOR ANY 74 // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES 75 // AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS 76 // AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. 77 // 78 // UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, 79 // WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR 80 // OTHERWISE, SHALL HEWLETT-PACKARD BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, 81 // INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE 82 // SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK 83 // STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL 84 // OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF HEWLETT-PACKARD SHALL 85 // HAVE BEEN INFORMED OF THE POSSIBILITY OF SUCH DAMAGES. 86 // THIS LIMITATION OF LIABILITY SHALL NOT APPLY TO LIABILITY RESULTING 87 // FROM HEWLETT-PACKARD's NEGLIGENCE TO THE EXTENT APPLICABLE 88 // LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE 89 // EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT 90 // EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. 91 // 92 93 // Added valueOf(string, radix), fixed some documentation comments. 94 // Hans_Boehm@hp.com 1/12/2001 95 // Fixed a serious typo in inv_CR(): For negative arguments it produced 96 // the wrong sign. This affected the sign of divisions. 97 // Added byteValue and fixed some comments. Hans.Boehm@hp.com 12/17/2002 98 // Added toStringFloatRep. Hans.Boehm@hp.com 4/1/2004 99 // Added get_appr() synchronization to allow access from multiple threads 100 // hboehm@google.com 4/25/2014 101 // Changed cos() prescaling to avoid logarithmic depth tree. 102 // hboehm@google.com 6/30/2014 103 // Added explicit asin() implementation. Remove one. Add ZERO and ONE and 104 // make them public. hboehm@google.com 5/21/2015 105 // Added Gauss-Legendre PI implementation. Removed two. 106 // hboehm@google.com 4/12/2016 107 // Fix shift operation in doubleValue. That produced incorrect values for 108 // large negative exponents. 109 // Don't negate argument and compute inverse for exp(). That causes severe 110 // performance problems for (-huge).exp() 111 // hboehm@google.com 8/21/2017 112 // Have comparison check for interruption. hboehm@google.com 10/31/2017 113 114 package com.hp.creals; 115 116 import java.math.BigInteger; 117 import java.util.ArrayList; 118 119 /** 120 * Constructive real numbers, also known as recursive, or computable reals. 121 * Each recursive real number is represented as an object that provides an 122 * approximation function for the real number. 123 * The approximation function guarantees that the generated approximation 124 * is accurate to the specified precision. 125 * Arithmetic operations on constructive reals produce new such objects; 126 * they typically do not perform any real computation. 127 * In this sense, arithmetic computations are exact: They produce 128 * a description which describes the exact answer, and can be used to 129 * later approximate it to arbitrary precision. 130 * <P> 131 * When approximations are generated, <I>e.g.</i> for output, they are 132 * accurate to the requested precision; no cumulative rounding errors 133 * are visible. 134 * In order to achieve this precision, the approximation function will often 135 * need to approximate subexpressions to greater precision than was originally 136 * demanded. Thus the approximation of a constructive real number 137 * generated through a complex sequence of operations may eventually require 138 * evaluation to very high precision. This usually makes such computations 139 * prohibitively expensive for large numerical problems. 140 * But it is perfectly appropriate for use in a desk calculator, 141 * for small numerical problems, for the evaluation of expressions 142 * computated by a symbolic algebra system, for testing of accuracy claims 143 * for floating point code on small inputs, or the like. 144 * <P> 145 * We expect that the vast majority of uses will ignore the particular 146 * implementation, and the member functons <TT>approximate</tt> 147 * and <TT>get_appr</tt>. Such applications will treat <TT>CR</tt> as 148 * a conventional numerical type, with an interface modelled on 149 * <TT>java.math.BigInteger</tt>. No subclasses of <TT>CR</tt> 150 * will be explicitly mentioned by such a program. 151 * <P> 152 * All standard arithmetic operations, as well as a few algebraic 153 * and transcendal functions are provided. Constructive reals are 154 * immutable; thus all of these operations return a new constructive real. 155 * <P> 156 * A few uses will require explicit construction of approximation functions. 157 * The requires the construction of a subclass of <TT>CR</tt> with 158 * an overridden <TT>approximate</tt> function. Note that <TT>approximate</tt> 159 * should only be defined, but never called. <TT>get_appr</tt> 160 * provides the same functionality, but adds the caching necessary to obtain 161 * reasonable performance. 162 * <P> 163 * Any operation may throw <TT>com.hp.creals.AbortedException</tt> if the thread 164 * in which it is executing is interrupted. (<TT>InterruptedException</tt> 165 * cannot be used for this purpose, since CR inherits from <TT>Number</tt>.) 166 * <P> 167 * Any operation may also throw <TT>com.hp.creals.PrecisionOverflowException</tt> 168 * If the precision request generated during any subcalculation overflows 169 * a 28-bit integer. (This should be extremely unlikely, except as an 170 * outcome of a division by zero, or other erroneous computation.) 171 * 172 */ 173 public abstract class CR extends Number { 174 // CR is the basic representation of a number. 175 // Abstractly this is a function for computing an approximation 176 // plus the current best approximation. 177 // We could do without the latter, but that would 178 // be atrociously slow. 179 180 /** 181 * Indicates a constructive real operation was interrupted. 182 * Most constructive real operations may throw such an exception. 183 * This is unchecked, since Number methods may not raise checked 184 * exceptions. 185 */ 186 public static class AbortedException extends RuntimeException { AbortedException()187 public AbortedException() { super(); } AbortedException(String s)188 public AbortedException(String s) { super(s); } 189 } 190 191 /** 192 * Indicates that the number of bits of precision requested by 193 * a computation on constructive reals required more than 28 bits, 194 * and was thus in danger of overflowing an int. 195 * This is likely to be a symptom of a diverging computation, 196 * <I>e.g.</i> division by zero. 197 */ 198 public static class PrecisionOverflowException extends RuntimeException { PrecisionOverflowException()199 public PrecisionOverflowException() { super(); } PrecisionOverflowException(String s)200 public PrecisionOverflowException(String s) { super(s); } 201 } 202 203 // First some frequently used constants, so we don't have to 204 // recompute these all over the place. 205 static final BigInteger big0 = BigInteger.ZERO; 206 static final BigInteger big1 = BigInteger.ONE; 207 static final BigInteger bigm1 = BigInteger.valueOf(-1); 208 static final BigInteger big2 = BigInteger.valueOf(2); 209 static final BigInteger bigm2 = BigInteger.valueOf(-2); 210 static final BigInteger big3 = BigInteger.valueOf(3); 211 static final BigInteger big6 = BigInteger.valueOf(6); 212 static final BigInteger big8 = BigInteger.valueOf(8); 213 static final BigInteger big10 = BigInteger.TEN; 214 static final BigInteger big750 = BigInteger.valueOf(750); 215 static final BigInteger bigm750 = BigInteger.valueOf(-750); 216 217 /** 218 * Setting this to true requests that all computations be aborted by 219 * throwing AbortedException. Must be rest to false before any further 220 * computation. Ideally Thread.interrupt() should be used instead, but 221 * that doesn't appear to be consistently supported by browser VMs. 222 */ 223 public volatile static boolean please_stop = false; 224 225 /** 226 * Must be defined in subclasses of <TT>CR</tt>. 227 * Most users can ignore the existence of this method, and will 228 * not ever need to define a <TT>CR</tt> subclass. 229 * Returns value / 2 ** precision rounded to an integer. 230 * The error in the result is strictly < 1. 231 * Informally, approximate(n) gives a scaled approximation 232 * accurate to 2**n. 233 * Implementations may safely assume that precision is 234 * at least a factor of 8 away from overflow. 235 * Called only with the lock on the <TT>CR</tt> object 236 * already held. 237 */ approximate(int precision)238 protected abstract BigInteger approximate(int precision); 239 transient int min_prec; 240 // The smallest precision value with which the above 241 // has been called. 242 transient BigInteger max_appr; 243 // The scaled approximation corresponding to min_prec. 244 transient boolean appr_valid = false; 245 // min_prec and max_val are valid. 246 247 // Helper functions bound_log2(int n)248 static int bound_log2(int n) { 249 int abs_n = Math.abs(n); 250 return (int)Math.ceil(Math.log((double)(abs_n + 1))/Math.log(2.0)); 251 } 252 // Check that a precision is at least a factor of 8 away from 253 // overflowng the integer used to hold a precision spec. 254 // We generally perform this check early on, and then convince 255 // ourselves that none of the operations performed on precisions 256 // inside a function can generate an overflow. check_prec(int n)257 static void check_prec(int n) { 258 int high = n >> 28; 259 // if n is not in danger of overflowing, then the 4 high order 260 // bits should be identical. Thus high is either 0 or -1. 261 // The rest of this is to test for either of those in a way 262 // that should be as cheap as possible. 263 int high_shifted = n >> 29; 264 if (0 != (high ^ high_shifted)) { 265 throw new PrecisionOverflowException(); 266 } 267 } 268 269 /** 270 * The constructive real number corresponding to a 271 * <TT>BigInteger</tt>. 272 */ valueOf(BigInteger n)273 public static CR valueOf(BigInteger n) { 274 return new int_CR(n); 275 } 276 277 /** 278 * The constructive real number corresponding to a 279 * Java <TT>int</tt>. 280 */ valueOf(int n)281 public static CR valueOf(int n) { 282 return valueOf(BigInteger.valueOf(n)); 283 } 284 285 /** 286 * The constructive real number corresponding to a 287 * Java <TT>long</tt>. 288 */ valueOf(long n)289 public static CR valueOf(long n) { 290 return valueOf(BigInteger.valueOf(n)); 291 } 292 293 /** 294 * The constructive real number corresponding to a 295 * Java <TT>double</tt>. 296 * The result is undefined if argument is infinite or NaN. 297 */ valueOf(double n)298 public static CR valueOf(double n) { 299 if (Double.isNaN(n)) throw new ArithmeticException("Nan argument"); 300 if (Double.isInfinite(n)) { 301 throw new ArithmeticException("Infinite argument"); 302 } 303 boolean negative = (n < 0.0); 304 long bits = Double.doubleToLongBits(Math.abs(n)); 305 long mantissa = (bits & 0xfffffffffffffL); 306 int biased_exp = (int)(bits >> 52); 307 int exp = biased_exp - 1075; 308 if (biased_exp != 0) { 309 mantissa += (1L << 52); 310 } else { 311 mantissa <<= 1; 312 } 313 CR result = valueOf(mantissa).shiftLeft(exp); 314 if (negative) result = result.negate(); 315 return result; 316 } 317 318 /** 319 * The constructive real number corresponding to a 320 * Java <TT>float</tt>. 321 * The result is undefined if argument is infinite or NaN. 322 */ valueOf(float n)323 public static CR valueOf(float n) { 324 return valueOf((double) n); 325 } 326 327 public static CR ZERO = valueOf(0); 328 public static CR ONE = valueOf(1); 329 330 // Multiply k by 2**n. shift(BigInteger k, int n)331 static BigInteger shift(BigInteger k, int n) { 332 if (n == 0) return k; 333 if (n < 0) return k.shiftRight(-n); 334 return k.shiftLeft(n); 335 } 336 337 // Multiply by 2**n, rounding result scale(BigInteger k, int n)338 static BigInteger scale(BigInteger k, int n) { 339 if (n >= 0) { 340 return k.shiftLeft(n); 341 } else { 342 BigInteger adj_k = shift(k, n+1).add(big1); 343 return adj_k.shiftRight(1); 344 } 345 } 346 347 // Identical to approximate(), but maintain and update cache. 348 /** 349 * Returns value / 2 ** prec rounded to an integer. 350 * The error in the result is strictly < 1. 351 * Produces the same answer as <TT>approximate</tt>, but uses and 352 * maintains a cached approximation. 353 * Normally not overridden, and called only from <TT>approximate</tt> 354 * methods in subclasses. Not needed if the provided operations 355 * on constructive reals suffice. 356 */ get_appr(int precision)357 public synchronized BigInteger get_appr(int precision) { 358 check_prec(precision); 359 if (appr_valid && precision >= min_prec) { 360 return scale(max_appr, min_prec - precision); 361 } else { 362 BigInteger result = approximate(precision); 363 min_prec = precision; 364 max_appr = result; 365 appr_valid = true; 366 return result; 367 } 368 } 369 370 // Return the position of the msd. 371 // If x.msd() == n then 372 // 2**(n-1) < abs(x) < 2**(n+1) 373 // This initial version assumes that max_appr is valid 374 // and sufficiently removed from zero 375 // that the msd is determined. known_msd()376 int known_msd() { 377 int first_digit; 378 int length; 379 if (max_appr.signum() >= 0) { 380 length = max_appr.bitLength(); 381 } else { 382 length = max_appr.negate().bitLength(); 383 } 384 first_digit = min_prec + length - 1; 385 return first_digit; 386 } 387 388 // This version may return Integer.MIN_VALUE if the correct 389 // answer is < n. msd(int n)390 int msd(int n) { 391 if (!appr_valid || 392 max_appr.compareTo(big1) <= 0 393 && max_appr.compareTo(bigm1) >= 0) { 394 get_appr(n - 1); 395 if (max_appr.abs().compareTo(big1) <= 0) { 396 // msd could still be arbitrarily far to the right. 397 return Integer.MIN_VALUE; 398 } 399 } 400 return known_msd(); 401 } 402 403 404 // Functionally equivalent, but iteratively evaluates to higher 405 // precision. iter_msd(int n)406 int iter_msd(int n) 407 { 408 int prec = 0; 409 410 for (;prec > n + 30; prec = (prec * 3)/2 - 16) { 411 int msd = msd(prec); 412 if (msd != Integer.MIN_VALUE) return msd; 413 check_prec(prec); 414 if (Thread.interrupted() || please_stop) { 415 throw new AbortedException(); 416 } 417 } 418 return msd(n); 419 } 420 421 // This version returns a correct answer eventually, except 422 // that it loops forever (or throws an exception when the 423 // requested precision overflows) if this constructive real is zero. msd()424 int msd() { 425 return iter_msd(Integer.MIN_VALUE); 426 } 427 428 // A helper function for toString. 429 // Generate a String containing n zeroes. zeroes(int n)430 private static String zeroes(int n) { 431 char[] a = new char[n]; 432 for (int i = 0; i < n; ++i) { 433 a[i] = '0'; 434 } 435 return new String(a); 436 } 437 438 // Natural log of 2. Needed for some prescaling below. 439 // ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) simple_ln()440 CR simple_ln() { 441 return new prescaled_ln_CR(this.subtract(ONE)); 442 } 443 static CR ten_ninths = valueOf(10).divide(valueOf(9)); 444 static CR twentyfive_twentyfourths = valueOf(25).divide(valueOf(24)); 445 static CR eightyone_eightyeths = valueOf(81).divide(valueOf(80)); 446 static CR ln2_1 = valueOf(7).multiply(ten_ninths.simple_ln()); 447 static CR ln2_2 = 448 valueOf(2).multiply(twentyfive_twentyfourths.simple_ln()); 449 static CR ln2_3 = valueOf(3).multiply(eightyone_eightyeths.simple_ln()); 450 static CR ln2 = ln2_1.subtract(ln2_2).add(ln2_3); 451 452 // Atan of integer reciprocal. Used for atan_PI. Could perhaps be made 453 // public. atan_reciprocal(int n)454 static CR atan_reciprocal(int n) { 455 return new integral_atan_CR(n); 456 } 457 // Other constants used for PI computation. 458 static CR four = valueOf(4); 459 460 // Public operations. 461 /** 462 * Return 0 if x = y to within the indicated tolerance, 463 * -1 if x < y, and +1 if x > y. If x and y are indeed 464 * equal, it is guaranteed that 0 will be returned. If 465 * they differ by less than the tolerance, anything 466 * may happen. The tolerance allowed is 467 * the maximum of (abs(this)+abs(x))*(2**r) and 2**a 468 * @param x The other constructive real 469 * @param r Relative tolerance in bits 470 * @param a Absolute tolerance in bits 471 */ compareTo(CR x, int r, int a)472 public int compareTo(CR x, int r, int a) { 473 int this_msd = iter_msd(a); 474 int x_msd = x.iter_msd(this_msd > a? this_msd : a); 475 int max_msd = (x_msd > this_msd? x_msd : this_msd); 476 int rel = max_msd + r; 477 // This can't approach overflow, since r and a are 478 // effectively divided by 2, and msds are checked. 479 int abs_prec = (rel > a? rel : a); 480 return compareTo(x, abs_prec); 481 } 482 483 /** 484 * Approximate comparison with only an absolute tolerance. 485 * Identical to the three argument version, but without a relative 486 * tolerance. 487 * Result is 0 if both constructive reals are equal, indeterminate 488 * if they differ by less than 2**a. 489 * 490 * @param x The other constructive real 491 * @param a Absolute tolerance in bits 492 */ compareTo(CR x, int a)493 public int compareTo(CR x, int a) { 494 int needed_prec = a - 1; 495 BigInteger this_appr = get_appr(needed_prec); 496 BigInteger x_appr = x.get_appr(needed_prec); 497 int comp1 = this_appr.compareTo(x_appr.add(big1)); 498 if (comp1 > 0) return 1; 499 int comp2 = this_appr.compareTo(x_appr.subtract(big1)); 500 if (comp2 < 0) return -1; 501 return 0; 502 } 503 504 /** 505 * Return -1 if <TT>this < x</tt>, or +1 if <TT>this > x</tt>. 506 * Should be called only if <TT>this != x</tt>. 507 * If <TT>this == x</tt>, this will not terminate correctly; typically it 508 * will run until it exhausts memory. 509 * If the two constructive reals may be equal, the two or 3 argument 510 * version of compareTo should be used. 511 */ compareTo(CR x)512 public int compareTo(CR x) { 513 for (int a = -20; ; a *= 2) { 514 check_prec(a); 515 int result = compareTo(x, a); 516 if (0 != result) return result; 517 if (Thread.interrupted() || please_stop) { 518 throw new AbortedException(); 519 } 520 } 521 } 522 523 /** 524 * Equivalent to <TT>compareTo(CR.valueOf(0), a)</tt> 525 */ signum(int a)526 public int signum(int a) { 527 if (appr_valid) { 528 int quick_try = max_appr.signum(); 529 if (0 != quick_try) return quick_try; 530 } 531 int needed_prec = a - 1; 532 BigInteger this_appr = get_appr(needed_prec); 533 return this_appr.signum(); 534 } 535 536 /** 537 * Return -1 if negative, +1 if positive. 538 * Should be called only if <TT>this != 0</tt>. 539 * In the 0 case, this will not terminate correctly; typically it 540 * will run until it exhausts memory. 541 * If the two constructive reals may be equal, the one or two argument 542 * version of signum should be used. 543 */ signum()544 public int signum() { 545 for (int a = -20; ; a *= 2) { 546 check_prec(a); 547 int result = signum(a); 548 if (0 != result) return result; 549 if (Thread.interrupted() || please_stop) { 550 throw new AbortedException(); 551 } 552 } 553 } 554 555 /** 556 * Return the constructive real number corresponding to the given 557 * textual representation and radix. 558 * 559 * @param s [-] digit* [. digit*] 560 * @param radix 561 */ 562 valueOf(String s, int radix)563 public static CR valueOf(String s, int radix) 564 throws NumberFormatException { 565 int len = s.length(); 566 int start_pos = 0, point_pos; 567 String fraction; 568 while (s.charAt(start_pos) == ' ') ++start_pos; 569 while (s.charAt(len - 1) == ' ') --len; 570 point_pos = s.indexOf('.', start_pos); 571 if (point_pos == -1) { 572 point_pos = len; 573 fraction = "0"; 574 } else { 575 fraction = s.substring(point_pos + 1, len); 576 } 577 String whole = s.substring(start_pos, point_pos); 578 BigInteger scaled_result = new BigInteger(whole + fraction, radix); 579 BigInteger divisor = BigInteger.valueOf(radix).pow(fraction.length()); 580 return CR.valueOf(scaled_result).divide(CR.valueOf(divisor)); 581 } 582 583 /** 584 * Return a textual representation accurate to <TT>n</tt> places 585 * to the right of the decimal point. <TT>n</tt> must be nonnegative. 586 * 587 * @param n Number of digits (>= 0) included to the right of decimal point 588 * @param radix Base ( >= 2, <= 16) for the resulting representation. 589 */ toString(int n, int radix)590 public String toString(int n, int radix) { 591 CR scaled_CR; 592 if (16 == radix) { 593 scaled_CR = shiftLeft(4*n); 594 } else { 595 BigInteger scale_factor = BigInteger.valueOf(radix).pow(n); 596 scaled_CR = multiply(new int_CR(scale_factor)); 597 } 598 BigInteger scaled_int = scaled_CR.get_appr(0); 599 String scaled_string = scaled_int.abs().toString(radix); 600 String result; 601 if (0 == n) { 602 result = scaled_string; 603 } else { 604 int len = scaled_string.length(); 605 if (len <= n) { 606 // Add sufficient leading zeroes 607 String z = zeroes(n + 1 - len); 608 scaled_string = z + scaled_string; 609 len = n + 1; 610 } 611 String whole = scaled_string.substring(0, len - n); 612 String fraction = scaled_string.substring(len - n); 613 result = whole + "." + fraction; 614 } 615 if (scaled_int.signum() < 0) { 616 result = "-" + result; 617 } 618 return result; 619 } 620 621 622 /** 623 * Equivalent to <TT>toString(n,10)</tt> 624 * 625 * @param n Number of digits included to the right of decimal point 626 */ toString(int n)627 public String toString(int n) { 628 return toString(n, 10); 629 } 630 631 /** 632 * Equivalent to <TT>toString(10, 10)</tt> 633 */ toString()634 public String toString() { 635 return toString(10); 636 } 637 638 static double doubleLog2 = Math.log(2.0); 639 /** 640 * Return a textual scientific notation representation accurate 641 * to <TT>n</tt> places to the right of the decimal point. 642 * <TT>n</tt> must be nonnegative. A value smaller than 643 * <TT>radix</tt>**-<TT>m</tt> may be displayed as 0. 644 * The <TT>mantissa</tt> component of the result is either "0" 645 * or exactly <TT>n</tt> digits long. The <TT>sign</tt> 646 * component is zero exactly when the mantissa is "0". 647 * 648 * @param n Number of digits (> 0) included to the right of decimal point. 649 * @param radix Base ( ≥ 2, ≤ 16) for the resulting representation. 650 * @param m Precision used to distinguish number from zero. 651 * Expressed as a power of m. 652 */ toStringFloatRep(int n, int radix, int m)653 public StringFloatRep toStringFloatRep(int n, int radix, int m) { 654 if (n <= 0) throw new ArithmeticException("Bad precision argument"); 655 double log2_radix = Math.log((double)radix)/doubleLog2; 656 BigInteger big_radix = BigInteger.valueOf(radix); 657 long long_msd_prec = (long)(log2_radix * (double)m); 658 if (long_msd_prec > (long)Integer.MAX_VALUE 659 || long_msd_prec < (long)Integer.MIN_VALUE) 660 throw new PrecisionOverflowException(); 661 int msd_prec = (int)long_msd_prec; 662 check_prec(msd_prec); 663 int msd = iter_msd(msd_prec - 2); 664 if (msd == Integer.MIN_VALUE) 665 return new StringFloatRep(0, "0", radix, 0); 666 int exponent = (int)Math.ceil((double)msd / log2_radix); 667 // Guess for the exponent. Try to get it usually right. 668 int scale_exp = exponent - n; 669 CR scale; 670 if (scale_exp > 0) { 671 scale = CR.valueOf(big_radix.pow(scale_exp)).inverse(); 672 } else { 673 scale = CR.valueOf(big_radix.pow(-scale_exp)); 674 } 675 CR scaled_res = multiply(scale); 676 BigInteger scaled_int = scaled_res.get_appr(0); 677 int sign = scaled_int.signum(); 678 String scaled_string = scaled_int.abs().toString(radix); 679 while (scaled_string.length() < n) { 680 // exponent was too large. Adjust. 681 scaled_res = scaled_res.multiply(CR.valueOf(big_radix)); 682 exponent -= 1; 683 scaled_int = scaled_res.get_appr(0); 684 sign = scaled_int.signum(); 685 scaled_string = scaled_int.abs().toString(radix); 686 } 687 if (scaled_string.length() > n) { 688 // exponent was too small. Adjust by truncating. 689 exponent += (scaled_string.length() - n); 690 scaled_string = scaled_string.substring(0, n); 691 } 692 return new StringFloatRep(sign, scaled_string, radix, exponent); 693 } 694 695 /** 696 * Return a BigInteger which differs by less than one from the 697 * constructive real. 698 */ BigIntegerValue()699 public BigInteger BigIntegerValue() { 700 return get_appr(0); 701 } 702 703 /** 704 * Return an int which differs by less than one from the 705 * constructive real. Behavior on overflow is undefined. 706 */ intValue()707 public int intValue() { 708 return BigIntegerValue().intValue(); 709 } 710 711 /** 712 * Return an int which differs by less than one from the 713 * constructive real. Behavior on overflow is undefined. 714 */ byteValue()715 public byte byteValue() { 716 return BigIntegerValue().byteValue(); 717 } 718 719 /** 720 * Return a long which differs by less than one from the 721 * constructive real. Behavior on overflow is undefined. 722 */ longValue()723 public long longValue() { 724 return BigIntegerValue().longValue(); 725 } 726 727 /** 728 * Return a double which differs by less than one in the least 729 * represented bit from the constructive real. 730 * (We're in fact closer to round-to-nearest than that, but we can't and 731 * don't promise correct rounding.) 732 */ doubleValue()733 public double doubleValue() { 734 int my_msd = iter_msd(-1080 /* slightly > exp. range */); 735 if (Integer.MIN_VALUE == my_msd) return 0.0; 736 int needed_prec = my_msd - 60; 737 double scaled_int = get_appr(needed_prec).doubleValue(); 738 boolean may_underflow = (needed_prec < -1000); 739 long scaled_int_rep = Double.doubleToLongBits(scaled_int); 740 long exp_adj = may_underflow? needed_prec + 96 : needed_prec; 741 long orig_exp = (scaled_int_rep >> 52) & 0x7ff; 742 if (((orig_exp + exp_adj) & ~0x7ff) != 0) { 743 // Original unbiased exponent is > 50. Exp_adj > -1050. 744 // Thus this can overflow the 11 bit exponent only if the result 745 // itself overflows. 746 if (scaled_int < 0.0) { 747 return Double.NEGATIVE_INFINITY; 748 } else { 749 return Double.POSITIVE_INFINITY; 750 } 751 } 752 scaled_int_rep += exp_adj << 52; 753 double result = Double.longBitsToDouble(scaled_int_rep); 754 if (may_underflow) { 755 double two48 = (double)(1L << 48); 756 return result/two48/two48; 757 } else { 758 return result; 759 } 760 } 761 762 /** 763 * Return a float which differs by less than one in the least 764 * represented bit from the constructive real. 765 */ floatValue()766 public float floatValue() { 767 return (float)doubleValue(); 768 // Note that double-rounding is not a problem here, since we 769 // cannot, and do not, guarantee correct rounding. 770 } 771 772 /** 773 * Add two constructive reals. 774 */ add(CR x)775 public CR add(CR x) { 776 return new add_CR(this, x); 777 } 778 779 /** 780 * Multiply a constructive real by 2**n. 781 * @param n shift count, may be negative 782 */ shiftLeft(int n)783 public CR shiftLeft(int n) { 784 check_prec(n); 785 return new shifted_CR(this, n); 786 } 787 788 /** 789 * Multiply a constructive real by 2**(-n). 790 * @param n shift count, may be negative 791 */ shiftRight(int n)792 public CR shiftRight(int n) { 793 check_prec(n); 794 return new shifted_CR(this, -n); 795 } 796 797 /** 798 * Produce a constructive real equivalent to the original, assuming 799 * the original was an integer. Undefined results if the original 800 * was not an integer. Prevents evaluation of digits to the right 801 * of the decimal point, and may thus improve performance. 802 */ assumeInt()803 public CR assumeInt() { 804 return new assumed_int_CR(this); 805 } 806 807 /** 808 * The additive inverse of a constructive real 809 */ negate()810 public CR negate() { 811 return new neg_CR(this); 812 } 813 814 /** 815 * The difference between two constructive reals 816 */ subtract(CR x)817 public CR subtract(CR x) { 818 return new add_CR(this, x.negate()); 819 } 820 821 /** 822 * The product of two constructive reals 823 */ multiply(CR x)824 public CR multiply(CR x) { 825 return new mult_CR(this, x); 826 } 827 828 /** 829 * The multiplicative inverse of a constructive real. 830 * <TT>x.inverse()</tt> is equivalent to <TT>CR.valueOf(1).divide(x)</tt>. 831 */ inverse()832 public CR inverse() { 833 return new inv_CR(this); 834 } 835 836 /** 837 * The quotient of two constructive reals. 838 */ divide(CR x)839 public CR divide(CR x) { 840 return new mult_CR(this, x.inverse()); 841 } 842 843 /** 844 * The real number <TT>x</tt> if <TT>this</tt> < 0, or <TT>y</tt> otherwise. 845 * Requires <TT>x</tt> = <TT>y</tt> if <TT>this</tt> = 0. 846 * Since comparisons may diverge, this is often 847 * a useful alternative to conditionals. 848 */ select(CR x, CR y)849 public CR select(CR x, CR y) { 850 return new select_CR(this, x, y); 851 } 852 853 /** 854 * The maximum of two constructive reals. 855 */ max(CR x)856 public CR max(CR x) { 857 return subtract(x).select(x, this); 858 } 859 860 /** 861 * The minimum of two constructive reals. 862 */ min(CR x)863 public CR min(CR x) { 864 return subtract(x).select(this, x); 865 } 866 867 /** 868 * The absolute value of a constructive reals. 869 * Note that this cannot be written as a conditional. 870 */ abs()871 public CR abs() { 872 return select(negate(), this); 873 } 874 875 /** 876 * The exponential function, that is e**<TT>this</tt>. 877 */ exp()878 public CR exp() { 879 final int low_prec = -10; 880 BigInteger rough_appr = get_appr(low_prec); 881 // Handle negative arguments directly; negating and computing inverse 882 // can be very expensive. 883 if (rough_appr.compareTo(big2) > 0 || rough_appr.compareTo(bigm2) < 0) { 884 CR square_root = shiftRight(1).exp(); 885 return square_root.multiply(square_root); 886 } else { 887 return new prescaled_exp_CR(this); 888 } 889 } 890 891 /** 892 * The ratio of a circle's circumference to its diameter. 893 */ 894 public static CR PI = new gl_pi_CR(); 895 896 // Our old PI implementation. Keep this around for now to allow checking. 897 // This implementation may also be faster for BigInteger implementations 898 // that support only quadratic multiplication, but exhibit high performance 899 // for small computations. (The standard Android 6 implementation supports 900 // subquadratic multiplication, but has high constant overhead.) Many other 901 // atan-based formulas are possible, but based on superficial 902 // experimentation, this is roughly as good as the more complex formulas. 903 public static CR atan_PI = four.multiply(four.multiply(atan_reciprocal(5)) 904 .subtract(atan_reciprocal(239))); 905 // pi/4 = 4*atan(1/5) - atan(1/239) 906 static CR half_pi = PI.shiftRight(1); 907 908 /** 909 * The trigonometric cosine function. 910 */ cos()911 public CR cos() { 912 BigInteger halfpi_multiples = divide(PI).get_appr(-1); 913 BigInteger abs_halfpi_multiples = halfpi_multiples.abs(); 914 if (abs_halfpi_multiples.compareTo(big2) >= 0) { 915 // Subtract multiples of PI 916 BigInteger pi_multiples = scale(halfpi_multiples, -1); 917 CR adjustment = PI.multiply(CR.valueOf(pi_multiples)); 918 if (pi_multiples.and(big1).signum() != 0) { 919 return subtract(adjustment).cos().negate(); 920 } else { 921 return subtract(adjustment).cos(); 922 } 923 } else if (get_appr(-1).abs().compareTo(big2) >= 0) { 924 // Scale further with double angle formula 925 CR cos_half = shiftRight(1).cos(); 926 return cos_half.multiply(cos_half).shiftLeft(1).subtract(ONE); 927 } else { 928 return new prescaled_cos_CR(this); 929 } 930 } 931 932 /** 933 * The trigonometric sine function. 934 */ sin()935 public CR sin() { 936 return half_pi.subtract(this).cos(); 937 } 938 939 /** 940 * The trignonometric arc (inverse) sine function. 941 */ asin()942 public CR asin() { 943 BigInteger rough_appr = get_appr(-10); 944 if (rough_appr.compareTo(big750) /* 1/sqrt(2) + a bit */ > 0){ 945 CR new_arg = ONE.subtract(multiply(this)).sqrt(); 946 return new_arg.acos(); 947 } else if (rough_appr.compareTo(bigm750) < 0) { 948 return negate().asin().negate(); 949 } else { 950 return new prescaled_asin_CR(this); 951 } 952 } 953 954 /** 955 * The trignonometric arc (inverse) cosine function. 956 */ acos()957 public CR acos() { 958 return half_pi.subtract(asin()); 959 } 960 961 static final BigInteger low_ln_limit = big8; /* sixteenths, i.e. 1/2 */ 962 static final BigInteger high_ln_limit = 963 BigInteger.valueOf(16 + 8 /* 1.5 */); 964 static final BigInteger scaled_4 = 965 BigInteger.valueOf(4*16); 966 967 /** 968 * The natural (base e) logarithm. 969 */ ln()970 public CR ln() { 971 final int low_prec = -4; 972 BigInteger rough_appr = get_appr(low_prec); /* In sixteenths */ 973 if (rough_appr.compareTo(big0) < 0) { 974 throw new ArithmeticException("ln(negative)"); 975 } 976 if (rough_appr.compareTo(low_ln_limit) <= 0) { 977 return inverse().ln().negate(); 978 } 979 if (rough_appr.compareTo(high_ln_limit) >= 0) { 980 if (rough_appr.compareTo(scaled_4) <= 0) { 981 CR quarter = sqrt().sqrt().ln(); 982 return quarter.shiftLeft(2); 983 } else { 984 int extra_bits = rough_appr.bitLength() - 3; 985 CR scaled_result = shiftRight(extra_bits).ln(); 986 return scaled_result.add(CR.valueOf(extra_bits).multiply(ln2)); 987 } 988 } 989 return simple_ln(); 990 } 991 992 /** 993 * The square root of a constructive real. 994 */ sqrt()995 public CR sqrt() { 996 return new sqrt_CR(this); 997 } 998 999 } // end of CR 1000 1001 1002 // 1003 // A specialization of CR for cases in which approximate() calls 1004 // to increase evaluation precision are somewhat expensive. 1005 // If we need to (re)evaluate, we speculatively evaluate to slightly 1006 // higher precision, miminimizing reevaluations. 1007 // Note that this requires any arguments to be evaluated to higher 1008 // precision than absolutely necessary. It can thus potentially 1009 // result in lots of wasted effort, and should be used judiciously. 1010 // This assumes that the order of magnitude of the number is roughly one. 1011 // 1012 abstract class slow_CR extends CR { 1013 static int max_prec = -64; 1014 static int prec_incr = 32; get_appr(int precision)1015 public synchronized BigInteger get_appr(int precision) { 1016 check_prec(precision); 1017 if (appr_valid && precision >= min_prec) { 1018 return scale(max_appr, min_prec - precision); 1019 } else { 1020 int eval_prec = (precision >= max_prec? max_prec : 1021 (precision - prec_incr + 1) & ~(prec_incr - 1)); 1022 BigInteger result = approximate(eval_prec); 1023 min_prec = eval_prec; 1024 max_appr = result; 1025 appr_valid = true; 1026 return scale(result, eval_prec - precision); 1027 } 1028 } 1029 } 1030 1031 1032 // Representation of an integer constant. Private. 1033 class int_CR extends CR { 1034 BigInteger value; int_CR(BigInteger n)1035 int_CR(BigInteger n) { 1036 value = n; 1037 } approximate(int p)1038 protected BigInteger approximate(int p) { 1039 return scale(value, -p) ; 1040 } 1041 } 1042 1043 // Representation of a number that may not have been completely 1044 // evaluated, but is assumed to be an integer. Hence we never 1045 // evaluate beyond the decimal point. 1046 class assumed_int_CR extends CR { 1047 CR value; assumed_int_CR(CR x)1048 assumed_int_CR(CR x) { 1049 value = x; 1050 } approximate(int p)1051 protected BigInteger approximate(int p) { 1052 if (p >= 0) { 1053 return value.get_appr(p); 1054 } else { 1055 return scale(value.get_appr(0), -p) ; 1056 } 1057 } 1058 } 1059 1060 // Representation of the sum of 2 constructive reals. Private. 1061 class add_CR extends CR { 1062 CR op1; 1063 CR op2; add_CR(CR x, CR y)1064 add_CR(CR x, CR y) { 1065 op1 = x; 1066 op2 = y; 1067 } approximate(int p)1068 protected BigInteger approximate(int p) { 1069 // Args need to be evaluated so that each error is < 1/4 ulp. 1070 // Rounding error from the cale call is <= 1/2 ulp, so that 1071 // final error is < 1 ulp. 1072 return scale(op1.get_appr(p-2).add(op2.get_appr(p-2)), -2); 1073 } 1074 } 1075 1076 // Representation of a CR multiplied by 2**n 1077 class shifted_CR extends CR { 1078 CR op; 1079 int count; shifted_CR(CR x, int n)1080 shifted_CR(CR x, int n) { 1081 op = x; 1082 count = n; 1083 } approximate(int p)1084 protected BigInteger approximate(int p) { 1085 return op.get_appr(p - count); 1086 } 1087 } 1088 1089 // Representation of the negation of a constructive real. Private. 1090 class neg_CR extends CR { 1091 CR op; neg_CR(CR x)1092 neg_CR(CR x) { 1093 op = x; 1094 } approximate(int p)1095 protected BigInteger approximate(int p) { 1096 return op.get_appr(p).negate(); 1097 } 1098 } 1099 1100 // Representation of: 1101 // op1 if selector < 0 1102 // op2 if selector >= 0 1103 // Assumes x = y if s = 0 1104 class select_CR extends CR { 1105 CR selector; 1106 int selector_sign; 1107 CR op1; 1108 CR op2; select_CR(CR s, CR x, CR y)1109 select_CR(CR s, CR x, CR y) { 1110 selector = s; 1111 int selector_sign = selector.get_appr(-20).signum(); 1112 op1 = x; 1113 op2 = y; 1114 } approximate(int p)1115 protected BigInteger approximate(int p) { 1116 if (selector_sign < 0) return op1.get_appr(p); 1117 if (selector_sign > 0) return op2.get_appr(p); 1118 BigInteger op1_appr = op1.get_appr(p-1); 1119 BigInteger op2_appr = op2.get_appr(p-1); 1120 BigInteger diff = op1_appr.subtract(op2_appr).abs(); 1121 if (diff.compareTo(big1) <= 0) { 1122 // close enough; use either 1123 return scale(op1_appr, -1); 1124 } 1125 // op1 and op2 are different; selector != 0; 1126 // safe to get sign of selector. 1127 if (selector.signum() < 0) { 1128 selector_sign = -1; 1129 return scale(op1_appr, -1); 1130 } else { 1131 selector_sign = 1; 1132 return scale(op2_appr, -1); 1133 } 1134 } 1135 } 1136 1137 // Representation of the product of 2 constructive reals. Private. 1138 class mult_CR extends CR { 1139 CR op1; 1140 CR op2; mult_CR(CR x, CR y)1141 mult_CR(CR x, CR y) { 1142 op1 = x; 1143 op2 = y; 1144 } approximate(int p)1145 protected BigInteger approximate(int p) { 1146 int half_prec = (p >> 1) - 1; 1147 int msd_op1 = op1.msd(half_prec); 1148 int msd_op2; 1149 1150 if (msd_op1 == Integer.MIN_VALUE) { 1151 msd_op2 = op2.msd(half_prec); 1152 if (msd_op2 == Integer.MIN_VALUE) { 1153 // Product is small enough that zero will do as an 1154 // approximation. 1155 return big0; 1156 } else { 1157 // Swap them, so the larger operand (in absolute value) 1158 // is first. 1159 CR tmp; 1160 tmp = op1; 1161 op1 = op2; 1162 op2 = tmp; 1163 msd_op1 = msd_op2; 1164 } 1165 } 1166 // msd_op1 is valid at this point. 1167 int prec2 = p - msd_op1 - 3; // Precision needed for op2. 1168 // The appr. error is multiplied by at most 1169 // 2 ** (msd_op1 + 1) 1170 // Thus each approximation contributes 1/4 ulp 1171 // to the rounding error, and the final rounding adds 1172 // another 1/2 ulp. 1173 BigInteger appr2 = op2.get_appr(prec2); 1174 if (appr2.signum() == 0) return big0; 1175 msd_op2 = op2.known_msd(); 1176 int prec1 = p - msd_op2 - 3; // Precision needed for op1. 1177 BigInteger appr1 = op1.get_appr(prec1); 1178 int scale_digits = prec1 + prec2 - p; 1179 return scale(appr1.multiply(appr2), scale_digits); 1180 } 1181 } 1182 1183 // Representation of the multiplicative inverse of a constructive 1184 // real. Private. Should use Newton iteration to refine estimates. 1185 class inv_CR extends CR { 1186 CR op; inv_CR(CR x)1187 inv_CR(CR x) { op = x; } approximate(int p)1188 protected BigInteger approximate(int p) { 1189 int msd = op.msd(); 1190 int inv_msd = 1 - msd; 1191 int digits_needed = inv_msd - p + 3; 1192 // Number of SIGNIFICANT digits needed for 1193 // argument, excl. msd position, which may 1194 // be fictitious, since msd routine can be 1195 // off by 1. Roughly 1 extra digit is 1196 // needed since the relative error is the 1197 // same in the argument and result, but 1198 // this isn't quite the same as the number 1199 // of significant digits. Another digit 1200 // is needed to compensate for slop in the 1201 // calculation. 1202 // One further bit is required, since the 1203 // final rounding introduces a 0.5 ulp 1204 // error. 1205 int prec_needed = msd - digits_needed; 1206 int log_scale_factor = -p - prec_needed; 1207 if (log_scale_factor < 0) return big0; 1208 BigInteger dividend = big1.shiftLeft(log_scale_factor); 1209 BigInteger scaled_divisor = op.get_appr(prec_needed); 1210 BigInteger abs_scaled_divisor = scaled_divisor.abs(); 1211 BigInteger adj_dividend = dividend.add( 1212 abs_scaled_divisor.shiftRight(1)); 1213 // Adjustment so that final result is rounded. 1214 BigInteger result = adj_dividend.divide(abs_scaled_divisor); 1215 if (scaled_divisor.signum() < 0) { 1216 return result.negate(); 1217 } else { 1218 return result; 1219 } 1220 } 1221 } 1222 1223 1224 // Representation of the exponential of a constructive real. Private. 1225 // Uses a Taylor series expansion. Assumes |x| < 1/2. 1226 // Note: this is known to be a bad algorithm for 1227 // floating point. Unfortunately, other alternatives 1228 // appear to require precomputed information. 1229 class prescaled_exp_CR extends CR { 1230 CR op; prescaled_exp_CR(CR x)1231 prescaled_exp_CR(CR x) { op = x; } approximate(int p)1232 protected BigInteger approximate(int p) { 1233 if (p >= 1) return big0; 1234 int iterations_needed = -p/2 + 2; // conservative estimate > 0. 1235 // Claim: each intermediate term is accurate 1236 // to 2*2^calc_precision. 1237 // Total rounding error in series computation is 1238 // 2*iterations_needed*2^calc_precision, 1239 // exclusive of error in op. 1240 int calc_precision = p - bound_log2(2*iterations_needed) 1241 - 4; // for error in op, truncation. 1242 int op_prec = p - 3; 1243 BigInteger op_appr = op.get_appr(op_prec); 1244 // Error in argument results in error of < 3/8 ulp. 1245 // Sum of term eval. rounding error is < 1/16 ulp. 1246 // Series truncation error < 1/16 ulp. 1247 // Final rounding error is <= 1/2 ulp. 1248 // Thus final error is < 1 ulp. 1249 BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 1250 BigInteger current_term = scaled_1; 1251 BigInteger current_sum = scaled_1; 1252 int n = 0; 1253 BigInteger max_trunc_error = 1254 big1.shiftLeft(p - 4 - calc_precision); 1255 while (current_term.abs().compareTo(max_trunc_error) >= 0) { 1256 if (Thread.interrupted() || please_stop) throw new AbortedException(); 1257 n += 1; 1258 /* current_term = current_term * op / n */ 1259 current_term = scale(current_term.multiply(op_appr), op_prec); 1260 current_term = current_term.divide(BigInteger.valueOf(n)); 1261 current_sum = current_sum.add(current_term); 1262 } 1263 return scale(current_sum, calc_precision - p); 1264 } 1265 } 1266 1267 // Representation of the cosine of a constructive real. Private. 1268 // Uses a Taylor series expansion. Assumes |x| < 1. 1269 class prescaled_cos_CR extends slow_CR { 1270 CR op; prescaled_cos_CR(CR x)1271 prescaled_cos_CR(CR x) { 1272 op = x; 1273 } approximate(int p)1274 protected BigInteger approximate(int p) { 1275 if (p >= 1) return big0; 1276 int iterations_needed = -p/2 + 4; // conservative estimate > 0. 1277 // Claim: each intermediate term is accurate 1278 // to 2*2^calc_precision. 1279 // Total rounding error in series computation is 1280 // 2*iterations_needed*2^calc_precision, 1281 // exclusive of error in op. 1282 int calc_precision = p - bound_log2(2*iterations_needed) 1283 - 4; // for error in op, truncation. 1284 int op_prec = p - 2; 1285 BigInteger op_appr = op.get_appr(op_prec); 1286 // Error in argument results in error of < 1/4 ulp. 1287 // Cumulative arithmetic rounding error is < 1/16 ulp. 1288 // Series truncation error < 1/16 ulp. 1289 // Final rounding error is <= 1/2 ulp. 1290 // Thus final error is < 1 ulp. 1291 BigInteger current_term; 1292 int n; 1293 BigInteger max_trunc_error = 1294 big1.shiftLeft(p - 4 - calc_precision); 1295 n = 0; 1296 current_term = big1.shiftLeft(-calc_precision); 1297 BigInteger current_sum = current_term; 1298 while (current_term.abs().compareTo(max_trunc_error) >= 0) { 1299 if (Thread.interrupted() || please_stop) throw new AbortedException(); 1300 n += 2; 1301 /* current_term = - current_term * op * op / n * (n - 1) */ 1302 current_term = scale(current_term.multiply(op_appr), op_prec); 1303 current_term = scale(current_term.multiply(op_appr), op_prec); 1304 BigInteger divisor = BigInteger.valueOf(-n) 1305 .multiply(BigInteger.valueOf(n-1)); 1306 current_term = current_term.divide(divisor); 1307 current_sum = current_sum.add(current_term); 1308 } 1309 return scale(current_sum, calc_precision - p); 1310 } 1311 } 1312 1313 // The constructive real atan(1/n), where n is a small integer 1314 // > base. 1315 // This gives a simple and moderately fast way to compute PI. 1316 class integral_atan_CR extends slow_CR { 1317 int op; integral_atan_CR(int x)1318 integral_atan_CR(int x) { op = x; } approximate(int p)1319 protected BigInteger approximate(int p) { 1320 if (p >= 1) return big0; 1321 int iterations_needed = -p/2 + 2; // conservative estimate > 0. 1322 // Claim: each intermediate term is accurate 1323 // to 2*base^calc_precision. 1324 // Total rounding error in series computation is 1325 // 2*iterations_needed*base^calc_precision, 1326 // exclusive of error in op. 1327 int calc_precision = p - bound_log2(2*iterations_needed) 1328 - 2; // for error in op, truncation. 1329 // Error in argument results in error of < 3/8 ulp. 1330 // Cumulative arithmetic rounding error is < 1/4 ulp. 1331 // Series truncation error < 1/4 ulp. 1332 // Final rounding error is <= 1/2 ulp. 1333 // Thus final error is < 1 ulp. 1334 BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 1335 BigInteger big_op = BigInteger.valueOf(op); 1336 BigInteger big_op_squared = BigInteger.valueOf(op*op); 1337 BigInteger op_inverse = scaled_1.divide(big_op); 1338 BigInteger current_power = op_inverse; 1339 BigInteger current_term = op_inverse; 1340 BigInteger current_sum = op_inverse; 1341 int current_sign = 1; 1342 int n = 1; 1343 BigInteger max_trunc_error = 1344 big1.shiftLeft(p - 2 - calc_precision); 1345 while (current_term.abs().compareTo(max_trunc_error) >= 0) { 1346 if (Thread.interrupted() || please_stop) throw new AbortedException(); 1347 n += 2; 1348 current_power = current_power.divide(big_op_squared); 1349 current_sign = -current_sign; 1350 current_term = 1351 current_power.divide(BigInteger.valueOf(current_sign*n)); 1352 current_sum = current_sum.add(current_term); 1353 } 1354 return scale(current_sum, calc_precision - p); 1355 } 1356 } 1357 1358 // Representation for ln(1 + op) 1359 class prescaled_ln_CR extends slow_CR { 1360 CR op; prescaled_ln_CR(CR x)1361 prescaled_ln_CR(CR x) { op = x; } 1362 // Compute an approximation of ln(1+x) to precision 1363 // prec. This assumes |x| < 1/2. 1364 // It uses a Taylor series expansion. 1365 // Unfortunately there appears to be no way to take 1366 // advantage of old information. 1367 // Note: this is known to be a bad algorithm for 1368 // floating point. Unfortunately, other alternatives 1369 // appear to require precomputed tabular information. approximate(int p)1370 protected BigInteger approximate(int p) { 1371 if (p >= 0) return big0; 1372 int iterations_needed = -p; // conservative estimate > 0. 1373 // Claim: each intermediate term is accurate 1374 // to 2*2^calc_precision. Total error is 1375 // 2*iterations_needed*2^calc_precision 1376 // exclusive of error in op. 1377 int calc_precision = p - bound_log2(2*iterations_needed) 1378 - 4; // for error in op, truncation. 1379 int op_prec = p - 3; 1380 BigInteger op_appr = op.get_appr(op_prec); 1381 // Error analysis as for exponential. 1382 BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 1383 BigInteger x_nth = scale(op_appr, op_prec - calc_precision); 1384 BigInteger current_term = x_nth; // x**n 1385 BigInteger current_sum = current_term; 1386 int n = 1; 1387 int current_sign = 1; // (-1)^(n-1) 1388 BigInteger max_trunc_error = 1389 big1.shiftLeft(p - 4 - calc_precision); 1390 while (current_term.abs().compareTo(max_trunc_error) >= 0) { 1391 if (Thread.interrupted() || please_stop) throw new AbortedException(); 1392 n += 1; 1393 current_sign = -current_sign; 1394 x_nth = scale(x_nth.multiply(op_appr), op_prec); 1395 current_term = x_nth.divide(BigInteger.valueOf(n * current_sign)); 1396 // x**n / (n * (-1)**(n-1)) 1397 current_sum = current_sum.add(current_term); 1398 } 1399 return scale(current_sum, calc_precision - p); 1400 } 1401 } 1402 1403 // Representation of the arcsine of a constructive real. Private. 1404 // Uses a Taylor series expansion. Assumes |x| < (1/2)^(1/3). 1405 class prescaled_asin_CR extends slow_CR { 1406 CR op; prescaled_asin_CR(CR x)1407 prescaled_asin_CR(CR x) { 1408 op = x; 1409 } approximate(int p)1410 protected BigInteger approximate(int p) { 1411 // The Taylor series is the sum of x^(2n+1) * (2n)!/(4^n n!^2 (2n+1)) 1412 // Note that (2n)!/(4^n n!^2) is always less than one. 1413 // (The denominator is effectively 2n*2n*(2n-2)*(2n-2)*...*2*2 1414 // which is clearly > (2n)!) 1415 // Thus all terms are bounded by x^(2n+1). 1416 // Unfortunately, there's no easy way to prescale the argument 1417 // to less than 1/sqrt(2), and we can only approximate that. 1418 // Thus the worst case iteration count is fairly high. 1419 // But it doesn't make much difference. 1420 if (p >= 2) return big0; // Never bigger than 4. 1421 int iterations_needed = -3 * p / 2 + 4; 1422 // conservative estimate > 0. 1423 // Follows from assumed bound on x and 1424 // the fact that only every other Taylor 1425 // Series term is present. 1426 // Claim: each intermediate term is accurate 1427 // to 2*2^calc_precision. 1428 // Total rounding error in series computation is 1429 // 2*iterations_needed*2^calc_precision, 1430 // exclusive of error in op. 1431 int calc_precision = p - bound_log2(2*iterations_needed) 1432 - 4; // for error in op, truncation. 1433 int op_prec = p - 3; // always <= -2 1434 BigInteger op_appr = op.get_appr(op_prec); 1435 // Error in argument results in error of < 1/4 ulp. 1436 // (Derivative is bounded by 2 in the specified range and we use 1437 // 3 extra digits.) 1438 // Ignoring the argument error, each term has an error of 1439 // < 3ulps relative to calc_precision, which is more precise than p. 1440 // Cumulative arithmetic rounding error is < 3/16 ulp (relative to p). 1441 // Series truncation error < 2/16 ulp. (Each computed term 1442 // is at most 2/3 of last one, so some of remaining series < 1443 // 3/2 * current term.) 1444 // Final rounding error is <= 1/2 ulp. 1445 // Thus final error is < 1 ulp (relative to p). 1446 BigInteger max_last_term = 1447 big1.shiftLeft(p - 4 - calc_precision); 1448 int exp = 1; // Current exponent, = 2n+1 in above expression 1449 BigInteger current_term = op_appr.shiftLeft(op_prec - calc_precision); 1450 BigInteger current_sum = current_term; 1451 BigInteger current_factor = current_term; 1452 // Current scaled Taylor series term 1453 // before division by the exponent. 1454 // Accurate to 3 ulp at calc_precision. 1455 while (current_term.abs().compareTo(max_last_term) >= 0) { 1456 if (Thread.interrupted() || please_stop) throw new AbortedException(); 1457 exp += 2; 1458 // current_factor = current_factor * op * op * (exp-1) * (exp-2) / 1459 // (exp-1) * (exp-1), with the two exp-1 factors cancelling, 1460 // giving 1461 // current_factor = current_factor * op * op * (exp-2) / (exp-1) 1462 // Thus the error any in the previous term is multiplied by 1463 // op^2, adding an error of < (1/2)^(2/3) < 2/3 the original 1464 // error. 1465 current_factor = current_factor.multiply(BigInteger.valueOf(exp - 2)); 1466 current_factor = scale(current_factor.multiply(op_appr), op_prec + 2); 1467 // Carry 2 extra bits of precision forward; thus 1468 // this effectively introduces 1/8 ulp error. 1469 current_factor = current_factor.multiply(op_appr); 1470 BigInteger divisor = BigInteger.valueOf(exp - 1); 1471 current_factor = current_factor.divide(divisor); 1472 // Another 1/4 ulp error here. 1473 current_factor = scale(current_factor, op_prec - 2); 1474 // Remove extra 2 bits. 1/2 ulp rounding error. 1475 // Current_factor has original 3 ulp rounding error, which we 1476 // reduced by 1, plus < 1 ulp new rounding error. 1477 current_term = current_factor.divide(BigInteger.valueOf(exp)); 1478 // Contributes 1 ulp error to sum plus at most 3 ulp 1479 // from current_factor. 1480 current_sum = current_sum.add(current_term); 1481 } 1482 return scale(current_sum, calc_precision - p); 1483 } 1484 } 1485 1486 1487 class sqrt_CR extends CR { 1488 CR op; sqrt_CR(CR x)1489 sqrt_CR(CR x) { op = x; } 1490 // Explicitly provide an initial approximation. 1491 // Useful for arithmetic geometric mean algorithms, where we've previously 1492 // computed a very similar square root. sqrt_CR(CR x, int min_p, BigInteger max_a)1493 sqrt_CR(CR x, int min_p, BigInteger max_a) { 1494 op = x; 1495 min_prec = min_p; 1496 max_appr = max_a; 1497 appr_valid = true; 1498 } 1499 final int fp_prec = 50; // Conservative estimate of number of 1500 // significant bits in double precision 1501 // computation. 1502 final int fp_op_prec = 60; approximate(int p)1503 protected BigInteger approximate(int p) { 1504 int max_op_prec_needed = 2*p - 1; 1505 int msd = op.iter_msd(max_op_prec_needed); 1506 if (msd <= max_op_prec_needed) return big0; 1507 int result_msd = msd/2; // +- 1 1508 int result_digits = result_msd - p; // +- 2 1509 if (result_digits > fp_prec) { 1510 // Compute less precise approximation and use a Newton iter. 1511 int appr_digits = result_digits/2 + 6; 1512 // This should be conservative. Is fewer enough? 1513 int appr_prec = result_msd - appr_digits; 1514 int prod_prec = 2*appr_prec; 1515 // First compute the argument to maximal precision, so we don't end up 1516 // reevaluating it incrementally. 1517 BigInteger op_appr = op.get_appr(prod_prec); 1518 BigInteger last_appr = get_appr(appr_prec); 1519 // Compute (last_appr * last_appr + op_appr) / last_appr / 2 1520 // while adjusting the scaling to make everything work 1521 BigInteger prod_prec_scaled_numerator = 1522 last_appr.multiply(last_appr).add(op_appr); 1523 BigInteger scaled_numerator = 1524 scale(prod_prec_scaled_numerator, appr_prec - p); 1525 BigInteger shifted_result = scaled_numerator.divide(last_appr); 1526 return shifted_result.add(big1).shiftRight(1); 1527 } else { 1528 // Use a double precision floating point approximation. 1529 // Make sure all precisions are even 1530 int op_prec = (msd - fp_op_prec) & ~1; 1531 int working_prec = op_prec - fp_op_prec; 1532 BigInteger scaled_bi_appr = op.get_appr(op_prec) 1533 .shiftLeft(fp_op_prec); 1534 double scaled_appr = scaled_bi_appr.doubleValue(); 1535 if (scaled_appr < 0.0) 1536 throw new ArithmeticException("sqrt(negative)"); 1537 double scaled_fp_sqrt = Math.sqrt(scaled_appr); 1538 BigInteger scaled_sqrt = BigInteger.valueOf((long)scaled_fp_sqrt); 1539 int shift_count = working_prec/2 - p; 1540 return shift(scaled_sqrt, shift_count); 1541 } 1542 } 1543 } 1544 1545 // The constant PI, computed using the Gauss-Legendre alternating 1546 // arithmetic-geometric mean algorithm: 1547 // a[0] = 1 1548 // b[0] = 1/sqrt(2) 1549 // t[0] = 1/4 1550 // p[0] = 1 1551 // 1552 // a[n+1] = (a[n] + b[n])/2 (arithmetic mean, between 0.8 and 1) 1553 // b[n+1] = sqrt(a[n] * b[n]) (geometric mean, between 0.7 and 1) 1554 // t[n+1] = t[n] - (2^n)(a[n]-a[n+1])^2, (always between 0.2 and 0.25) 1555 // 1556 // pi is then approximated as (a[n+1]+b[n+1])^2 / 4*t[n+1]. 1557 // 1558 class gl_pi_CR extends slow_CR { 1559 // In addition to the best approximation kept by the CR base class, we keep 1560 // the entire sequence b[n], to the extent we've needed it so far. Each 1561 // reevaluation leads to slightly different sqrt arguments, but the 1562 // previous result can be used to avoid repeating low precision Newton 1563 // iterations for the sqrt approximation. 1564 ArrayList<Integer> b_prec = new ArrayList<Integer>(); 1565 ArrayList<BigInteger> b_val = new ArrayList<BigInteger>(); gl_pi_CR()1566 gl_pi_CR() { 1567 b_prec.add(null); // Zeroth entry unused. 1568 b_val.add(null); 1569 } 1570 private static BigInteger TOLERANCE = BigInteger.valueOf(4); 1571 // sqrt(1/2) 1572 private static CR SQRT_HALF = new sqrt_CR(ONE.shiftRight(1)); 1573 approximate(int p)1574 protected BigInteger approximate(int p) { 1575 // Rough approximations are easy. 1576 if (p >= 0) return scale(BigInteger.valueOf(3), -p); 1577 // We need roughly log2(p) iterations. Each iteration should 1578 // contribute no more than 2 ulps to the error in the corresponding 1579 // term (a[n], b[n], or t[n]). Thus 2log2(n) bits plus a few for the 1580 // final calulation and rounding suffice. 1581 final int extra_eval_prec = 1582 (int)Math.ceil(Math.log(-p) / Math.log(2)) + 10; 1583 // All our terms are implicitly scaled by eval_prec. 1584 final int eval_prec = p - extra_eval_prec; 1585 BigInteger a = BigInteger.ONE.shiftLeft(-eval_prec); 1586 BigInteger b = SQRT_HALF.get_appr(eval_prec); 1587 BigInteger t = BigInteger.ONE.shiftLeft(-eval_prec - 2); 1588 int n = 0; 1589 while (a.subtract(b).subtract(TOLERANCE).signum() > 0) { 1590 // Current values correspond to n, next_ values to n + 1 1591 // b_prec.size() == b_val.size() >= n + 1 1592 final BigInteger next_a = a.add(b).shiftRight(1); 1593 final BigInteger a_diff = a.subtract(next_a); 1594 CR next_b_as_CR; 1595 final BigInteger b_prod = a.multiply(b).shiftRight(-eval_prec); 1596 // We the compute square root approximations using a nested 1597 // temporary CR computation, to avoid implementing BigInteger 1598 // square roots separately. 1599 final CR b_prod_as_CR = CR.valueOf(b_prod).shiftRight(-eval_prec); 1600 if (b_prec.size() == n + 1) { 1601 // Need an n+1st slot. 1602 b_prec.add(null); 1603 b_val.add(null); 1604 next_b_as_CR = b_prod_as_CR.sqrt(); 1605 } else { 1606 // Reuse previous approximation to reduce sqrt iterations, 1607 // hopefully to one. 1608 next_b_as_CR = new sqrt_CR(b_prod_as_CR, b_prec.get(n + 1), 1609 b_val.get(n + 1)); 1610 } 1611 // b_prec.size() == b_val.size() >= n + 2 1612 final BigInteger next_b = next_b_as_CR.get_appr(eval_prec); 1613 b_prec.set(n + 1, Integer.valueOf(p)); 1614 b_val.set(n + 1, scale(next_b, -extra_eval_prec)); 1615 final BigInteger next_t = 1616 t.subtract(a_diff.multiply(a_diff) 1617 .shiftLeft(n + eval_prec)); // shift dist. usually neg. 1618 a = next_a; 1619 b = next_b; 1620 t = next_t; 1621 ++n; 1622 } 1623 final BigInteger sum = a.add(b); 1624 final BigInteger result = sum.multiply(sum).divide(t).shiftRight(2); 1625 return scale(result, -extra_eval_prec); 1626 } 1627 } 1628