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