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