Home | History | Annotate | Download | only in creals
      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 (at) 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 (at) hp.com 12/17/2002
     98 // Added toStringFloatRep.      Hans.Boehm (at) hp.com 4/1/2004
     99 // Added get_appr() synchronization to allow access from multiple threads
    100 // hboehm (at) google.com 4/25/2014
    101 // Changed cos() prescaling to avoid logarithmic depth tree.
    102 // hboehm (at) google.com 6/30/2014
    103 // Added explicit asin() implementation.  Remove one.  Add ZERO and ONE and
    104 // make them public.  hboehm (at) google.com 5/21/2015
    105 
    106 package com.hp.creals;
    107 
    108 import 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 */
    164 public 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 */
    177 public 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 */
    189 public 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 */
    213 public 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 &lt; x</tt>, or +1 if <TT>this &gt; 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 (&gt; 0) included to the right of decimal point.
    629 *       @param  radix   Base ( &ge; 2, &le; 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)(1 << 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 //
    978 abstract 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.
    999 class 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.
   1012 class 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.
   1027 class 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
   1043 class 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.
   1056 class 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
   1070 class 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.
   1104 class 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.
   1151 class 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.
   1195 class 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.
   1235 class 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.
   1282 class 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)
   1325 class 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).
   1371 class 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 
   1453 class 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