Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 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 | // |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 23 | // Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED |
| 24 | // |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 28 | // |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 38 | // |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 51 | // |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 59 | // venue lying exclusively in Santa Clara County, California. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 60 | |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 93 | // Added valueOf(string, radix), fixed some documentation comments. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 94 | // Hans_Boehm@hp.com 1/12/2001 |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 95 | // Fixed a serious typo in inv_CR(): For negative arguments it produced |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 96 | // the wrong sign. This affected the sign of divisions. |
| 97 | // Added byteValue and fixed some comments. Hans.Boehm@hp.com 12/17/2002 |
| 98 | // Added toStringFloatRep. Hans.Boehm@hp.com 4/1/2004 |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 99 | // Added get_appr() synchronization to allow access from multiple threads |
| 100 | // hboehm@google.com 4/25/2014 |
| 101 | // Changed cos() prescaling to avoid logarithmic depth tree. |
| 102 | // hboehm@google.com 6/30/2014 |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 103 | // Added explicit asin() implementation. Remove one. Add ZERO and ONE and |
| 104 | // make them public. hboehm@google.com 5/21/2015 |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 105 | |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 106 | package com.hp.creals; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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> |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 154 | * Any operation may throw <TT>com.hp.creals.AbortedException</tt> if the thread in |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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> |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 158 | * Any operation may also throw <TT>com.hp.creals.PrecisionOverflowException</tt> |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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.) |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 162 | * |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 | |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 194 | // First some frequently used constants, so we don't have to |
| 195 | // recompute these all over the place. |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 196 | static final BigInteger big0 = BigInteger.ZERO; |
| 197 | static final BigInteger big1 = BigInteger.ONE; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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); |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 203 | static final BigInteger big10 = BigInteger.TEN; |
| 204 | static final BigInteger big750 = BigInteger.valueOf(750); |
| 205 | static final BigInteger bigm750 = BigInteger.valueOf(-750); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 206 | |
| 207 | /** |
| 208 | * Setting this to true requests that all computations be aborted by |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 209 | * throwing AbortedException. Must be rest to false before any further |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 219 | * Returns value / 2 ** precision rounded to an integer. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 225 | * Called only with the lock on the <TT>CR</tt> object |
| 226 | * already held. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 227 | */ |
| 228 | protected abstract BigInteger approximate(int precision); |
| 229 | transient int min_prec; |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 230 | // The smallest precision value with which the above |
| 231 | // has been called. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 232 | transient BigInteger max_appr; |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 233 | // The scaled approximation corresponding to min_prec. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 234 | transient boolean appr_valid = false; |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 235 | // min_prec and max_val are valid. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 236 | |
| 237 | // Helper functions |
| 238 | static int bound_log2(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 239 | int abs_n = Math.abs(n); |
| 240 | return (int)Math.ceil(Math.log((double)(abs_n + 1))/Math.log(2.0)); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 246 | // inside a function can generate an overflow. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 247 | static void check_prec(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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)) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 255 | throw new PrecisionOverflowException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 256 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 257 | } |
| 258 | |
| 259 | /** |
| 260 | * The constructive real number corresponding to a |
| 261 | * <TT>BigInteger</tt>. |
| 262 | */ |
| 263 | public static CR valueOf(BigInteger n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 264 | return new int_CR(n); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 265 | } |
| 266 | |
| 267 | /** |
| 268 | * The constructive real number corresponding to a |
| 269 | * Java <TT>int</tt>. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 270 | */ |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 271 | public static CR valueOf(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 272 | return valueOf(BigInteger.valueOf(n)); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 273 | } |
| 274 | |
| 275 | /** |
| 276 | * The constructive real number corresponding to a |
| 277 | * Java <TT>long</tt>. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 278 | */ |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 279 | public static CR valueOf(long n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 280 | return valueOf(BigInteger.valueOf(n)); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 287 | */ |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 288 | public static CR valueOf(double n) { |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 289 | if (Double.isNaN(n)) throw new ArithmeticException("Nan argument"); |
| 290 | if (Double.isInfinite(n)) throw new ArithmeticException("Infinite argument"); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 310 | */ |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 311 | public static CR valueOf(float n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 312 | return valueOf((double) n); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 313 | } |
| 314 | |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 315 | public static CR ZERO = valueOf(0); |
| 316 | public static CR ONE = valueOf(1); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 317 | |
| 318 | // Multiply k by 2**n. |
| 319 | static BigInteger shift(BigInteger k, int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 320 | if (n == 0) return k; |
| 321 | if (n < 0) return k.shiftRight(-n); |
| 322 | return k.shiftLeft(n); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 323 | } |
| 324 | |
| 325 | // Multiply by 2**n, rounding result |
| 326 | static BigInteger scale(BigInteger k, int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 327 | if (n >= 0) { |
| 328 | return k.shiftLeft(n); |
| 329 | } else { |
| 330 | BigInteger adj_k = shift(k, n+1).add(big1); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 331 | return adj_k.shiftRight(1); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 332 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 344 | */ |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 345 | public synchronized BigInteger get_appr(int precision) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 356 | } |
| 357 | |
| 358 | // Return the position of the msd. |
| 359 | // If x.msd() == n then |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 360 | // 2**(n-1) < abs(x) < 2**(n+1) |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 365 | int first_digit; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 375 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 376 | // This version may return Integer.MIN_VALUE if the correct |
| 377 | // answer is < n. |
| 378 | int msd(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 389 | } |
| 390 | |
| 391 | |
| 392 | // Functionally equivalent, but iteratively evaluates to higher |
| 393 | // precision. |
| 394 | int iter_msd(int n) |
| 395 | { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 396 | int prec = 0; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 397 | |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 402 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 403 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 411 | return iter_msd(Integer.MIN_VALUE); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 412 | } |
| 413 | |
| 414 | // A helper function for toString. |
| 415 | // Generate a String containing n zeroes. |
| 416 | private static String zeroes(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 417 | char[] a = new char[n]; |
| 418 | for (int i = 0; i < n; ++i) { |
| 419 | a[i] = '0'; |
| 420 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 421 | return new String(a); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 422 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 423 | |
| 424 | // Natural log of 2. Needed for some prescaling below. |
| 425 | // ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 426 | CR simple_ln() { |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 427 | return new prescaled_ln_CR(this.subtract(ONE)); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 437 | |
| 438 | // Atan of integer reciprocal. Used for PI. Could perhaps |
| 439 | // be made public. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 440 | static CR atan_reciprocal(int n) { |
| 441 | return new integral_atan_CR(n); |
| 442 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 443 | // Other constants used for PI computation. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 444 | static CR four = valueOf(4); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 454 | * @param x The other constructive real |
| 455 | * @param r Relative tolerance in bits |
| 456 | * @param a Absolute tolerance in bits |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 457 | */ |
| 458 | public int compareTo(CR x, int r, int a) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 | * |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 476 | * @param x The other constructive real |
| 477 | * @param a Absolute tolerance in bits |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 478 | */ |
| 479 | public int compareTo(CR x, int a) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 488 | } |
| 489 | |
| 490 | /** |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 491 | * Return -1 if <TT>this < x</tt>, or +1 if <TT>this > x</tt>. |
| 492 | * Should be called only if <TT>this != x</tt>. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 499 | for (int a = -20; ; a *= 2) { |
| 500 | check_prec(a); |
| 501 | int result = compareTo(x, a); |
| 502 | if (0 != result) return result; |
| 503 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 504 | } |
| 505 | |
| 506 | /** |
| 507 | * Equivalent to <TT>compareTo(CR.valueOf(0), a)</tt> |
| 508 | */ |
| 509 | public int signum(int a) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 515 | BigInteger this_appr = get_appr(needed_prec); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 516 | return this_appr.signum(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 517 | } |
| 518 | |
| 519 | /** |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 520 | * Return -1 if negative, +1 if positive. |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 521 | * Should be called only if <TT>this != 0</tt>. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 528 | for (int a = -20; ; a *= 2) { |
| 529 | check_prec(a); |
| 530 | int result = signum(a); |
| 531 | if (0 != result) return result; |
| 532 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 533 | } |
| 534 | |
| 535 | /** |
| 536 | * Return the constructive real number corresponding to the given |
| 537 | * textual representation and radix. |
| 538 | * |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 539 | * @param s [-] digit* [. digit*] |
| 540 | * @param radix |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 541 | */ |
| 542 | |
| 543 | public static CR valueOf(String s, int radix) |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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)); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 561 | } |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 562 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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 | * |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 567 | * @param n Number of digits (>= 0) included to the right of decimal point |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 568 | * @param radix Base ( >= 2, <= 16) for the resulting representation. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 569 | */ |
| 570 | public String toString(int n, int radix) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 599 | } |
| 600 | |
| 601 | |
| 602 | /** |
| 603 | * Equivalent to <TT>toString(n,10)</tt> |
| 604 | * |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 605 | * @param n Number of digits included to the right of decimal point |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 606 | */ |
| 607 | public String toString(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 608 | return toString(n, 10); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 609 | } |
| 610 | |
| 611 | /** |
| 612 | * Equivalent to <TT>toString(10, 10)</tt> |
| 613 | */ |
| 614 | public String toString() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 615 | return toString(10); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 616 | } |
| 617 | |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 618 | static double doubleLog2 = Math.log(2.0); |
| 619 | /** |
| 620 | * Return a textual scientific notation representation accurate |
| 621 | * to <TT>n</tt> places to the right of the decimal point. |
| 622 | * <TT>n</tt> must be nonnegative. A value smaller than |
| 623 | * <TT>radix</tt>**-<TT>m</tt> may be displayed as 0. |
| 624 | * The <TT>mantissa</tt> component of the result is either "0" |
| 625 | * or exactly <TT>n</tt> digits long. The <TT>sign</tt> |
| 626 | * component is zero exactly when the mantissa is "0". |
| 627 | * |
| 628 | * @param n Number of digits (> 0) included to the right of decimal point. |
| 629 | * @param radix Base ( ≥ 2, ≤ 16) for the resulting representation. |
| 630 | * @param m Precision used to distinguish number from zero. |
| 631 | * Expressed as a power of m. |
| 632 | */ |
| 633 | public StringFloatRep toStringFloatRep(int n, int radix, int m) { |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 634 | if (n <= 0) throw new ArithmeticException("Bad precision argument"); |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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) |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 640 | throw new PrecisionOverflowException(); |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 675 | /** |
| 676 | * Return a BigInteger which differs by less than one from the |
| 677 | * constructive real. |
| 678 | */ |
| 679 | public BigInteger BigIntegerValue() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 680 | return get_appr(0); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 688 | return BigIntegerValue().intValue(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 689 | } |
| 690 | |
| 691 | /** |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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 | /** |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 704 | return BigIntegerValue().longValue(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 720 | if (((orig_exp + exp_adj) & ~0x7ff) != 0) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 743 | return (float)doubleValue(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 755 | * @param n shift count, may be negative |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 756 | */ |
| 757 | public CR shiftLeft(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 758 | check_prec(n); |
| 759 | return new shifted_CR(this, n); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 760 | } |
| 761 | |
| 762 | /** |
| 763 | * Multiply a constructive real by 2**(-n). |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 764 | * @param n shift count, may be negative |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 765 | */ |
| 766 | public CR shiftRight(int n) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 767 | check_prec(n); |
| 768 | return new shifted_CR(this, -n); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 769 | } |
| 770 | |
| 771 | /** |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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 | /** |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 824 | return new select_CR(this, x, y); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 825 | } |
| 826 | |
| 827 | /** |
| 828 | * The maximum of two constructive reals. |
| 829 | */ |
| 830 | public CR max(CR x) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 831 | return subtract(x).select(x, this); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 832 | } |
| 833 | |
| 834 | /** |
| 835 | * The minimum of two constructive reals. |
| 836 | */ |
| 837 | public CR min(CR x) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 838 | return subtract(x).select(this, x); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 846 | return select(negate(), this); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 847 | } |
| 848 | |
| 849 | /** |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 850 | * The exponential function, that is e**<TT>this</tt>. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 851 | */ |
| 852 | public CR exp() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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)) |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 870 | .subtract(atan_reciprocal(239))); |
| 871 | // pi/4 = 4*atan(1/5) - atan(1/239) |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 872 | static CR half_pi = PI.shiftRight(1); |
| 873 | |
| 874 | /** |
| 875 | * The trigonometric cosine function. |
| 876 | */ |
| 877 | public CR cos() { |
Hans Boehm | a8e45a3 | 2015-02-25 14:53:13 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 881 | // Subtract multiples of PI |
Hans Boehm | a8e45a3 | 2015-02-25 14:53:13 -0800 | [diff] [blame] | 882 | BigInteger pi_multiples = scale(halfpi_multiples, -1); |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 883 | CR adjustment = PI.multiply(CR.valueOf(pi_multiples)); |
Hans Boehm | a8e45a3 | 2015-02-25 14:53:13 -0800 | [diff] [blame] | 884 | if (pi_multiples.and(big1).signum() != 0) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 885 | return subtract(adjustment).cos().negate(); |
| 886 | } else { |
| 887 | return subtract(adjustment).cos(); |
| 888 | } |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 889 | } else if (get_appr(-1).abs().compareTo(big2) >= 0) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 890 | // Scale further with double angle formula |
| 891 | CR cos_half = shiftRight(1).cos(); |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 892 | return cos_half.multiply(cos_half).shiftLeft(1).subtract(ONE); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 893 | } else { |
| 894 | return new prescaled_cos_CR(this); |
| 895 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 896 | } |
| 897 | |
| 898 | /** |
| 899 | * The trigonometric sine function. |
| 900 | */ |
| 901 | public CR sin() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 902 | return half_pi.subtract(this).cos(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 903 | } |
| 904 | |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 927 | static final BigInteger low_ln_limit = big8; /* sixteenths, i.e. 1/2 */ |
| 928 | static final BigInteger high_ln_limit = |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 929 | BigInteger.valueOf(16 + 8 /* 1.5 */); |
| 930 | static final BigInteger scaled_4 = |
| 931 | BigInteger.valueOf(4*16); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 932 | |
| 933 | /** |
| 934 | * The natural (base e) logarithm. |
| 935 | */ |
| 936 | public CR ln() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 937 | final int low_prec = -4; |
| 938 | BigInteger rough_appr = get_appr(low_prec); /* In sixteenths */ |
| 939 | if (rough_appr.compareTo(big0) < 0) { |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 940 | throw new ArithmeticException("ln(negative)"); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 956 | } |
| 957 | |
| 958 | /** |
| 959 | * The square root of a constructive real. |
| 960 | */ |
| 961 | public CR sqrt() { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 962 | return new sqrt_CR(this); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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; |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 981 | public synchronized BigInteger get_appr(int precision) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1002 | value = n; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1003 | } |
| 1004 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1005 | return scale(value, -p) ; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1006 | } |
| 1007 | } |
| 1008 | |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1031 | op1 = x; |
| 1032 | op2 = y; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1033 | } |
| 1034 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1047 | op = x; |
| 1048 | count = n; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1049 | } |
| 1050 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1051 | return op.get_appr(p - count); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1059 | op = x; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1060 | } |
| 1061 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1062 | return op.get_appr(p).negate(); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1063 | } |
| 1064 | } |
| 1065 | |
| 1066 | // Representation of: |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1067 | // op1 if selector < 0 |
| 1068 | // op2 if selector >= 0 |
| 1069 | // Assumes x = y if s = 0 |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1076 | selector = s; |
| 1077 | int selector_sign = selector.get_appr(-20).signum(); |
| 1078 | op1 = x; |
| 1079 | op2 = y; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1080 | } |
| 1081 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1108 | op1 = x; |
| 1109 | op2 = y; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1110 | } |
| 1111 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1112 | int half_prec = (p >> 1) - 1; |
| 1113 | int msd_op1 = op1.msd(half_prec); |
| 1114 | int msd_op2; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1115 | |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1140 | if (appr2.signum() == 0) return big0; |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1146 | } |
| 1147 | } |
| 1148 | |
Hans Boehm | 97767fa | 2014-11-26 17:16:40 -0800 | [diff] [blame] | 1149 | // Representation of the multiplicative inverse of a constructive |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1167 | // calculation. |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 1222 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1238 | op = x; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1239 | } |
| 1240 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 1265 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1301 | BigInteger big_op = BigInteger.valueOf(op); |
| 1302 | BigInteger big_op_squared = BigInteger.valueOf(op*op); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 1312 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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)); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1318 | current_sum = current_sum.add(current_term); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1319 | } |
| 1320 | return scale(current_sum, calc_precision - p); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 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; } |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1334 | // floating point. Unfortunately, other alternatives |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1335 | // appear to require precomputed tabular information. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1336 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 1357 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1358 | n += 1; |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1359 | current_sign = -current_sign; |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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); |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1366 | } |
| 1367 | } |
| 1368 | |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 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) { |
Hans Boehm | 9666c57 | 2015-06-19 18:27:45 -0700 | [diff] [blame] | 1422 | if (Thread.interrupted() || please_stop) throw new AbortedException(); |
Hans Boehm | b849a8e | 2015-05-20 17:53:27 -0700 | [diff] [blame] | 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 | |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1453 | class sqrt_CR extends CR { |
| 1454 | CR op; |
| 1455 | sqrt_CR(CR x) { op = x; } |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 1456 | final int fp_prec = 50; // Conservative estimate of number of |
| 1457 | // significant bits in double precision |
| 1458 | // computation. |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1459 | final int fp_op_prec = 60; |
| 1460 | protected BigInteger approximate(int p) { |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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(); |
Hans Boehm | 7b2383f | 2014-11-26 18:07:25 -0800 | [diff] [blame] | 1491 | if (scaled_appr < 0.0) |
| 1492 | throw new ArithmeticException("sqrt(negative)"); |
Hans Boehm | 36dc537 | 2014-11-26 13:35:59 -0800 | [diff] [blame] | 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 | } |
Hans Boehm | 349dbd7 | 2014-11-25 15:28:42 -0800 | [diff] [blame] | 1498 | } |
| 1499 | } |