/* * Copyright (C) 2015 The Android Open Source Project * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. */ /* * The above license covers additions and changes by AOSP authors. * The original code is licensed as follows: */ // // Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED // // Permission is granted free of charge to copy, modify, use and distribute // this software provided you include the entirety of this notice in all // copies made. // // THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY // KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, // WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT // FOR A PARTICULAR PURPOSE OR NON-INFRINGING. SGI ASSUMES NO RISK AS TO THE // QUALITY AND PERFORMANCE OF THE SOFTWARE. SHOULD THE SOFTWARE PROVE // DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES // AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS // AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. // // UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, // WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR // OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, // INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE // SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK // STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL // OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF // THE POSSIBILITY OF SUCH DAMAGES. THIS LIMITATION OF LIABILITY SHALL NOT // APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE // LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE // EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT // EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. // // These license terms shall be governed by and construed in accordance with // the laws of the United States and the State of California as applied to // agreements entered into and to be performed entirely within California // between California residents. Any litigation relating to these license // terms shall be subject to the exclusive jurisdiction of the Federal Courts // of the Northern District of California (or, absent subject matter // jurisdiction in such courts, the courts of the State of California), with // venue lying exclusively in Santa Clara County, California. // Copyright (c) 2001-2004, Hewlett-Packard Development Company, L.P. // // Permission is granted free of charge to copy, modify, use and distribute // this software provided you include the entirety of this notice in all // copies made. // // THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY // KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, // WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT // FOR A PARTICULAR PURPOSE OR NON-INFRINGING. HEWLETT-PACKARD ASSUMES // NO RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE. // SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, // HEWLETT-PACKARD ASSUMES NO COST OR LIABILITY FOR ANY // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES // AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS // AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. // // UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, // WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR // OTHERWISE, SHALL HEWLETT-PACKARD BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, // INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE // SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK // STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL // OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF HEWLETT-PACKARD SHALL // HAVE BEEN INFORMED OF THE POSSIBILITY OF SUCH DAMAGES. // THIS LIMITATION OF LIABILITY SHALL NOT APPLY TO LIABILITY RESULTING // FROM HEWLETT-PACKARD's NEGLIGENCE TO THE EXTENT APPLICABLE // LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE // EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT // EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. // // Added valueOf(string, radix), fixed some documentation comments. // Hans_Boehm@hp.com 1/12/2001 // Fixed a serious typo in inv_CR(): For negative arguments it produced // the wrong sign. This affected the sign of divisions. // Added byteValue and fixed some comments. Hans.Boehm@hp.com 12/17/2002 // Added toStringFloatRep. Hans.Boehm@hp.com 4/1/2004 // Added get_appr() synchronization to allow access from multiple threads // hboehm@google.com 4/25/2014 // Changed cos() prescaling to avoid logarithmic depth tree. // hboehm@google.com 6/30/2014 // Added explicit asin() implementation. Remove one. Add ZERO and ONE and // make them public. hboehm@google.com 5/21/2015 // Added Gauss-Legendre PI implementation. Removed two. // hboehm@google.com 4/12/2016 // Fix shift operation in doubleValue. That produced incorrect values for // large negative exponents. // Don't negate argument and compute inverse for exp(). That causes severe // performance problems for (-huge).exp() // hboehm@google.com 8/21/2017 // Have comparison check for interruption. hboehm@google.com 10/31/2017 // Fix precision overflow issue in most general compareTo function. // Fix a couple of unused variable bugs. Notably selector_sign was // accidentally locally redeclared. (This turns out to be safe but useless.) // hboehm@google.com 11/20/2018. // Fix an exception-safety issue in gl_pi_CR.approximate. // hboehm@google.com 3/3/2019. package com.hp.creals; import java.math.BigInteger; import java.util.ArrayList; /** * Constructive real numbers, also known as recursive, or computable reals. * Each recursive real number is represented as an object that provides an * approximation function for the real number. * The approximation function guarantees that the generated approximation * is accurate to the specified precision. * Arithmetic operations on constructive reals produce new such objects; * they typically do not perform any real computation. * In this sense, arithmetic computations are exact: They produce * a description which describes the exact answer, and can be used to * later approximate it to arbitrary precision. *

* When approximations are generated, e.g. for output, they are * accurate to the requested precision; no cumulative rounding errors * are visible. * In order to achieve this precision, the approximation function will often * need to approximate subexpressions to greater precision than was originally * demanded. Thus the approximation of a constructive real number * generated through a complex sequence of operations may eventually require * evaluation to very high precision. This usually makes such computations * prohibitively expensive for large numerical problems. * But it is perfectly appropriate for use in a desk calculator, * for small numerical problems, for the evaluation of expressions * computated by a symbolic algebra system, for testing of accuracy claims * for floating point code on small inputs, or the like. *

* We expect that the vast majority of uses will ignore the particular * implementation, and the member functons approximate * and get_appr. Such applications will treat CR as * a conventional numerical type, with an interface modelled on * java.math.BigInteger. No subclasses of CR * will be explicitly mentioned by such a program. *

* All standard arithmetic operations, as well as a few algebraic * and transcendal functions are provided. Constructive reals are * immutable; thus all of these operations return a new constructive real. *

* A few uses will require explicit construction of approximation functions. * The requires the construction of a subclass of CR with * an overridden approximate function. Note that approximate * should only be defined, but never called. get_appr * provides the same functionality, but adds the caching necessary to obtain * reasonable performance. *

* Any operation may throw com.hp.creals.AbortedException if the thread * in which it is executing is interrupted. (InterruptedException * cannot be used for this purpose, since CR inherits from Number.) *

* Any operation may also throw com.hp.creals.PrecisionOverflowException * If the precision request generated during any subcalculation overflows * a 28-bit integer. (This should be extremely unlikely, except as an * outcome of a division by zero, or other erroneous computation.) * */ public abstract class CR extends Number { // CR is the basic representation of a number. // Abstractly this is a function for computing an approximation // plus the current best approximation. // We could do without the latter, but that would // be atrociously slow. /** * Indicates a constructive real operation was interrupted. * Most constructive real operations may throw such an exception. * This is unchecked, since Number methods may not raise checked * exceptions. */ public static class AbortedException extends RuntimeException { public AbortedException() { super(); } public AbortedException(String s) { super(s); } } /** * Indicates that the number of bits of precision requested by * a computation on constructive reals required more than 28 bits, * and was thus in danger of overflowing an int. * This is likely to be a symptom of a diverging computation, * e.g. division by zero. */ public static class PrecisionOverflowException extends RuntimeException { public PrecisionOverflowException() { super(); } public PrecisionOverflowException(String s) { super(s); } } // First some frequently used constants, so we don't have to // recompute these all over the place. static final BigInteger big0 = BigInteger.ZERO; static final BigInteger big1 = BigInteger.ONE; static final BigInteger bigm1 = BigInteger.valueOf(-1); static final BigInteger big2 = BigInteger.valueOf(2); static final BigInteger bigm2 = BigInteger.valueOf(-2); static final BigInteger big3 = BigInteger.valueOf(3); static final BigInteger big6 = BigInteger.valueOf(6); static final BigInteger big8 = BigInteger.valueOf(8); static final BigInteger big10 = BigInteger.TEN; static final BigInteger big750 = BigInteger.valueOf(750); static final BigInteger bigm750 = BigInteger.valueOf(-750); /** * Setting this to true requests that all computations be aborted by * throwing AbortedException. Must be rest to false before any further * computation. Ideally Thread.interrupt() should be used instead, but * that doesn't appear to be consistently supported by browser VMs. */ public volatile static boolean please_stop = false; /** * Must be defined in subclasses of CR. * Most users can ignore the existence of this method, and will * not ever need to define a CR subclass. * Returns value / 2 ** precision rounded to an integer. * The error in the result is strictly < 1. * Informally, approximate(n) gives a scaled approximation * accurate to 2**n. * Implementations may safely assume that precision is * at least a factor of 8 away from overflow. * Called only with the lock on the CR object * already held. */ protected abstract BigInteger approximate(int precision); transient int min_prec; // The smallest precision value with which the above // has been called. transient BigInteger max_appr; // The scaled approximation corresponding to min_prec. transient boolean appr_valid = false; // min_prec and max_val are valid. // Helper functions static int bound_log2(int n) { int abs_n = Math.abs(n); return (int)Math.ceil(Math.log((double)(abs_n + 1))/Math.log(2.0)); } // Check that a precision is at least a factor of 8 away from // overflowng the integer used to hold a precision spec. // We generally perform this check early on, and then convince // ourselves that none of the operations performed on precisions // inside a function can generate an overflow. static void check_prec(int n) { int high = n >> 28; // if n is not in danger of overflowing, then the 4 high order // bits should be identical. Thus high is either 0 or -1. // The rest of this is to test for either of those in a way // that should be as cheap as possible. int high_shifted = n >> 29; if (0 != (high ^ high_shifted)) { throw new PrecisionOverflowException(); } } /** * The constructive real number corresponding to a * BigInteger. */ public static CR valueOf(BigInteger n) { return new int_CR(n); } /** * The constructive real number corresponding to a * Java int. */ public static CR valueOf(int n) { return valueOf(BigInteger.valueOf(n)); } /** * The constructive real number corresponding to a * Java long. */ public static CR valueOf(long n) { return valueOf(BigInteger.valueOf(n)); } /** * The constructive real number corresponding to a * Java double. * The result is undefined if argument is infinite or NaN. */ public static CR valueOf(double n) { if (Double.isNaN(n)) throw new ArithmeticException("Nan argument"); if (Double.isInfinite(n)) { throw new ArithmeticException("Infinite argument"); } boolean negative = (n < 0.0); long bits = Double.doubleToLongBits(Math.abs(n)); long mantissa = (bits & 0xfffffffffffffL); int biased_exp = (int)(bits >> 52); int exp = biased_exp - 1075; if (biased_exp != 0) { mantissa += (1L << 52); } else { mantissa <<= 1; } CR result = valueOf(mantissa).shiftLeft(exp); if (negative) result = result.negate(); return result; } /** * The constructive real number corresponding to a * Java float. * The result is undefined if argument is infinite or NaN. */ public static CR valueOf(float n) { return valueOf((double) n); } public static CR ZERO = valueOf(0); public static CR ONE = valueOf(1); // Multiply k by 2**n. static BigInteger shift(BigInteger k, int n) { if (n == 0) return k; if (n < 0) return k.shiftRight(-n); return k.shiftLeft(n); } // Multiply by 2**n, rounding result static BigInteger scale(BigInteger k, int n) { if (n >= 0) { return k.shiftLeft(n); } else { BigInteger adj_k = shift(k, n+1).add(big1); return adj_k.shiftRight(1); } } // Identical to approximate(), but maintain and update cache. /** * Returns value / 2 ** prec rounded to an integer. * The error in the result is strictly < 1. * Produces the same answer as approximate, but uses and * maintains a cached approximation. * Normally not overridden, and called only from approximate * methods in subclasses. Not needed if the provided operations * on constructive reals suffice. */ public synchronized BigInteger get_appr(int precision) { check_prec(precision); if (appr_valid && precision >= min_prec) { return scale(max_appr, min_prec - precision); } else { BigInteger result = approximate(precision); min_prec = precision; max_appr = result; appr_valid = true; return result; } } // Return the position of the msd. // If x.msd() == n then // 2**(n-1) < abs(x) < 2**(n+1) // This initial version assumes that max_appr is valid // and sufficiently removed from zero // that the msd is determined. int known_msd() { int first_digit; int length; if (max_appr.signum() >= 0) { length = max_appr.bitLength(); } else { length = max_appr.negate().bitLength(); } first_digit = min_prec + length - 1; return first_digit; } // This version may return Integer.MIN_VALUE if the correct // answer is < n. int msd(int n) { if (!appr_valid || max_appr.compareTo(big1) <= 0 && max_appr.compareTo(bigm1) >= 0) { get_appr(n - 1); if (max_appr.abs().compareTo(big1) <= 0) { // msd could still be arbitrarily far to the right. return Integer.MIN_VALUE; } } return known_msd(); } // Functionally equivalent, but iteratively evaluates to higher // precision. int iter_msd(int n) { int prec = 0; for (;prec > n + 30; prec = (prec * 3)/2 - 16) { int msd = msd(prec); if (msd != Integer.MIN_VALUE) return msd; check_prec(prec); if (Thread.interrupted() || please_stop) { throw new AbortedException(); } } return msd(n); } // This version returns a correct answer eventually, except // that it loops forever (or throws an exception when the // requested precision overflows) if this constructive real is zero. int msd() { return iter_msd(Integer.MIN_VALUE); } // A helper function for toString. // Generate a String containing n zeroes. private static String zeroes(int n) { char[] a = new char[n]; for (int i = 0; i < n; ++i) { a[i] = '0'; } return new String(a); } // Natural log of 2. Needed for some prescaling below. // ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) CR simple_ln() { return new prescaled_ln_CR(this.subtract(ONE)); } static CR ten_ninths = valueOf(10).divide(valueOf(9)); static CR twentyfive_twentyfourths = valueOf(25).divide(valueOf(24)); static CR eightyone_eightyeths = valueOf(81).divide(valueOf(80)); static CR ln2_1 = valueOf(7).multiply(ten_ninths.simple_ln()); static CR ln2_2 = valueOf(2).multiply(twentyfive_twentyfourths.simple_ln()); static CR ln2_3 = valueOf(3).multiply(eightyone_eightyeths.simple_ln()); static CR ln2 = ln2_1.subtract(ln2_2).add(ln2_3); // Atan of integer reciprocal. Used for atan_PI. Could perhaps be made // public. static CR atan_reciprocal(int n) { return new integral_atan_CR(n); } // Other constants used for PI computation. static CR four = valueOf(4); // Public operations. /** * Return 0 if x = y to within the indicated tolerance, * -1 if x < y, and +1 if x > y. If x and y are indeed * equal, it is guaranteed that 0 will be returned. If * they differ by less than the tolerance, anything * may happen. The tolerance allowed is * the maximum of (abs(this)+abs(x))*(2**r) and 2**a * @param x The other constructive real * @param r Relative tolerance in bits * @param a Absolute tolerance in bits */ public int compareTo(CR x, int r, int a) { int this_msd = iter_msd(a); int x_msd = x.iter_msd(this_msd > a? this_msd : a); int max_msd = (x_msd > this_msd? x_msd : this_msd); if (max_msd == Integer.MIN_VALUE) { return 0; } check_prec(r); int rel = max_msd + r; int abs_prec = (rel > a? rel : a); return compareTo(x, abs_prec); } /** * Approximate comparison with only an absolute tolerance. * Identical to the three argument version, but without a relative * tolerance. * Result is 0 if both constructive reals are equal, indeterminate * if they differ by less than 2**a. * * @param x The other constructive real * @param a Absolute tolerance in bits */ public int compareTo(CR x, int a) { int needed_prec = a - 1; BigInteger this_appr = get_appr(needed_prec); BigInteger x_appr = x.get_appr(needed_prec); int comp1 = this_appr.compareTo(x_appr.add(big1)); if (comp1 > 0) return 1; int comp2 = this_appr.compareTo(x_appr.subtract(big1)); if (comp2 < 0) return -1; return 0; } /** * Return -1 if this < x, or +1 if this > x. * Should be called only if this != x. * If this == x, this will not terminate correctly; typically it * will run until it exhausts memory. * If the two constructive reals may be equal, the two or 3 argument * version of compareTo should be used. */ public int compareTo(CR x) { for (int a = -20; ; a *= 2) { check_prec(a); int result = compareTo(x, a); if (0 != result) return result; if (Thread.interrupted() || please_stop) { throw new AbortedException(); } } } /** * Equivalent to compareTo(CR.valueOf(0), a) */ public int signum(int a) { if (appr_valid) { int quick_try = max_appr.signum(); if (0 != quick_try) return quick_try; } int needed_prec = a - 1; BigInteger this_appr = get_appr(needed_prec); return this_appr.signum(); } /** * Return -1 if negative, +1 if positive. * Should be called only if this != 0. * In the 0 case, this will not terminate correctly; typically it * will run until it exhausts memory. * If the two constructive reals may be equal, the one or two argument * version of signum should be used. */ public int signum() { for (int a = -20; ; a *= 2) { check_prec(a); int result = signum(a); if (0 != result) return result; if (Thread.interrupted() || please_stop) { throw new AbortedException(); } } } /** * Return the constructive real number corresponding to the given * textual representation and radix. * * @param s [-] digit* [. digit*] * @param radix */ public static CR valueOf(String s, int radix) throws NumberFormatException { int len = s.length(); int start_pos = 0, point_pos; String fraction; while (s.charAt(start_pos) == ' ') ++start_pos; while (s.charAt(len - 1) == ' ') --len; point_pos = s.indexOf('.', start_pos); if (point_pos == -1) { point_pos = len; fraction = "0"; } else { fraction = s.substring(point_pos + 1, len); } String whole = s.substring(start_pos, point_pos); BigInteger scaled_result = new BigInteger(whole + fraction, radix); BigInteger divisor = BigInteger.valueOf(radix).pow(fraction.length()); return CR.valueOf(scaled_result).divide(CR.valueOf(divisor)); } /** * Return a textual representation accurate to n places * to the right of the decimal point. n must be nonnegative. * * @param n Number of digits (>= 0) included to the right of decimal point * @param radix Base ( >= 2, <= 16) for the resulting representation. */ public String toString(int n, int radix) { CR scaled_CR; if (16 == radix) { scaled_CR = shiftLeft(4*n); } else { BigInteger scale_factor = BigInteger.valueOf(radix).pow(n); scaled_CR = multiply(new int_CR(scale_factor)); } BigInteger scaled_int = scaled_CR.get_appr(0); String scaled_string = scaled_int.abs().toString(radix); String result; if (0 == n) { result = scaled_string; } else { int len = scaled_string.length(); if (len <= n) { // Add sufficient leading zeroes String z = zeroes(n + 1 - len); scaled_string = z + scaled_string; len = n + 1; } String whole = scaled_string.substring(0, len - n); String fraction = scaled_string.substring(len - n); result = whole + "." + fraction; } if (scaled_int.signum() < 0) { result = "-" + result; } return result; } /** * Equivalent to toString(n,10) * * @param n Number of digits included to the right of decimal point */ public String toString(int n) { return toString(n, 10); } /** * Equivalent to toString(10, 10) */ public String toString() { return toString(10); } static double doubleLog2 = Math.log(2.0); /** * Return a textual scientific notation representation accurate * to n places to the right of the decimal point. * n must be nonnegative. A value smaller than * radix**-m may be displayed as 0. * The mantissa component of the result is either "0" * or exactly n digits long. The sign * component is zero exactly when the mantissa is "0". * * @param n Number of digits (> 0) included to the right of decimal point. * @param radix Base ( ≥ 2, ≤ 16) for the resulting representation. * @param m Precision used to distinguish number from zero. * Expressed as a power of m. */ public StringFloatRep toStringFloatRep(int n, int radix, int m) { if (n <= 0) throw new ArithmeticException("Bad precision argument"); double log2_radix = Math.log((double)radix)/doubleLog2; BigInteger big_radix = BigInteger.valueOf(radix); long long_msd_prec = (long)(log2_radix * (double)m); if (long_msd_prec > (long)Integer.MAX_VALUE || long_msd_prec < (long)Integer.MIN_VALUE) throw new PrecisionOverflowException(); int msd_prec = (int)long_msd_prec; check_prec(msd_prec); int msd = iter_msd(msd_prec - 2); if (msd == Integer.MIN_VALUE) return new StringFloatRep(0, "0", radix, 0); int exponent = (int)Math.ceil((double)msd / log2_radix); // Guess for the exponent. Try to get it usually right. int scale_exp = exponent - n; CR scale; if (scale_exp > 0) { scale = CR.valueOf(big_radix.pow(scale_exp)).inverse(); } else { scale = CR.valueOf(big_radix.pow(-scale_exp)); } CR scaled_res = multiply(scale); BigInteger scaled_int = scaled_res.get_appr(0); int sign = scaled_int.signum(); String scaled_string = scaled_int.abs().toString(radix); while (scaled_string.length() < n) { // exponent was too large. Adjust. scaled_res = scaled_res.multiply(CR.valueOf(big_radix)); exponent -= 1; scaled_int = scaled_res.get_appr(0); sign = scaled_int.signum(); scaled_string = scaled_int.abs().toString(radix); } if (scaled_string.length() > n) { // exponent was too small. Adjust by truncating. exponent += (scaled_string.length() - n); scaled_string = scaled_string.substring(0, n); } return new StringFloatRep(sign, scaled_string, radix, exponent); } /** * Return a BigInteger which differs by less than one from the * constructive real. */ public BigInteger BigIntegerValue() { return get_appr(0); } /** * Return an int which differs by less than one from the * constructive real. Behavior on overflow is undefined. */ public int intValue() { return BigIntegerValue().intValue(); } /** * Return an int which differs by less than one from the * constructive real. Behavior on overflow is undefined. */ public byte byteValue() { return BigIntegerValue().byteValue(); } /** * Return a long which differs by less than one from the * constructive real. Behavior on overflow is undefined. */ public long longValue() { return BigIntegerValue().longValue(); } /** * Return a double which differs by less than one in the least * represented bit from the constructive real. * (We're in fact closer to round-to-nearest than that, but we can't and * don't promise correct rounding.) */ public double doubleValue() { int my_msd = iter_msd(-1080 /* slightly > exp. range */); if (Integer.MIN_VALUE == my_msd) return 0.0; int needed_prec = my_msd - 60; double scaled_int = get_appr(needed_prec).doubleValue(); boolean may_underflow = (needed_prec < -1000); long scaled_int_rep = Double.doubleToLongBits(scaled_int); long exp_adj = may_underflow? needed_prec + 96 : needed_prec; long orig_exp = (scaled_int_rep >> 52) & 0x7ff; if (((orig_exp + exp_adj) & ~0x7ff) != 0) { // Original unbiased exponent is > 50. Exp_adj > -1050. // Thus this can overflow the 11 bit exponent only if the result // itself overflows. if (scaled_int < 0.0) { return Double.NEGATIVE_INFINITY; } else { return Double.POSITIVE_INFINITY; } } scaled_int_rep += exp_adj << 52; double result = Double.longBitsToDouble(scaled_int_rep); if (may_underflow) { double two48 = (double)(1L << 48); return result/two48/two48; } else { return result; } } /** * Return a float which differs by less than one in the least * represented bit from the constructive real. */ public float floatValue() { return (float)doubleValue(); // Note that double-rounding is not a problem here, since we // cannot, and do not, guarantee correct rounding. } /** * Add two constructive reals. */ public CR add(CR x) { return new add_CR(this, x); } /** * Multiply a constructive real by 2**n. * @param n shift count, may be negative */ public CR shiftLeft(int n) { check_prec(n); return new shifted_CR(this, n); } /** * Multiply a constructive real by 2**(-n). * @param n shift count, may be negative */ public CR shiftRight(int n) { check_prec(n); return new shifted_CR(this, -n); } /** * Produce a constructive real equivalent to the original, assuming * the original was an integer. Undefined results if the original * was not an integer. Prevents evaluation of digits to the right * of the decimal point, and may thus improve performance. */ public CR assumeInt() { return new assumed_int_CR(this); } /** * The additive inverse of a constructive real */ public CR negate() { return new neg_CR(this); } /** * The difference between two constructive reals */ public CR subtract(CR x) { return new add_CR(this, x.negate()); } /** * The product of two constructive reals */ public CR multiply(CR x) { return new mult_CR(this, x); } /** * The multiplicative inverse of a constructive real. * x.inverse() is equivalent to CR.valueOf(1).divide(x). */ public CR inverse() { return new inv_CR(this); } /** * The quotient of two constructive reals. */ public CR divide(CR x) { return new mult_CR(this, x.inverse()); } /** * The real number x if this < 0, or y otherwise. * Requires x = y if this = 0. * Since comparisons may diverge, this is often * a useful alternative to conditionals. */ public CR select(CR x, CR y) { return new select_CR(this, x, y); } /** * The maximum of two constructive reals. */ public CR max(CR x) { return subtract(x).select(x, this); } /** * The minimum of two constructive reals. */ public CR min(CR x) { return subtract(x).select(this, x); } /** * The absolute value of a constructive reals. * Note that this cannot be written as a conditional. */ public CR abs() { return select(negate(), this); } /** * The exponential function, that is e**this. */ public CR exp() { final int low_prec = -10; BigInteger rough_appr = get_appr(low_prec); // Handle negative arguments directly; negating and computing inverse // can be very expensive. if (rough_appr.compareTo(big2) > 0 || rough_appr.compareTo(bigm2) < 0) { CR square_root = shiftRight(1).exp(); return square_root.multiply(square_root); } else { return new prescaled_exp_CR(this); } } /** * The ratio of a circle's circumference to its diameter. */ public static CR PI = new gl_pi_CR(); // Our old PI implementation. Keep this around for now to allow checking. // This implementation may also be faster for BigInteger implementations // that support only quadratic multiplication, but exhibit high performance // for small computations. (The standard Android 6 implementation supports // subquadratic multiplication, but has high constant overhead.) Many other // atan-based formulas are possible, but based on superficial // experimentation, this is roughly as good as the more complex formulas. public static CR atan_PI = four.multiply(four.multiply(atan_reciprocal(5)) .subtract(atan_reciprocal(239))); // pi/4 = 4*atan(1/5) - atan(1/239) static CR half_pi = PI.shiftRight(1); /** * The trigonometric cosine function. */ public CR cos() { BigInteger halfpi_multiples = divide(PI).get_appr(-1); BigInteger abs_halfpi_multiples = halfpi_multiples.abs(); if (abs_halfpi_multiples.compareTo(big2) >= 0) { // Subtract multiples of PI BigInteger pi_multiples = scale(halfpi_multiples, -1); CR adjustment = PI.multiply(CR.valueOf(pi_multiples)); if (pi_multiples.and(big1).signum() != 0) { return subtract(adjustment).cos().negate(); } else { return subtract(adjustment).cos(); } } else if (get_appr(-1).abs().compareTo(big2) >= 0) { // Scale further with double angle formula CR cos_half = shiftRight(1).cos(); return cos_half.multiply(cos_half).shiftLeft(1).subtract(ONE); } else { return new prescaled_cos_CR(this); } } /** * The trigonometric sine function. */ public CR sin() { return half_pi.subtract(this).cos(); } /** * The trignonometric arc (inverse) sine function. */ public CR asin() { BigInteger rough_appr = get_appr(-10); if (rough_appr.compareTo(big750) /* 1/sqrt(2) + a bit */ > 0){ CR new_arg = ONE.subtract(multiply(this)).sqrt(); return new_arg.acos(); } else if (rough_appr.compareTo(bigm750) < 0) { return negate().asin().negate(); } else { return new prescaled_asin_CR(this); } } /** * The trignonometric arc (inverse) cosine function. */ public CR acos() { return half_pi.subtract(asin()); } static final BigInteger low_ln_limit = big8; /* sixteenths, i.e. 1/2 */ static final BigInteger high_ln_limit = BigInteger.valueOf(16 + 8 /* 1.5 */); static final BigInteger scaled_4 = BigInteger.valueOf(4*16); /** * The natural (base e) logarithm. */ public CR ln() { final int low_prec = -4; BigInteger rough_appr = get_appr(low_prec); /* In sixteenths */ if (rough_appr.compareTo(big0) < 0) { throw new ArithmeticException("ln(negative)"); } if (rough_appr.compareTo(low_ln_limit) <= 0) { return inverse().ln().negate(); } if (rough_appr.compareTo(high_ln_limit) >= 0) { if (rough_appr.compareTo(scaled_4) <= 0) { CR quarter = sqrt().sqrt().ln(); return quarter.shiftLeft(2); } else { int extra_bits = rough_appr.bitLength() - 3; CR scaled_result = shiftRight(extra_bits).ln(); return scaled_result.add(CR.valueOf(extra_bits).multiply(ln2)); } } return simple_ln(); } /** * The square root of a constructive real. */ public CR sqrt() { return new sqrt_CR(this); } } // end of CR // // A specialization of CR for cases in which approximate() calls // to increase evaluation precision are somewhat expensive. // If we need to (re)evaluate, we speculatively evaluate to slightly // higher precision, miminimizing reevaluations. // Note that this requires any arguments to be evaluated to higher // precision than absolutely necessary. It can thus potentially // result in lots of wasted effort, and should be used judiciously. // This assumes that the order of magnitude of the number is roughly one. // abstract class slow_CR extends CR { static int max_prec = -64; static int prec_incr = 32; public synchronized BigInteger get_appr(int precision) { check_prec(precision); if (appr_valid && precision >= min_prec) { return scale(max_appr, min_prec - precision); } else { int eval_prec = (precision >= max_prec? max_prec : (precision - prec_incr + 1) & ~(prec_incr - 1)); BigInteger result = approximate(eval_prec); min_prec = eval_prec; max_appr = result; appr_valid = true; return scale(result, eval_prec - precision); } } } // Representation of an integer constant. Private. class int_CR extends CR { BigInteger value; int_CR(BigInteger n) { value = n; } protected BigInteger approximate(int p) { return scale(value, -p) ; } } // Representation of a number that may not have been completely // evaluated, but is assumed to be an integer. Hence we never // evaluate beyond the decimal point. class assumed_int_CR extends CR { CR value; assumed_int_CR(CR x) { value = x; } protected BigInteger approximate(int p) { if (p >= 0) { return value.get_appr(p); } else { return scale(value.get_appr(0), -p) ; } } } // Representation of the sum of 2 constructive reals. Private. class add_CR extends CR { CR op1; CR op2; add_CR(CR x, CR y) { op1 = x; op2 = y; } protected BigInteger approximate(int p) { // Args need to be evaluated so that each error is < 1/4 ulp. // Rounding error from the cale call is <= 1/2 ulp, so that // final error is < 1 ulp. return scale(op1.get_appr(p-2).add(op2.get_appr(p-2)), -2); } } // Representation of a CR multiplied by 2**n class shifted_CR extends CR { CR op; int count; shifted_CR(CR x, int n) { op = x; count = n; } protected BigInteger approximate(int p) { return op.get_appr(p - count); } } // Representation of the negation of a constructive real. Private. class neg_CR extends CR { CR op; neg_CR(CR x) { op = x; } protected BigInteger approximate(int p) { return op.get_appr(p).negate(); } } // Representation of: // op1 if selector < 0 // op2 if selector >= 0 // Assumes x = y if s = 0 class select_CR extends CR { CR selector; int selector_sign; CR op1; CR op2; select_CR(CR s, CR x, CR y) { selector = s; selector_sign = selector.get_appr(-20).signum(); op1 = x; op2 = y; } protected BigInteger approximate(int p) { if (selector_sign < 0) return op1.get_appr(p); if (selector_sign > 0) return op2.get_appr(p); BigInteger op1_appr = op1.get_appr(p-1); BigInteger op2_appr = op2.get_appr(p-1); BigInteger diff = op1_appr.subtract(op2_appr).abs(); if (diff.compareTo(big1) <= 0) { // close enough; use either return scale(op1_appr, -1); } // op1 and op2 are different; selector != 0; // safe to get sign of selector. if (selector.signum() < 0) { selector_sign = -1; return scale(op1_appr, -1); } else { selector_sign = 1; return scale(op2_appr, -1); } } } // Representation of the product of 2 constructive reals. Private. class mult_CR extends CR { CR op1; CR op2; mult_CR(CR x, CR y) { op1 = x; op2 = y; } protected BigInteger approximate(int p) { int half_prec = (p >> 1) - 1; int msd_op1 = op1.msd(half_prec); int msd_op2; if (msd_op1 == Integer.MIN_VALUE) { msd_op2 = op2.msd(half_prec); if (msd_op2 == Integer.MIN_VALUE) { // Product is small enough that zero will do as an // approximation. return big0; } else { // Swap them, so the larger operand (in absolute value) // is first. CR tmp; tmp = op1; op1 = op2; op2 = tmp; msd_op1 = msd_op2; } } // msd_op1 is valid at this point. int prec2 = p - msd_op1 - 3; // Precision needed for op2. // The appr. error is multiplied by at most // 2 ** (msd_op1 + 1) // Thus each approximation contributes 1/4 ulp // to the rounding error, and the final rounding adds // another 1/2 ulp. BigInteger appr2 = op2.get_appr(prec2); if (appr2.signum() == 0) return big0; msd_op2 = op2.known_msd(); int prec1 = p - msd_op2 - 3; // Precision needed for op1. BigInteger appr1 = op1.get_appr(prec1); int scale_digits = prec1 + prec2 - p; return scale(appr1.multiply(appr2), scale_digits); } } // Representation of the multiplicative inverse of a constructive // real. Private. Should use Newton iteration to refine estimates. class inv_CR extends CR { CR op; inv_CR(CR x) { op = x; } protected BigInteger approximate(int p) { int msd = op.msd(); int inv_msd = 1 - msd; int digits_needed = inv_msd - p + 3; // Number of SIGNIFICANT digits needed for // argument, excl. msd position, which may // be fictitious, since msd routine can be // off by 1. Roughly 1 extra digit is // needed since the relative error is the // same in the argument and result, but // this isn't quite the same as the number // of significant digits. Another digit // is needed to compensate for slop in the // calculation. // One further bit is required, since the // final rounding introduces a 0.5 ulp // error. int prec_needed = msd - digits_needed; int log_scale_factor = -p - prec_needed; if (log_scale_factor < 0) return big0; BigInteger dividend = big1.shiftLeft(log_scale_factor); BigInteger scaled_divisor = op.get_appr(prec_needed); BigInteger abs_scaled_divisor = scaled_divisor.abs(); BigInteger adj_dividend = dividend.add( abs_scaled_divisor.shiftRight(1)); // Adjustment so that final result is rounded. BigInteger result = adj_dividend.divide(abs_scaled_divisor); if (scaled_divisor.signum() < 0) { return result.negate(); } else { return result; } } } // Representation of the exponential of a constructive real. Private. // Uses a Taylor series expansion. Assumes |x| < 1/2. // Note: this is known to be a bad algorithm for // floating point. Unfortunately, other alternatives // appear to require precomputed information. class prescaled_exp_CR extends CR { CR op; prescaled_exp_CR(CR x) { op = x; } protected BigInteger approximate(int p) { if (p >= 1) return big0; int iterations_needed = -p/2 + 2; // conservative estimate > 0. // Claim: each intermediate term is accurate // to 2*2^calc_precision. // Total rounding error in series computation is // 2*iterations_needed*2^calc_precision, // exclusive of error in op. int calc_precision = p - bound_log2(2*iterations_needed) - 4; // for error in op, truncation. int op_prec = p - 3; BigInteger op_appr = op.get_appr(op_prec); // Error in argument results in error of < 3/8 ulp. // Sum of term eval. rounding error is < 1/16 ulp. // Series truncation error < 1/16 ulp. // Final rounding error is <= 1/2 ulp. // Thus final error is < 1 ulp. BigInteger scaled_1 = big1.shiftLeft(-calc_precision); BigInteger current_term = scaled_1; BigInteger current_sum = scaled_1; int n = 0; BigInteger max_trunc_error = big1.shiftLeft(p - 4 - calc_precision); while (current_term.abs().compareTo(max_trunc_error) >= 0) { if (Thread.interrupted() || please_stop) throw new AbortedException(); n += 1; /* current_term = current_term * op / n */ current_term = scale(current_term.multiply(op_appr), op_prec); current_term = current_term.divide(BigInteger.valueOf(n)); current_sum = current_sum.add(current_term); } return scale(current_sum, calc_precision - p); } } // Representation of the cosine of a constructive real. Private. // Uses a Taylor series expansion. Assumes |x| < 1. class prescaled_cos_CR extends slow_CR { CR op; prescaled_cos_CR(CR x) { op = x; } protected BigInteger approximate(int p) { if (p >= 1) return big0; int iterations_needed = -p/2 + 4; // conservative estimate > 0. // Claim: each intermediate term is accurate // to 2*2^calc_precision. // Total rounding error in series computation is // 2*iterations_needed*2^calc_precision, // exclusive of error in op. int calc_precision = p - bound_log2(2*iterations_needed) - 4; // for error in op, truncation. int op_prec = p - 2; BigInteger op_appr = op.get_appr(op_prec); // Error in argument results in error of < 1/4 ulp. // Cumulative arithmetic rounding error is < 1/16 ulp. // Series truncation error < 1/16 ulp. // Final rounding error is <= 1/2 ulp. // Thus final error is < 1 ulp. BigInteger current_term; int n; BigInteger max_trunc_error = big1.shiftLeft(p - 4 - calc_precision); n = 0; current_term = big1.shiftLeft(-calc_precision); BigInteger current_sum = current_term; while (current_term.abs().compareTo(max_trunc_error) >= 0) { if (Thread.interrupted() || please_stop) throw new AbortedException(); n += 2; /* current_term = - current_term * op * op / n * (n - 1) */ current_term = scale(current_term.multiply(op_appr), op_prec); current_term = scale(current_term.multiply(op_appr), op_prec); BigInteger divisor = BigInteger.valueOf(-n) .multiply(BigInteger.valueOf(n-1)); current_term = current_term.divide(divisor); current_sum = current_sum.add(current_term); } return scale(current_sum, calc_precision - p); } } // The constructive real atan(1/n), where n is a small integer // > base. // This gives a simple and moderately fast way to compute PI. class integral_atan_CR extends slow_CR { int op; integral_atan_CR(int x) { op = x; } protected BigInteger approximate(int p) { if (p >= 1) return big0; int iterations_needed = -p/2 + 2; // conservative estimate > 0. // Claim: each intermediate term is accurate // to 2*base^calc_precision. // Total rounding error in series computation is // 2*iterations_needed*base^calc_precision, // exclusive of error in op. int calc_precision = p - bound_log2(2*iterations_needed) - 2; // for error in op, truncation. // Error in argument results in error of < 3/8 ulp. // Cumulative arithmetic rounding error is < 1/4 ulp. // Series truncation error < 1/4 ulp. // Final rounding error is <= 1/2 ulp. // Thus final error is < 1 ulp. BigInteger scaled_1 = big1.shiftLeft(-calc_precision); BigInteger big_op = BigInteger.valueOf(op); BigInteger big_op_squared = BigInteger.valueOf(op*op); BigInteger op_inverse = scaled_1.divide(big_op); BigInteger current_power = op_inverse; BigInteger current_term = op_inverse; BigInteger current_sum = op_inverse; int current_sign = 1; int n = 1; BigInteger max_trunc_error = big1.shiftLeft(p - 2 - calc_precision); while (current_term.abs().compareTo(max_trunc_error) >= 0) { if (Thread.interrupted() || please_stop) throw new AbortedException(); n += 2; current_power = current_power.divide(big_op_squared); current_sign = -current_sign; current_term = current_power.divide(BigInteger.valueOf(current_sign*n)); current_sum = current_sum.add(current_term); } return scale(current_sum, calc_precision - p); } } // Representation for ln(1 + op) class prescaled_ln_CR extends slow_CR { CR op; prescaled_ln_CR(CR x) { op = x; } // Compute an approximation of ln(1+x) to precision // prec. This assumes |x| < 1/2. // It uses a Taylor series expansion. // Unfortunately there appears to be no way to take // advantage of old information. // Note: this is known to be a bad algorithm for // floating point. Unfortunately, other alternatives // appear to require precomputed tabular information. protected BigInteger approximate(int p) { if (p >= 0) return big0; int iterations_needed = -p; // conservative estimate > 0. // Claim: each intermediate term is accurate // to 2*2^calc_precision. Total error is // 2*iterations_needed*2^calc_precision // exclusive of error in op. int calc_precision = p - bound_log2(2*iterations_needed) - 4; // for error in op, truncation. int op_prec = p - 3; BigInteger op_appr = op.get_appr(op_prec); // Error analysis as for exponential. BigInteger x_nth = scale(op_appr, op_prec - calc_precision); BigInteger current_term = x_nth; // x**n BigInteger current_sum = current_term; int n = 1; int current_sign = 1; // (-1)^(n-1) BigInteger max_trunc_error = big1.shiftLeft(p - 4 - calc_precision); while (current_term.abs().compareTo(max_trunc_error) >= 0) { if (Thread.interrupted() || please_stop) throw new AbortedException(); n += 1; current_sign = -current_sign; x_nth = scale(x_nth.multiply(op_appr), op_prec); current_term = x_nth.divide(BigInteger.valueOf(n * current_sign)); // x**n / (n * (-1)**(n-1)) current_sum = current_sum.add(current_term); } return scale(current_sum, calc_precision - p); } } // Representation of the arcsine of a constructive real. Private. // Uses a Taylor series expansion. Assumes |x| < (1/2)^(1/3). class prescaled_asin_CR extends slow_CR { CR op; prescaled_asin_CR(CR x) { op = x; } protected BigInteger approximate(int p) { // The Taylor series is the sum of x^(2n+1) * (2n)!/(4^n n!^2 (2n+1)) // Note that (2n)!/(4^n n!^2) is always less than one. // (The denominator is effectively 2n*2n*(2n-2)*(2n-2)*...*2*2 // which is clearly > (2n)!) // Thus all terms are bounded by x^(2n+1). // Unfortunately, there's no easy way to prescale the argument // to less than 1/sqrt(2), and we can only approximate that. // Thus the worst case iteration count is fairly high. // But it doesn't make much difference. if (p >= 2) return big0; // Never bigger than 4. int iterations_needed = -3 * p / 2 + 4; // conservative estimate > 0. // Follows from assumed bound on x and // the fact that only every other Taylor // Series term is present. // Claim: each intermediate term is accurate // to 2*2^calc_precision. // Total rounding error in series computation is // 2*iterations_needed*2^calc_precision, // exclusive of error in op. int calc_precision = p - bound_log2(2*iterations_needed) - 4; // for error in op, truncation. int op_prec = p - 3; // always <= -2 BigInteger op_appr = op.get_appr(op_prec); // Error in argument results in error of < 1/4 ulp. // (Derivative is bounded by 2 in the specified range and we use // 3 extra digits.) // Ignoring the argument error, each term has an error of // < 3ulps relative to calc_precision, which is more precise than p. // Cumulative arithmetic rounding error is < 3/16 ulp (relative to p). // Series truncation error < 2/16 ulp. (Each computed term // is at most 2/3 of last one, so some of remaining series < // 3/2 * current term.) // Final rounding error is <= 1/2 ulp. // Thus final error is < 1 ulp (relative to p). BigInteger max_last_term = big1.shiftLeft(p - 4 - calc_precision); int exp = 1; // Current exponent, = 2n+1 in above expression BigInteger current_term = op_appr.shiftLeft(op_prec - calc_precision); BigInteger current_sum = current_term; BigInteger current_factor = current_term; // Current scaled Taylor series term // before division by the exponent. // Accurate to 3 ulp at calc_precision. while (current_term.abs().compareTo(max_last_term) >= 0) { if (Thread.interrupted() || please_stop) throw new AbortedException(); exp += 2; // current_factor = current_factor * op * op * (exp-1) * (exp-2) / // (exp-1) * (exp-1), with the two exp-1 factors cancelling, // giving // current_factor = current_factor * op * op * (exp-2) / (exp-1) // Thus the error any in the previous term is multiplied by // op^2, adding an error of < (1/2)^(2/3) < 2/3 the original // error. current_factor = current_factor.multiply(BigInteger.valueOf(exp - 2)); current_factor = scale(current_factor.multiply(op_appr), op_prec + 2); // Carry 2 extra bits of precision forward; thus // this effectively introduces 1/8 ulp error. current_factor = current_factor.multiply(op_appr); BigInteger divisor = BigInteger.valueOf(exp - 1); current_factor = current_factor.divide(divisor); // Another 1/4 ulp error here. current_factor = scale(current_factor, op_prec - 2); // Remove extra 2 bits. 1/2 ulp rounding error. // Current_factor has original 3 ulp rounding error, which we // reduced by 1, plus < 1 ulp new rounding error. current_term = current_factor.divide(BigInteger.valueOf(exp)); // Contributes 1 ulp error to sum plus at most 3 ulp // from current_factor. current_sum = current_sum.add(current_term); } return scale(current_sum, calc_precision - p); } } class sqrt_CR extends CR { CR op; sqrt_CR(CR x) { op = x; } // Explicitly provide an initial approximation. // Useful for arithmetic geometric mean algorithms, where we've previously // computed a very similar square root. sqrt_CR(CR x, int min_p, BigInteger max_a) { op = x; min_prec = min_p; max_appr = max_a; appr_valid = true; } final int fp_prec = 50; // Conservative estimate of number of // significant bits in double precision // computation. final int fp_op_prec = 60; protected BigInteger approximate(int p) { int max_op_prec_needed = 2*p - 1; int msd = op.iter_msd(max_op_prec_needed); if (msd <= max_op_prec_needed) return big0; int result_msd = msd/2; // +- 1 int result_digits = result_msd - p; // +- 2 if (result_digits > fp_prec) { // Compute less precise approximation and use a Newton iter. int appr_digits = result_digits/2 + 6; // This should be conservative. Is fewer enough? int appr_prec = result_msd - appr_digits; int prod_prec = 2*appr_prec; // First compute the argument to maximal precision, so we don't end up // reevaluating it incrementally. BigInteger op_appr = op.get_appr(prod_prec); BigInteger last_appr = get_appr(appr_prec); // Compute (last_appr * last_appr + op_appr) / last_appr / 2 // while adjusting the scaling to make everything work BigInteger prod_prec_scaled_numerator = last_appr.multiply(last_appr).add(op_appr); BigInteger scaled_numerator = scale(prod_prec_scaled_numerator, appr_prec - p); BigInteger shifted_result = scaled_numerator.divide(last_appr); return shifted_result.add(big1).shiftRight(1); } else { // Use a double precision floating point approximation. // Make sure all precisions are even int op_prec = (msd - fp_op_prec) & ~1; int working_prec = op_prec - fp_op_prec; BigInteger scaled_bi_appr = op.get_appr(op_prec) .shiftLeft(fp_op_prec); double scaled_appr = scaled_bi_appr.doubleValue(); if (scaled_appr < 0.0) throw new ArithmeticException("sqrt(negative)"); double scaled_fp_sqrt = Math.sqrt(scaled_appr); BigInteger scaled_sqrt = BigInteger.valueOf((long)scaled_fp_sqrt); int shift_count = working_prec/2 - p; return shift(scaled_sqrt, shift_count); } } } // The constant PI, computed using the Gauss-Legendre alternating // arithmetic-geometric mean algorithm: // a[0] = 1 // b[0] = 1/sqrt(2) // t[0] = 1/4 // p[0] = 1 // // a[n+1] = (a[n] + b[n])/2 (arithmetic mean, between 0.8 and 1) // b[n+1] = sqrt(a[n] * b[n]) (geometric mean, between 0.7 and 1) // t[n+1] = t[n] - (2^n)(a[n]-a[n+1])^2, (always between 0.2 and 0.25) // // pi is then approximated as (a[n+1]+b[n+1])^2 / 4*t[n+1]. // class gl_pi_CR extends slow_CR { // In addition to the best approximation kept by the CR base class, we keep // the entire sequence b[n], to the extent we've needed it so far. Each // reevaluation leads to slightly different sqrt arguments, but the // previous result can be used to avoid repeating low precision Newton // iterations for the sqrt approximation. ArrayList b_prec = new ArrayList(); ArrayList b_val = new ArrayList(); gl_pi_CR() { b_prec.add(null); // Zeroth entry unused. b_val.add(null); } private static BigInteger TOLERANCE = BigInteger.valueOf(4); // sqrt(1/2) private static CR SQRT_HALF = new sqrt_CR(ONE.shiftRight(1)); protected BigInteger approximate(int p) { // Get us back into a consistent state if the last computation // was interrupted after pushing onto b_prec. if (b_prec.size() > b_val.size()) { b_prec.remove(b_prec.size() - 1); } // Rough approximations are easy. if (p >= 0) return scale(BigInteger.valueOf(3), -p); // We need roughly log2(p) iterations. Each iteration should // contribute no more than 2 ulps to the error in the corresponding // term (a[n], b[n], or t[n]). Thus 2log2(n) bits plus a few for the // final calulation and rounding suffice. final int extra_eval_prec = (int)Math.ceil(Math.log(-p) / Math.log(2)) + 10; // All our terms are implicitly scaled by eval_prec. final int eval_prec = p - extra_eval_prec; BigInteger a = BigInteger.ONE.shiftLeft(-eval_prec); BigInteger b = SQRT_HALF.get_appr(eval_prec); BigInteger t = BigInteger.ONE.shiftLeft(-eval_prec - 2); int n = 0; while (a.subtract(b).subtract(TOLERANCE).signum() > 0) { // Current values correspond to n, next_ values to n + 1 // b_prec.size() == b_val.size() >= n + 1 final BigInteger next_a = a.add(b).shiftRight(1); final BigInteger next_b; final BigInteger a_diff = a.subtract(next_a); final BigInteger b_prod = a.multiply(b).shiftRight(-eval_prec); // We compute square root approximations using a nested // temporary CR computation, to avoid implementing BigInteger // square roots separately. final CR b_prod_as_CR = CR.valueOf(b_prod).shiftRight(-eval_prec); if (b_prec.size() == n + 1) { // Add an n+1st slot. // Take care to make this exception-safe; b_prec and b_val // must remain consistent, even if we are interrupted, or run // out of memory. It's OK to just push on b_prec in that case. final CR next_b_as_CR = b_prod_as_CR.sqrt(); next_b = next_b_as_CR.get_appr(eval_prec); final BigInteger scaled_next_b = scale(next_b, -extra_eval_prec); b_prec.add(p); b_val.add(scaled_next_b); } else { // Reuse previous approximation to reduce sqrt iterations, // hopefully to one. final CR next_b_as_CR = new sqrt_CR(b_prod_as_CR, b_prec.get(n + 1), b_val.get(n + 1)); next_b = next_b_as_CR.get_appr(eval_prec); // We assume that set() doesn't throw for any reason. b_prec.set(n + 1, p); b_val.set(n + 1, scale(next_b, -extra_eval_prec)); } // b_prec.size() == b_val.size() >= n + 2 final BigInteger next_t = t.subtract(a_diff.multiply(a_diff) .shiftLeft(n + eval_prec)); // shift dist. usually neg. a = next_a; b = next_b; t = next_t; ++n; } final BigInteger sum = a.add(b); final BigInteger result = sum.multiply(sum).divide(t).shiftRight(2); return scale(result, -extra_eval_prec); } }