1b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm/* 2b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * Copyright (C) 2015 The Android Open Source Project 3b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * 4b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * Licensed under the Apache License, Version 2.0 (the "License"); 5b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * you may not use this file except in compliance with the License. 6b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * You may obtain a copy of the License at 7b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * 8b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * http://www.apache.org/licenses/LICENSE-2.0 9b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * 10b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * Unless required by applicable law or agreed to in writing, software 11b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * distributed under the License is distributed on an "AS IS" BASIS, 12b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * See the License for the specific language governing permissions and 14b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * limitations under the License. 15b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm */ 16b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 17b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm/* 18b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * The above license covers additions and changes by AOSP authors. 19b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm * The original code is licensed as follows: 20b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm */ 21b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 22b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm// 2336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// Copyright (c) 1999, Silicon Graphics, Inc. -- ALL RIGHTS RESERVED 2436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// 25349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Permission is granted free of charge to copy, modify, use and distribute 26349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// this software provided you include the entirety of this notice in all 27349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// copies made. 2836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// 29349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY 30349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, 31349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT 32349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// FOR A PARTICULAR PURPOSE OR NON-INFRINGING. SGI ASSUMES NO RISK AS TO THE 33349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// QUALITY AND PERFORMANCE OF THE SOFTWARE. SHOULD THE SOFTWARE PROVE 34349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// DEFECTIVE IN ANY RESPECT, SGI ASSUMES NO COST OR LIABILITY FOR ANY 35349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES 36349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS 37349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. 3836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// 39349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, 40349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR 41349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// OTHERWISE, SHALL SGI BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, 42349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE 43349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK 44349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL 45349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF SGI SHALL HAVE BEEN INFORMED OF 46349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// THE POSSIBILITY OF SUCH DAMAGES. THIS LIMITATION OF LIABILITY SHALL NOT 47349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// APPLY TO LIABILITY RESULTING FROM SGI's NEGLIGENCE TO THE EXTENT APPLICABLE 48349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE 49349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT 50349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. 5136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// 52349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// These license terms shall be governed by and construed in accordance with 53349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// the laws of the United States and the State of California as applied to 54349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// agreements entered into and to be performed entirely within California 55349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// between California residents. Any litigation relating to these license 56349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// terms shall be subject to the exclusive jurisdiction of the Federal Courts 57349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// of the Northern District of California (or, absent subject matter 58349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// jurisdiction in such courts, the courts of the State of California), with 5936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// venue lying exclusively in Santa Clara County, California. 60349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 6197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Copyright (c) 2001-2004, Hewlett-Packard Development Company, L.P. 6297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// 6397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Permission is granted free of charge to copy, modify, use and distribute 6497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// this software provided you include the entirety of this notice in all 6597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// copies made. 6697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// 6797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// THIS SOFTWARE IS PROVIDED ON AN AS IS BASIS, WITHOUT WARRANTY OF ANY 6897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, 6997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// WARRANTIES THAT THE SUBJECT SOFTWARE IS FREE OF DEFECTS, MERCHANTABLE, FIT 7097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// FOR A PARTICULAR PURPOSE OR NON-INFRINGING. HEWLETT-PACKARD ASSUMES 7197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// NO RISK AS TO THE QUALITY AND PERFORMANCE OF THE SOFTWARE. 7297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// SHOULD THE SOFTWARE PROVE DEFECTIVE IN ANY RESPECT, 7397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// HEWLETT-PACKARD ASSUMES NO COST OR LIABILITY FOR ANY 7497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES 7597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// AN ESSENTIAL PART OF THIS LICENSE. NO USE OF ANY SUBJECT SOFTWARE IS 7697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// AUTHORIZED HEREUNDER EXCEPT UNDER THIS DISCLAIMER. 7797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// 7897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, WHETHER TORT (INCLUDING, 7997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// WITHOUT LIMITATION, NEGLIGENCE OR STRICT LIABILITY), CONTRACT, OR 8097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// OTHERWISE, SHALL HEWLETT-PACKARD BE LIABLE FOR ANY DIRECT, INDIRECT, SPECIAL, 8197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// INCIDENTAL, OR CONSEQUENTIAL DAMAGES OF ANY CHARACTER WITH RESPECT TO THE 8297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// SOFTWARE INCLUDING, WITHOUT LIMITATION, DAMAGES FOR LOSS OF GOODWILL, WORK 8397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// STOPPAGE, LOSS OF DATA, COMPUTER FAILURE OR MALFUNCTION, OR ANY AND ALL 8497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// OTHER COMMERCIAL DAMAGES OR LOSSES, EVEN IF HEWLETT-PACKARD SHALL 8597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// HAVE BEEN INFORMED OF THE POSSIBILITY OF SUCH DAMAGES. 8697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// THIS LIMITATION OF LIABILITY SHALL NOT APPLY TO LIABILITY RESULTING 8797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// FROM HEWLETT-PACKARD's NEGLIGENCE TO THE EXTENT APPLICABLE 8897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// LAW PROHIBITS SUCH LIMITATION. SOME JURISDICTIONS DO NOT ALLOW THE 8997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// EXCLUSION OR LIMITATION OF INCIDENTAL OR CONSEQUENTIAL DAMAGES, SO THAT 9097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// EXCLUSION AND LIMITATION MAY NOT APPLY TO YOU. 9197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// 9297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm 93349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Added valueOf(string, radix), fixed some documentation comments. 9436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// Hans_Boehm@hp.com 1/12/2001 95349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Fixed a serious typo in inv_CR(): For negative arguments it produced 9697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// the wrong sign. This affected the sign of divisions. 9797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Added byteValue and fixed some comments. Hans.Boehm@hp.com 12/17/2002 9897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Added toStringFloatRep. Hans.Boehm@hp.com 4/1/2004 997b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm// Added get_appr() synchronization to allow access from multiple threads 1007b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm// hboehm@google.com 4/25/2014 1017b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm// Changed cos() prescaling to avoid logarithmic depth tree. 1027b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm// hboehm@google.com 6/30/2014 103b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm// Added explicit asin() implementation. Remove one. Add ZERO and ONE and 104b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm// make them public. hboehm@google.com 5/21/2015 105349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 10697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehmpackage com.hp.creals; 107349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 108349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmimport java.math.BigInteger; 109349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 110349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 111349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Constructive real numbers, also known as recursive, or computable reals. 112349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Each recursive real number is represented as an object that provides an 113349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* approximation function for the real number. 114349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The approximation function guarantees that the generated approximation 115349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* is accurate to the specified precision. 116349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Arithmetic operations on constructive reals produce new such objects; 117349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* they typically do not perform any real computation. 118349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* In this sense, arithmetic computations are exact: They produce 119349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* a description which describes the exact answer, and can be used to 120349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* later approximate it to arbitrary precision. 121349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 122349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* When approximations are generated, <I>e.g.</i> for output, they are 123349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* accurate to the requested precision; no cumulative rounding errors 124349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* are visible. 125349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* In order to achieve this precision, the approximation function will often 126349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* need to approximate subexpressions to greater precision than was originally 127349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* demanded. Thus the approximation of a constructive real number 128349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* generated through a complex sequence of operations may eventually require 129349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* evaluation to very high precision. This usually makes such computations 130349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* prohibitively expensive for large numerical problems. 131349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* But it is perfectly appropriate for use in a desk calculator, 132349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* for small numerical problems, for the evaluation of expressions 133349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* computated by a symbolic algebra system, for testing of accuracy claims 134349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* for floating point code on small inputs, or the like. 135349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 136349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* We expect that the vast majority of uses will ignore the particular 137349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* implementation, and the member functons <TT>approximate</tt> 138349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* and <TT>get_appr</tt>. Such applications will treat <TT>CR</tt> as 139349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* a conventional numerical type, with an interface modelled on 140349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <TT>java.math.BigInteger</tt>. No subclasses of <TT>CR</tt> 141349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* will be explicitly mentioned by such a program. 142349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 143349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* All standard arithmetic operations, as well as a few algebraic 144349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* and transcendal functions are provided. Constructive reals are 145349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* immutable; thus all of these operations return a new constructive real. 146349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 147349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* A few uses will require explicit construction of approximation functions. 148349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The requires the construction of a subclass of <TT>CR</tt> with 149349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* an overridden <TT>approximate</tt> function. Note that <TT>approximate</tt> 150349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* should only be defined, but never called. <TT>get_appr</tt> 151349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* provides the same functionality, but adds the caching necessary to obtain 152349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* reasonable performance. 153349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 1549666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm* Any operation may throw <TT>com.hp.creals.AbortedException</tt> if the thread in 155349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* which it is executing is interrupted. (<TT>InterruptedException</tt> cannot 156349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* be used for this purpose, since CR inherits from <TT>Number</tt>.) 157349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <P> 1589666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm* Any operation may also throw <TT>com.hp.creals.PrecisionOverflowException</tt> 159349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* If the precision request generated during any subcalculation overflows 160349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* a 28-bit integer. (This should be extremely unlikely, except as an 161349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* outcome of a division by zero, or other erroneous computation.) 16236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* 163349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 164349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmpublic abstract class CR extends Number { 165349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // CR is the basic representation of a number. 166349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Abstractly this is a function for computing an approximation 167349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // plus the current best approximation. 168349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // We could do without the latter, but that would 169349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // be atrociously slow. 170349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1719666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm/** 1729666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * Indicates a constructive real operation was interrupted. 1739666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * Most constructive real operations may throw such an exception. 1749666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * This is unchecked, since Number methods may not raise checked 1759666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * exceptions. 1769666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm*/ 1779666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehmpublic static class AbortedException extends RuntimeException { 1789666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm public AbortedException() { super(); } 1799666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm public AbortedException(String s) { super(s); } 1809666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm} 1819666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm 1829666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm/** 1839666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * Indicates that the number of bits of precision requested by 1849666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * a computation on constructive reals required more than 28 bits, 1859666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * and was thus in danger of overflowing an int. 1869666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * This is likely to be a symptom of a diverging computation, 1879666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm * <I>e.g.</i> division by zero. 1889666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm*/ 1899666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehmpublic static class PrecisionOverflowException extends RuntimeException { 1909666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm public PrecisionOverflowException() { super(); } 1919666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm public PrecisionOverflowException(String s) { super(s); } 1929666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm} 1939666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm 194349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // First some frequently used constants, so we don't have to 195349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // recompute these all over the place. 196b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm static final BigInteger big0 = BigInteger.ZERO; 197b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm static final BigInteger big1 = BigInteger.ONE; 198349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger bigm1 = BigInteger.valueOf(-1); 199349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger big2 = BigInteger.valueOf(2); 200349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger big3 = BigInteger.valueOf(3); 201349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger big6 = BigInteger.valueOf(6); 202349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger big8 = BigInteger.valueOf(8); 203b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm static final BigInteger big10 = BigInteger.TEN; 204b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm static final BigInteger big750 = BigInteger.valueOf(750); 205b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm static final BigInteger bigm750 = BigInteger.valueOf(-750); 206349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 207349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 208349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Setting this to true requests that all computations be aborted by 2099666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm* throwing AbortedException. Must be rest to false before any further 210349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* computation. Ideally Thread.interrupt() should be used instead, but 211349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* that doesn't appear to be consistently supported by browser VMs. 212349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 213349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmpublic volatile static boolean please_stop = false; 214349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 215349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 216349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Must be defined in subclasses of <TT>CR</tt>. 217349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Most users can ignore the existence of this method, and will 218349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* not ever need to define a <TT>CR</tt> subclass. 21997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Returns value / 2 ** precision rounded to an integer. 220349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The error in the result is strictly < 1. 221349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Informally, approximate(n) gives a scaled approximation 222349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* accurate to 2**n. 223349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Implementations may safely assume that precision is 224349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* at least a factor of 8 away from overflow. 2257b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm* Called only with the lock on the <TT>CR</tt> object 2267b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm* already held. 227349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 228349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected abstract BigInteger approximate(int precision); 229349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm transient int min_prec; 23036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // The smallest precision value with which the above 23136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // has been called. 232349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm transient BigInteger max_appr; 23336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // The scaled approximation corresponding to min_prec. 234349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm transient boolean appr_valid = false; 23536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // min_prec and max_val are valid. 236349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 237349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Helper functions 238349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static int bound_log2(int n) { 23936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int abs_n = Math.abs(n); 24036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return (int)Math.ceil(Math.log((double)(abs_n + 1))/Math.log(2.0)); 241349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 242349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Check that a precision is at least a factor of 8 away from 243349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // overflowng the integer used to hold a precision spec. 244349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // We generally perform this check early on, and then convince 245349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // ourselves that none of the operations performed on precisions 246b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // inside a function can generate an overflow. 247349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static void check_prec(int n) { 24836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int high = n >> 28; 24936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // if n is not in danger of overflowing, then the 4 high order 25036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // bits should be identical. Thus high is either 0 or -1. 25136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // The rest of this is to test for either of those in a way 25236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // that should be as cheap as possible. 25336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int high_shifted = n >> 29; 25436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (0 != (high ^ high_shifted)) { 2559666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm throw new PrecisionOverflowException(); 25636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 257349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 258349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 259349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 260349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The constructive real number corresponding to a 261349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <TT>BigInteger</tt>. 262349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 263349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(BigInteger n) { 26436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new int_CR(n); 265349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 266349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 267349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 268349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The constructive real number corresponding to a 269349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Java <TT>int</tt>. 27036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm*/ 271349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(int n) { 27236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return valueOf(BigInteger.valueOf(n)); 273349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 274349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 275349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 276349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The constructive real number corresponding to a 277349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Java <TT>long</tt>. 27836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm*/ 279349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(long n) { 28036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return valueOf(BigInteger.valueOf(n)); 281349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 282349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 283349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 284349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The constructive real number corresponding to a 285349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Java <TT>double</tt>. 286349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The result is undefined if argument is infinite or NaN. 28736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm*/ 288349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(double n) { 2897b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm if (Double.isNaN(n)) throw new ArithmeticException("Nan argument"); 2907b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm if (Double.isInfinite(n)) throw new ArithmeticException("Infinite argument"); 29136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm boolean negative = (n < 0.0); 29236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm long bits = Double.doubleToLongBits(Math.abs(n)); 29336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm long mantissa = (bits & 0xfffffffffffffL); 29436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int biased_exp = (int)(bits >> 52); 29536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int exp = biased_exp - 1075; 29636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (biased_exp != 0) { 29736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm mantissa += (1L << 52); 29836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 29936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm mantissa <<= 1; 30036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 30136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR result = valueOf(mantissa).shiftLeft(exp); 30236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (negative) result = result.negate(); 30336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result; 304349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 305349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 306349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 307349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The constructive real number corresponding to a 308349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Java <TT>float</tt>. 309349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The result is undefined if argument is infinite or NaN. 31036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm*/ 311349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(float n) { 31236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return valueOf((double) n); 313349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 314349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 315b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm public static CR ZERO = valueOf(0); 316b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm public static CR ONE = valueOf(1); 317349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 318349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Multiply k by 2**n. 319349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static BigInteger shift(BigInteger k, int n) { 32036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (n == 0) return k; 32136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (n < 0) return k.shiftRight(-n); 32236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return k.shiftLeft(n); 323349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 324349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 325349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Multiply by 2**n, rounding result 326349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static BigInteger scale(BigInteger k, int n) { 32736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (n >= 0) { 32836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return k.shiftLeft(n); 32936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 33036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger adj_k = shift(k, n+1).add(big1); 331349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return adj_k.shiftRight(1); 33236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 333349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 334349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 335349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Identical to approximate(), but maintain and update cache. 336349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 337349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Returns value / 2 ** prec rounded to an integer. 338349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The error in the result is strictly < 1. 339349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Produces the same answer as <TT>approximate</tt>, but uses and 340349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* maintains a cached approximation. 341349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Normally not overridden, and called only from <TT>approximate</tt> 342349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* methods in subclasses. Not needed if the provided operations 343349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* on constructive reals suffice. 34436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm*/ 3457b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm public synchronized BigInteger get_appr(int precision) { 34636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(precision); 34736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (appr_valid && precision >= min_prec) { 34836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(max_appr, min_prec - precision); 34936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 35036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger result = approximate(precision); 35136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm min_prec = precision; 35236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm max_appr = result; 35336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm appr_valid = true; 35436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result; 35536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 356349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 357349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 358349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Return the position of the msd. 359349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // If x.msd() == n then 36036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2**(n-1) < abs(x) < 2**(n+1) 361349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // This initial version assumes that max_appr is valid 362349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // and sufficiently removed from zero 363349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // that the msd is determined. 364349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int known_msd() { 36536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int first_digit; 366349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int length; 367349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm if (max_appr.signum() >= 0) { 368349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm length = max_appr.bitLength(); 369349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } else { 370349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm length = max_appr.negate().bitLength(); 371349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 372349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm first_digit = min_prec + length - 1; 373349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return first_digit; 374349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 37536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm 376349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // This version may return Integer.MIN_VALUE if the correct 377349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // answer is < n. 378349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int msd(int n) { 37936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (!appr_valid || 38036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm max_appr.compareTo(big1) <= 0 38136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm && max_appr.compareTo(bigm1) >= 0) { 38236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm get_appr(n - 1); 38336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (max_appr.abs().compareTo(big1) <= 0) { 38436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // msd could still be arbitrarily far to the right. 38536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return Integer.MIN_VALUE; 38636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 38736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 38836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return known_msd(); 389349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 390349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 391349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 392349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Functionally equivalent, but iteratively evaluates to higher 393349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // precision. 394349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int iter_msd(int n) 395349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm { 39636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int prec = 0; 39736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm 39836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm for (;prec > n + 30; prec = (prec * 3)/2 - 16) { 39936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int msd = msd(prec); 40036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (msd != Integer.MIN_VALUE) return msd; 40136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(prec); 4029666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 40336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 404349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return msd(n); 405349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 406349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 407349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // This version returns a correct answer eventually, except 408349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // that it loops forever (or throws an exception when the 409349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // requested precision overflows) if this constructive real is zero. 410349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int msd() { 41136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return iter_msd(Integer.MIN_VALUE); 412349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 413349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 414349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // A helper function for toString. 415349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Generate a String containing n zeroes. 416349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm private static String zeroes(int n) { 41736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm char[] a = new char[n]; 41836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm for (int i = 0; i < n; ++i) { 41936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm a[i] = '0'; 42036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 421349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new String(a); 42236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 423349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 424349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Natural log of 2. Needed for some prescaling below. 425349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // ln(2) = 7ln(10/9) - 2ln(25/24) + 3ln(81/80) 42636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR simple_ln() { 427b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return new prescaled_ln_CR(this.subtract(ONE)); 42836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 42936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR ten_ninths = valueOf(10).divide(valueOf(9)); 43036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR twentyfive_twentyfourths = valueOf(25).divide(valueOf(24)); 43136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR eightyone_eightyeths = valueOf(81).divide(valueOf(80)); 43236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR ln2_1 = valueOf(7).multiply(ten_ninths.simple_ln()); 43336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR ln2_2 = 43436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm valueOf(2).multiply(twentyfive_twentyfourths.simple_ln()); 43536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR ln2_3 = valueOf(3).multiply(eightyone_eightyeths.simple_ln()); 43636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR ln2 = ln2_1.subtract(ln2_2).add(ln2_3); 437349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 438349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Atan of integer reciprocal. Used for PI. Could perhaps 439349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // be made public. 44036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR atan_reciprocal(int n) { 44136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new integral_atan_CR(n); 44236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 443349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Other constants used for PI computation. 44436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static CR four = valueOf(4); 445349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 446349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // Public operations. 447349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 448349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return 0 if x = y to within the indicated tolerance, 449349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* -1 if x < y, and +1 if x > y. If x and y are indeed 450349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* equal, it is guaranteed that 0 will be returned. If 451349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* they differ by less than the tolerance, anything 452349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* may happen. The tolerance allowed is 453349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* the maximum of (abs(this)+abs(x))*(2**r) and 2**a 45436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param x The other constructive real 45536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param r Relative tolerance in bits 45636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param a Absolute tolerance in bits 457349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 458349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int compareTo(CR x, int r, int a) { 45936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int this_msd = iter_msd(a); 46036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int x_msd = x.iter_msd(this_msd > a? this_msd : a); 46136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int max_msd = (x_msd > this_msd? x_msd : this_msd); 46236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int rel = max_msd + r; 46336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // This can't approach overflow, since r and a are 46436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // effectively divided by 2, and msds are checked. 46536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int abs_prec = (rel > a? rel : a); 46636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return compareTo(x, abs_prec); 467349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 468349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 469349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 470349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Approximate comparison with only an absolute tolerance. 471349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Identical to the three argument version, but without a relative 472349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* tolerance. 473349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Result is 0 if both constructive reals are equal, indeterminate 474349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* if they differ by less than 2**a. 475349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* 47636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param x The other constructive real 47736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param a Absolute tolerance in bits 478349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 479349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int compareTo(CR x, int a) { 48036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int needed_prec = a - 1; 48136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger this_appr = get_appr(needed_prec); 48236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger x_appr = x.get_appr(needed_prec); 48336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int comp1 = this_appr.compareTo(x_appr.add(big1)); 48436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (comp1 > 0) return 1; 48536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int comp2 = this_appr.compareTo(x_appr.subtract(big1)); 48636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (comp2 < 0) return -1; 48736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return 0; 488349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 489349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 490349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 49197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Return -1 if <TT>this < x</tt>, or +1 if <TT>this > x</tt>. 49297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Should be called only if <TT>this != x</tt>. 493349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* If <TT>this == x</tt>, this will not terminate correctly; typically it 494349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* will run until it exhausts memory. 495349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* If the two constructive reals may be equal, the two or 3 argument 496349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* version of compareTo should be used. 497349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 498349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int compareTo(CR x) { 49936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm for (int a = -20; ; a *= 2) { 50036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(a); 50136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int result = compareTo(x, a); 50236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (0 != result) return result; 50336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 504349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 505349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 506349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 507349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Equivalent to <TT>compareTo(CR.valueOf(0), a)</tt> 508349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 509349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int signum(int a) { 51036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (appr_valid) { 51136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int quick_try = max_appr.signum(); 51236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (0 != quick_try) return quick_try; 51336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 51436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int needed_prec = a - 1; 515349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm BigInteger this_appr = get_appr(needed_prec); 51636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return this_appr.signum(); 517349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 518349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 519349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 520349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return -1 if negative, +1 if positive. 52197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Should be called only if <TT>this != 0</tt>. 522349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* In the 0 case, this will not terminate correctly; typically it 523349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* will run until it exhausts memory. 524349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* If the two constructive reals may be equal, the one or two argument 525349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* version of signum should be used. 526349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 527349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int signum() { 52836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm for (int a = -20; ; a *= 2) { 52936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(a); 53036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int result = signum(a); 53136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (0 != result) return result; 53236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 533349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 534349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 535349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 536349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return the constructive real number corresponding to the given 537349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* textual representation and radix. 538349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* 53936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param s [-] digit* [. digit*] 54036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param radix 541349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 542349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 543349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR valueOf(String s, int radix) 54436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm throws NumberFormatException { 54536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int len = s.length(); 54636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int start_pos = 0, point_pos; 54736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String fraction; 54836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (s.charAt(start_pos) == ' ') ++start_pos; 54936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (s.charAt(len - 1) == ' ') --len; 55036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm point_pos = s.indexOf('.', start_pos); 55136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (point_pos == -1) { 55236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm point_pos = len; 55336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm fraction = "0"; 55436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 55536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm fraction = s.substring(point_pos + 1, len); 55636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 55736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String whole = s.substring(start_pos, point_pos); 55836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_result = new BigInteger(whole + fraction, radix); 55936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger divisor = BigInteger.valueOf(radix).pow(fraction.length()); 56036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return CR.valueOf(scaled_result).divide(CR.valueOf(divisor)); 561349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 56236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm 563349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 564349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return a textual representation accurate to <TT>n</tt> places 565349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* to the right of the decimal point. <TT>n</tt> must be nonnegative. 566349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* 56797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* @param n Number of digits (>= 0) included to the right of decimal point 56836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param radix Base ( >= 2, <= 16) for the resulting representation. 569349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 570349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public String toString(int n, int radix) { 57136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR scaled_CR; 57236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (16 == radix) { 57336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm scaled_CR = shiftLeft(4*n); 57436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 57536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scale_factor = BigInteger.valueOf(radix).pow(n); 57636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm scaled_CR = multiply(new int_CR(scale_factor)); 57736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 57836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_int = scaled_CR.get_appr(0); 57936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String scaled_string = scaled_int.abs().toString(radix); 58036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String result; 58136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (0 == n) { 58236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm result = scaled_string; 58336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 58436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int len = scaled_string.length(); 58536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (len <= n) { 58636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Add sufficient leading zeroes 58736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String z = zeroes(n + 1 - len); 58836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm scaled_string = z + scaled_string; 58936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm len = n + 1; 59036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 59136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String whole = scaled_string.substring(0, len - n); 59236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm String fraction = scaled_string.substring(len - n); 59336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm result = whole + "." + fraction; 59436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 59536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (scaled_int.signum() < 0) { 59636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm result = "-" + result; 59736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 59836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result; 599349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 600349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 601349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 602349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 603349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Equivalent to <TT>toString(n,10)</tt> 604349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* 60536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param n Number of digits included to the right of decimal point 606349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 607349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public String toString(int n) { 60836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return toString(n, 10); 609349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 610349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 611349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 612349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Equivalent to <TT>toString(10, 10)</tt> 613349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 614349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public String toString() { 61536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return toString(10); 616349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 617349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 61897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm static double doubleLog2 = Math.log(2.0); 61997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm/** 62097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Return a textual scientific notation representation accurate 62197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* to <TT>n</tt> places to the right of the decimal point. 62297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* <TT>n</tt> must be nonnegative. A value smaller than 62397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* <TT>radix</tt>**-<TT>m</tt> may be displayed as 0. 62497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* The <TT>mantissa</tt> component of the result is either "0" 62597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* or exactly <TT>n</tt> digits long. The <TT>sign</tt> 62697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* component is zero exactly when the mantissa is "0". 62797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* 62897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* @param n Number of digits (> 0) included to the right of decimal point. 62997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* @param radix Base ( ≥ 2, ≤ 16) for the resulting representation. 63097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* @param m Precision used to distinguish number from zero. 63197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Expressed as a power of m. 63297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm*/ 63397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm public StringFloatRep toStringFloatRep(int n, int radix, int m) { 6347b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm if (n <= 0) throw new ArithmeticException("Bad precision argument"); 63597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm double log2_radix = Math.log((double)radix)/doubleLog2; 63697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm BigInteger big_radix = BigInteger.valueOf(radix); 63797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm long long_msd_prec = (long)(log2_radix * (double)m); 63897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm if (long_msd_prec > (long)Integer.MAX_VALUE 63997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm || long_msd_prec < (long)Integer.MIN_VALUE) 6409666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm throw new PrecisionOverflowException(); 64197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm int msd_prec = (int)long_msd_prec; 64297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm check_prec(msd_prec); 64397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm int msd = iter_msd(msd_prec - 2); 64497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm if (msd == Integer.MIN_VALUE) 64597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return new StringFloatRep(0, "0", radix, 0); 64697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm int exponent = (int)Math.ceil((double)msd / log2_radix); 64797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm // Guess for the exponent. Try to get it usually right. 64897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm int scale_exp = exponent - n; 64997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm CR scale; 65097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm if (scale_exp > 0) { 65197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scale = CR.valueOf(big_radix.pow(scale_exp)).inverse(); 65297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } else { 65397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scale = CR.valueOf(big_radix.pow(-scale_exp)); 65497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 65597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm CR scaled_res = multiply(scale); 65697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm BigInteger scaled_int = scaled_res.get_appr(0); 65797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm int sign = scaled_int.signum(); 65897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm String scaled_string = scaled_int.abs().toString(radix); 65997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm while (scaled_string.length() < n) { 66097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm // exponent was too large. Adjust. 66197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scaled_res = scaled_res.multiply(CR.valueOf(big_radix)); 66297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm exponent -= 1; 66397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scaled_int = scaled_res.get_appr(0); 66497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm sign = scaled_int.signum(); 66597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scaled_string = scaled_int.abs().toString(radix); 66697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 66797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm if (scaled_string.length() > n) { 66897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm // exponent was too small. Adjust by truncating. 66997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm exponent += (scaled_string.length() - n); 67097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm scaled_string = scaled_string.substring(0, n); 67197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 67297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return new StringFloatRep(sign, scaled_string, radix, exponent); 67397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 67497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm 675349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 676349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return a BigInteger which differs by less than one from the 677349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* constructive real. 678349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 679349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public BigInteger BigIntegerValue() { 68036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return get_appr(0); 681349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 682349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 683349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 684349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return an int which differs by less than one from the 685349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* constructive real. Behavior on overflow is undefined. 686349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 687349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public int intValue() { 68836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return BigIntegerValue().intValue(); 689349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 690349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 691349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 69297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Return an int which differs by less than one from the 69397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* constructive real. Behavior on overflow is undefined. 69497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm*/ 69597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm public byte byteValue() { 69697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return BigIntegerValue().byteValue(); 69797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 69897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm 69997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm/** 700349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return a long which differs by less than one from the 701349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* constructive real. Behavior on overflow is undefined. 702349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 703349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public long longValue() { 70436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return BigIntegerValue().longValue(); 705349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 706349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 707349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 708349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return a double which differs by less than one in the least 709349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* represented bit from the constructive real. 710349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 711349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public double doubleValue() { 71236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int my_msd = iter_msd(-1080 /* slightly > exp. range */); 71336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (Integer.MIN_VALUE == my_msd) return 0.0; 71436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int needed_prec = my_msd - 60; 71536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm double scaled_int = get_appr(needed_prec).doubleValue(); 71636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm boolean may_underflow = (needed_prec < -1000); 71736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm long scaled_int_rep = Double.doubleToLongBits(scaled_int); 71836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm long exp_adj = may_underflow? needed_prec + 96 : needed_prec; 71936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm long orig_exp = (scaled_int_rep >> 52) & 0x7ff; 720349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm if (((orig_exp + exp_adj) & ~0x7ff) != 0) { 72136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // overflow 72236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (scaled_int < 0.0) { 72336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return Double.NEGATIVE_INFINITY; 72436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 72536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return Double.POSITIVE_INFINITY; 72636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 72736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 72836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm scaled_int_rep += exp_adj << 52; 72936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm double result = Double.longBitsToDouble(scaled_int_rep); 73036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (may_underflow) { 73136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm double two48 = (double)(1 << 48); 73236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result/two48/two48; 73336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 73436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result; 73536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 736349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 737349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 738349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 739349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Return a float which differs by less than one in the least 740349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* represented bit from the constructive real. 741349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 742349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public float floatValue() { 74336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return (float)doubleValue(); 744349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 745349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 746349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 747349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Add two constructive reals. 748349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 749349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR add(CR x) { 750349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new add_CR(this, x); 751349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 752349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 753349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 754349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Multiply a constructive real by 2**n. 75536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param n shift count, may be negative 756349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 757349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR shiftLeft(int n) { 75836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(n); 75936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new shifted_CR(this, n); 760349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 761349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 762349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 763349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Multiply a constructive real by 2**(-n). 76436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm* @param n shift count, may be negative 765349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 766349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR shiftRight(int n) { 76736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(n); 76836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new shifted_CR(this, -n); 769349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 770349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 771349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 77297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* Produce a constructive real equivalent to the original, assuming 77397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* the original was an integer. Undefined results if the original 77497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* was not an integer. Prevents evaluation of digits to the right 77597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* of the decimal point, and may thus improve performance. 77697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm*/ 77797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm public CR assumeInt() { 77897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return new assumed_int_CR(this); 77997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 78097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm 78197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm/** 782349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The additive inverse of a constructive real 783349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 784349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR negate() { 785349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new neg_CR(this); 786349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 787349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 788349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 789349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The difference between two constructive reals 790349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 791349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR subtract(CR x) { 792349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new add_CR(this, x.negate()); 793349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 794349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 795349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 796349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The product of two constructive reals 797349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 798349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR multiply(CR x) { 799349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new mult_CR(this, x); 800349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 801349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 802349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 803349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The multiplicative inverse of a constructive real. 804349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* <TT>x.inverse()</tt> is equivalent to <TT>CR.valueOf(1).divide(x)</tt>. 805349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 806349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR inverse() { 807349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new inv_CR(this); 808349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 809349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 810349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 811349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The quotient of two constructive reals. 812349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 813349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR divide(CR x) { 814349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm return new mult_CR(this, x.inverse()); 815349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 816349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 817349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 818349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The real number <TT>x</tt> if <TT>this</tt> < 0, or <TT>y</tt> otherwise. 819349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Requires <TT>x</tt> = <TT>y</tt> if <TT>this</tt> = 0. 820349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Since comparisons may diverge, this is often 821349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* a useful alternative to conditionals. 822349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 823349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR select(CR x, CR y) { 82436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new select_CR(this, x, y); 825349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 826349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 827349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 828349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The maximum of two constructive reals. 829349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 830349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR max(CR x) { 83136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return subtract(x).select(x, this); 832349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 833349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 834349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 835349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The minimum of two constructive reals. 836349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 837349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR min(CR x) { 83836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return subtract(x).select(this, x); 839349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 840349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 841349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 842349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The absolute value of a constructive reals. 843349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* Note that this cannot be written as a conditional. 844349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 845349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR abs() { 84636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return select(negate(), this); 847349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 848349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 849349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 85097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm* The exponential function, that is e**<TT>this</tt>. 851349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 852349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR exp() { 85336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm final int low_prec = -10; 85436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger rough_appr = get_appr(low_prec); 85536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.signum() < 0) return negate().exp().inverse(); 85636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.compareTo(big2) > 0) { 85736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR square_root = shiftRight(1).exp(); 85836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return square_root.multiply(square_root); 85936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 86036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new prescaled_exp_CR(this); 861349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 862349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 863349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 864349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static CR two = valueOf(2); 865349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 866349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 867349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The ratio of a circle's circumference to its diameter. 868349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 869349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public static CR PI = four.multiply(four.multiply(atan_reciprocal(5)) 87036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm .subtract(atan_reciprocal(239))); 87136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // pi/4 = 4*atan(1/5) - atan(1/239) 872349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static CR half_pi = PI.shiftRight(1); 873349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 874349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 875349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The trigonometric cosine function. 876349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 877349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR cos() { 878a8e45a3255089266ccf6500188a2d31af8634e4aHans Boehm BigInteger halfpi_multiples = divide(PI).get_appr(-1); 879a8e45a3255089266ccf6500188a2d31af8634e4aHans Boehm BigInteger abs_halfpi_multiples = halfpi_multiples.abs(); 880a8e45a3255089266ccf6500188a2d31af8634e4aHans Boehm if (abs_halfpi_multiples.compareTo(big2) >= 0) { 88136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Subtract multiples of PI 882a8e45a3255089266ccf6500188a2d31af8634e4aHans Boehm BigInteger pi_multiples = scale(halfpi_multiples, -1); 8837b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm CR adjustment = PI.multiply(CR.valueOf(pi_multiples)); 884a8e45a3255089266ccf6500188a2d31af8634e4aHans Boehm if (pi_multiples.and(big1).signum() != 0) { 88536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return subtract(adjustment).cos().negate(); 88636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 88736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return subtract(adjustment).cos(); 88836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 8897b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm } else if (get_appr(-1).abs().compareTo(big2) >= 0) { 89036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Scale further with double angle formula 89136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR cos_half = shiftRight(1).cos(); 892b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return cos_half.multiply(cos_half).shiftLeft(1).subtract(ONE); 89336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 89436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new prescaled_cos_CR(this); 89536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 896349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 897349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 898349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 899349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The trigonometric sine function. 900349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 901349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR sin() { 90236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return half_pi.subtract(this).cos(); 903349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 904349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 905b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm/** 906b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm* The trignonometric arc (inverse) sine function. 907b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm*/ 908b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm public CR asin() { 909b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger rough_appr = get_appr(-10); 910b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm if (rough_appr.compareTo(big750) /* 1/sqrt(2) + a bit */ > 0){ 911b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm CR new_arg = ONE.subtract(multiply(this)).sqrt(); 912b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return new_arg.acos(); 913b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } else if (rough_appr.compareTo(bigm750) < 0) { 914b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return negate().asin().negate(); 915b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } else { 916b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return new prescaled_asin_CR(this); 917b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 918b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 919b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 920b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm/** 921b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm* The trignonometric arc (inverse) cosine function. 922b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm*/ 923b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm public CR acos() { 924b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return half_pi.subtract(asin()); 925b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 926b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 927349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger low_ln_limit = big8; /* sixteenths, i.e. 1/2 */ 928349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static final BigInteger high_ln_limit = 92936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger.valueOf(16 + 8 /* 1.5 */); 93036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm static final BigInteger scaled_4 = 93136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger.valueOf(4*16); 932349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 933349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 934349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The natural (base e) logarithm. 935349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 936349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR ln() { 93736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm final int low_prec = -4; 93836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger rough_appr = get_appr(low_prec); /* In sixteenths */ 93936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.compareTo(big0) < 0) { 9407b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm throw new ArithmeticException("ln(negative)"); 94136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 94236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.compareTo(low_ln_limit) <= 0) { 94336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return inverse().ln().negate(); 94436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 94536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.compareTo(high_ln_limit) >= 0) { 94636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (rough_appr.compareTo(scaled_4) <= 0) { 94736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR quarter = sqrt().sqrt().ln(); 94836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return quarter.shiftLeft(2); 94936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 95036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int extra_bits = rough_appr.bitLength() - 3; 95136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR scaled_result = shiftRight(extra_bits).ln(); 95236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scaled_result.add(CR.valueOf(extra_bits).multiply(ln2)); 95336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 95436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 95536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return simple_ln(); 956349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 957349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 958349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm/** 959349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm* The square root of a constructive real. 960349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm*/ 961349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm public CR sqrt() { 96236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return new sqrt_CR(this); 963349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 964349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 965349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} // end of CR 966349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 967349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 968349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// 969349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// A specialization of CR for cases in which approximate() calls 970349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// to increase evaluation precision are somewhat expensive. 971349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// If we need to (re)evaluate, we speculatively evaluate to slightly 972349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// higher precision, miminimizing reevaluations. 973349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Note that this requires any arguments to be evaluated to higher 974349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// precision than absolutely necessary. It can thus potentially 975349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// result in lots of wasted effort, and should be used judiciously. 976349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// This assumes that the order of magnitude of the number is roughly one. 977349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// 978349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmabstract class slow_CR extends CR { 979349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static int max_prec = -64; 980349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm static int prec_incr = 32; 9817b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm public synchronized BigInteger get_appr(int precision) { 98236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm check_prec(precision); 98336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (appr_valid && precision >= min_prec) { 98436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(max_appr, min_prec - precision); 98536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 98636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int eval_prec = (precision >= max_prec? max_prec : 98736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm (precision - prec_incr + 1) & ~(prec_incr - 1)); 98836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger result = approximate(eval_prec); 98936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm min_prec = eval_prec; 99036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm max_appr = result; 99136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm appr_valid = true; 99236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(result, eval_prec - precision); 99336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 994349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 995349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 996349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 997349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 998349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of an integer constant. Private. 999349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass int_CR extends CR { 1000349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm BigInteger value; 1001349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int_CR(BigInteger n) { 100236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm value = n; 1003349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1004349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 100536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(value, -p) ; 1006349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1007349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1008349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 100997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Representation of a number that may not have been completely 101097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// evaluated, but is assumed to be an integer. Hence we never 101197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// evaluate beyond the decimal point. 101297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehmclass assumed_int_CR extends CR { 101397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm CR value; 101497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm assumed_int_CR(CR x) { 101597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm value = x; 101697767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 101797767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm protected BigInteger approximate(int p) { 101897767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm if (p >= 0) { 101997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return value.get_appr(p); 102097767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } else { 102197767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm return scale(value.get_appr(0), -p) ; 102297767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 102397767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm } 102497767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm} 102597767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm 1026349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of the sum of 2 constructive reals. Private. 1027349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass add_CR extends CR { 1028349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op1; 1029349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op2; 1030349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm add_CR(CR x, CR y) { 103136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op1 = x; 103236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op2 = y; 1033349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1034349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 103536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Args need to be evaluated so that each error is < 1/4 ulp. 103636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Rounding error from the cale call is <= 1/2 ulp, so that 103736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // final error is < 1 ulp. 103836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(op1.get_appr(p-2).add(op2.get_appr(p-2)), -2); 1039349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1040349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1041349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1042349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of a CR multiplied by 2**n 1043349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass shifted_CR extends CR { 1044349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1045349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int count; 1046349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm shifted_CR(CR x, int n) { 104736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op = x; 104836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm count = n; 1049349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1050349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 105136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return op.get_appr(p - count); 1052349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1053349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1054349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1055349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of the negation of a constructive real. Private. 1056349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass neg_CR extends CR { 1057349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1058349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm neg_CR(CR x) { 105936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op = x; 1060349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1061349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 106236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return op.get_appr(p).negate(); 1063349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1064349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1065349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1066349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of: 106736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// op1 if selector < 0 106836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// op2 if selector >= 0 106936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm// Assumes x = y if s = 0 1070349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass select_CR extends CR { 1071349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR selector; 1072349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int selector_sign; 1073349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op1; 1074349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op2; 1075349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm select_CR(CR s, CR x, CR y) { 107636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm selector = s; 107736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int selector_sign = selector.get_appr(-20).signum(); 107836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op1 = x; 107936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op2 = y; 1080349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1081349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 108236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (selector_sign < 0) return op1.get_appr(p); 108336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (selector_sign > 0) return op2.get_appr(p); 108436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op1_appr = op1.get_appr(p-1); 108536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op2_appr = op2.get_appr(p-1); 108636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger diff = op1_appr.subtract(op2_appr).abs(); 108736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (diff.compareTo(big1) <= 0) { 108836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // close enough; use either 108936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(op1_appr, -1); 109036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 109136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // op1 and op2 are different; selector != 0; 109236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // safe to get sign of selector. 109336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (selector.signum() < 0) { 109436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm selector_sign = -1; 109536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(op1_appr, -1); 109636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 109736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm selector_sign = 1; 109836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(op2_appr, -1); 109936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 1100349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1101349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1102349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1103349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of the product of 2 constructive reals. Private. 1104349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass mult_CR extends CR { 1105349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op1; 1106349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op2; 1107349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm mult_CR(CR x, CR y) { 110836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op1 = x; 110936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op2 = y; 1110349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1111349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 111236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int half_prec = (p >> 1) - 1; 111336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int msd_op1 = op1.msd(half_prec); 111436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int msd_op2; 111536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm 111636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (msd_op1 == Integer.MIN_VALUE) { 111736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm msd_op2 = op2.msd(half_prec); 111836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (msd_op2 == Integer.MIN_VALUE) { 111936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Product is small enough that zero will do as an 112036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // approximation. 112136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return big0; 112236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 112336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Swap them, so the larger operand (in absolute value) 112436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // is first. 112536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm CR tmp; 112636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm tmp = op1; 112736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op1 = op2; 112836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op2 = tmp; 112936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm msd_op1 = msd_op2; 113036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 113136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 113236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // msd_op1 is valid at this point. 113336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int prec2 = p - msd_op1 - 3; // Precision needed for op2. 113436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // The appr. error is multiplied by at most 113536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2 ** (msd_op1 + 1) 113636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Thus each approximation contributes 1/4 ulp 113736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // to the rounding error, and the final rounding adds 113836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // another 1/2 ulp. 113936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger appr2 = op2.get_appr(prec2); 1140349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm if (appr2.signum() == 0) return big0; 114136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm msd_op2 = op2.known_msd(); 114236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int prec1 = p - msd_op2 - 3; // Precision needed for op1. 114336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger appr1 = op1.get_appr(prec1); 114436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int scale_digits = prec1 + prec2 - p; 114536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(appr1.multiply(appr2), scale_digits); 1146349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1147349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1148349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 114997767fab3e7e9a458613a4deb580e73695b09e7eHans Boehm// Representation of the multiplicative inverse of a constructive 1150349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// real. Private. Should use Newton iteration to refine estimates. 1151349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass inv_CR extends CR { 1152349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1153349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm inv_CR(CR x) { op = x; } 1154349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 115536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int msd = op.msd(); 115636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int inv_msd = 1 - msd; 115736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int digits_needed = inv_msd - p + 3; 115836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Number of SIGNIFICANT digits needed for 115936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // argument, excl. msd position, which may 116036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // be fictitious, since msd routine can be 116136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // off by 1. Roughly 1 extra digit is 116236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // needed since the relative error is the 116336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // same in the argument and result, but 116436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // this isn't quite the same as the number 116536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // of significant digits. Another digit 116636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // is needed to compensate for slop in the 1167349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // calculation. 116836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // One further bit is required, since the 116936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // final rounding introduces a 0.5 ulp 117036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // error. 117136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int prec_needed = msd - digits_needed; 117236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int log_scale_factor = -p - prec_needed; 117336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (log_scale_factor < 0) return big0; 117436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger dividend = big1.shiftLeft(log_scale_factor); 117536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_divisor = op.get_appr(prec_needed); 117636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger abs_scaled_divisor = scaled_divisor.abs(); 117736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger adj_dividend = dividend.add( 117836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm abs_scaled_divisor.shiftRight(1)); 117936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Adjustment so that final result is rounded. 118036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger result = adj_dividend.divide(abs_scaled_divisor); 118136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (scaled_divisor.signum() < 0) { 118236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result.negate(); 118336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 118436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return result; 118536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 1186349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1187349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1188349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1189349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1190349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of the exponential of a constructive real. Private. 1191349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Uses a Taylor series expansion. Assumes x < 1/2. 1192349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Note: this is known to be a bad algorithm for 1193349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// floating point. Unfortunately, other alternatives 1194349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// appear to require precomputed information. 1195349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass prescaled_exp_CR extends CR { 1196349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1197349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm prescaled_exp_CR(CR x) { op = x; } 1198349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 119936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (p >= 1) return big0; 120036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int iterations_needed = -p/2 + 2; // conservative estimate > 0. 120136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Claim: each intermediate term is accurate 120236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // to 2*2^calc_precision. 120336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Total rounding error in series computation is 120436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2*iterations_needed*2^calc_precision, 120536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // exclusive of error in op. 120636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int calc_precision = p - bound_log2(2*iterations_needed) 120736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm - 4; // for error in op, truncation. 120836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int op_prec = p - 3; 120936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op_appr = op.get_appr(op_prec); 121036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Error in argument results in error of < 3/8 ulp. 121136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Sum of term eval. rounding error is < 1/16 ulp. 121236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Series truncation error < 1/16 ulp. 121336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Final rounding error is <= 1/2 ulp. 121436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Thus final error is < 1 ulp. 121536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 121636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_term = scaled_1; 121736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_sum = scaled_1; 121836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int n = 0; 121936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger max_trunc_error = 122036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm big1.shiftLeft(p - 4 - calc_precision); 122136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (current_term.abs().compareTo(max_trunc_error) >= 0) { 12229666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 122336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm n += 1; 122436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm /* current_term = current_term * op / n */ 122536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = scale(current_term.multiply(op_appr), op_prec); 122636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = current_term.divide(BigInteger.valueOf(n)); 122736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_sum = current_sum.add(current_term); 122836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 122936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(current_sum, calc_precision - p); 1230349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1231349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1232349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1233349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation of the cosine of a constructive real. Private. 1234349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Uses a Taylor series expansion. Assumes |x| < 1. 1235349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass prescaled_cos_CR extends slow_CR { 1236349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1237349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm prescaled_cos_CR(CR x) { 123836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm op = x; 1239349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1240349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 124136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (p >= 1) return big0; 124236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int iterations_needed = -p/2 + 4; // conservative estimate > 0. 124336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Claim: each intermediate term is accurate 124436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // to 2*2^calc_precision. 124536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Total rounding error in series computation is 124636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2*iterations_needed*2^calc_precision, 124736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // exclusive of error in op. 124836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int calc_precision = p - bound_log2(2*iterations_needed) 124936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm - 4; // for error in op, truncation. 125036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int op_prec = p - 2; 125136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op_appr = op.get_appr(op_prec); 125236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Error in argument results in error of < 1/4 ulp. 125336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Cumulative arithmetic rounding error is < 1/16 ulp. 125436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Series truncation error < 1/16 ulp. 125536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Final rounding error is <= 1/2 ulp. 125636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Thus final error is < 1 ulp. 125736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_term; 125836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int n; 125936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger max_trunc_error = 126036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm big1.shiftLeft(p - 4 - calc_precision); 126136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm n = 0; 126236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = big1.shiftLeft(-calc_precision); 126336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_sum = current_term; 126436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (current_term.abs().compareTo(max_trunc_error) >= 0) { 12659666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 126636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm n += 2; 126736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm /* current_term = - current_term * op * op / n * (n - 1) */ 126836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = scale(current_term.multiply(op_appr), op_prec); 126936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = scale(current_term.multiply(op_appr), op_prec); 127036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger divisor = BigInteger.valueOf(-n) 127136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm .multiply(BigInteger.valueOf(n-1)); 127236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = current_term.divide(divisor); 127336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_sum = current_sum.add(current_term); 127436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 127536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(current_sum, calc_precision - p); 1276349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1277349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1278349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1279349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// The constructive real atan(1/n), where n is a small integer 1280349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// > base. 1281349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// This gives a simple and moderately fast way to compute PI. 1282349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass integral_atan_CR extends slow_CR { 1283349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm int op; 1284349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm integral_atan_CR(int x) { op = x; } 1285349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 128636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (p >= 1) return big0; 128736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int iterations_needed = -p/2 + 2; // conservative estimate > 0. 128836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Claim: each intermediate term is accurate 128936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // to 2*base^calc_precision. 129036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Total rounding error in series computation is 129136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2*iterations_needed*base^calc_precision, 129236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // exclusive of error in op. 129336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int calc_precision = p - bound_log2(2*iterations_needed) 129436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm - 2; // for error in op, truncation. 129536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Error in argument results in error of < 3/8 ulp. 129636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Cumulative arithmetic rounding error is < 1/4 ulp. 129736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Series truncation error < 1/4 ulp. 129836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Final rounding error is <= 1/2 ulp. 129936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Thus final error is < 1 ulp. 130036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 1301349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm BigInteger big_op = BigInteger.valueOf(op); 1302349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm BigInteger big_op_squared = BigInteger.valueOf(op*op); 130336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op_inverse = scaled_1.divide(big_op); 130436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_power = op_inverse; 130536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_term = op_inverse; 130636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_sum = op_inverse; 130736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int current_sign = 1; 130836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int n = 1; 130936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger max_trunc_error = 131036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm big1.shiftLeft(p - 2 - calc_precision); 131136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (current_term.abs().compareTo(max_trunc_error) >= 0) { 13129666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 131336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm n += 2; 131436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_power = current_power.divide(big_op_squared); 131536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_sign = -current_sign; 131636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = 131736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_power.divide(BigInteger.valueOf(current_sign*n)); 1318349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm current_sum = current_sum.add(current_term); 131936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 132036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(current_sum, calc_precision - p); 1321349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1322349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1323349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1324349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm// Representation for ln(1 + op) 1325349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass prescaled_ln_CR extends slow_CR { 1326349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1327349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm prescaled_ln_CR(CR x) { op = x; } 132836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Compute an approximation of ln(1+x) to precision 132936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // prec. This assumes |x| < 1/2. 133036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // It uses a Taylor series expansion. 133136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Unfortunately there appears to be no way to take 133236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // advantage of old information. 133336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Note: this is known to be a bad algorithm for 1334349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm // floating point. Unfortunately, other alternatives 133536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // appear to require precomputed tabular information. 1336349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 133736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (p >= 0) return big0; 133836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int iterations_needed = -p; // conservative estimate > 0. 133936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Claim: each intermediate term is accurate 134036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // to 2*2^calc_precision. Total error is 134136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // 2*iterations_needed*2^calc_precision 134236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // exclusive of error in op. 134336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int calc_precision = p - bound_log2(2*iterations_needed) 134436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm - 4; // for error in op, truncation. 134536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int op_prec = p - 3; 134636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op_appr = op.get_appr(op_prec); 134736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Error analysis as for exponential. 134836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_1 = big1.shiftLeft(-calc_precision); 134936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger x_nth = scale(op_appr, op_prec - calc_precision); 135036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_term = x_nth; // x**n 135136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger current_sum = current_term; 135236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int n = 1; 135336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int current_sign = 1; // (-1)^(n-1) 135436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger max_trunc_error = 135536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm big1.shiftLeft(p - 4 - calc_precision); 135636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm while (current_term.abs().compareTo(max_trunc_error) >= 0) { 13579666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 135836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm n += 1; 1359349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm current_sign = -current_sign; 136036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm x_nth = scale(x_nth.multiply(op_appr), op_prec); 136136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_term = x_nth.divide(BigInteger.valueOf(n * current_sign)); 136236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // x**n / (n * (-1)**(n-1)) 136336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm current_sum = current_sum.add(current_term); 136436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 136536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return scale(current_sum, calc_precision - p); 1366349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1367349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1368349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm 1369b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm// Representation of the arcsine of a constructive real. Private. 1370b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm// Uses a Taylor series expansion. Assumes |x| < (1/2)^(1/3). 1371b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehmclass prescaled_asin_CR extends slow_CR { 1372b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm CR op; 1373b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm prescaled_asin_CR(CR x) { 1374b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm op = x; 1375b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 1376b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm protected BigInteger approximate(int p) { 1377b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // The Taylor series is the sum of x^(2n+1) * (2n)!/(4^n n!^2 (2n+1)) 1378b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Note that (2n)!/(4^n n!^2) is always less than one. 1379b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // (The denominator is effectively 2n*2n*(2n-2)*(2n-2)*...*2*2 1380b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // which is clearly > (2n)!) 1381b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Thus all terms are bounded by x^(2n+1). 1382b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Unfortunately, there's no easy way to prescale the argument 1383b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // to less than 1/sqrt(2), and we can only approximate that. 1384b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Thus the worst case iteration count is fairly high. 1385b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // But it doesn't make much difference. 1386b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm if (p >= 2) return big0; // Never bigger than 4. 1387b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm int iterations_needed = -3 * p / 2 + 4; 1388b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // conservative estimate > 0. 1389b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Follows from assumed bound on x and 1390b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // the fact that only every other Taylor 1391b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Series term is present. 1392b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Claim: each intermediate term is accurate 1393b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // to 2*2^calc_precision. 1394b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Total rounding error in series computation is 1395b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // 2*iterations_needed*2^calc_precision, 1396b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // exclusive of error in op. 1397b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm int calc_precision = p - bound_log2(2*iterations_needed) 1398b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm - 4; // for error in op, truncation. 1399b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm int op_prec = p - 3; // always <= -2 1400b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger op_appr = op.get_appr(op_prec); 1401b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Error in argument results in error of < 1/4 ulp. 1402b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // (Derivative is bounded by 2 in the specified range and we use 1403b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // 3 extra digits.) 1404b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Ignoring the argument error, each term has an error of 1405b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // < 3ulps relative to calc_precision, which is more precise than p. 1406b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Cumulative arithmetic rounding error is < 3/16 ulp (relative to p). 1407b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Series truncation error < 2/16 ulp. (Each computed term 1408b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // is at most 2/3 of last one, so some of remaining series < 1409b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // 3/2 * current term.) 1410b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Final rounding error is <= 1/2 ulp. 1411b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Thus final error is < 1 ulp (relative to p). 1412b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger max_last_term = 1413b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm big1.shiftLeft(p - 4 - calc_precision); 1414b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm int exp = 1; // Current exponent, = 2n+1 in above expression 1415b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger current_term = op_appr.shiftLeft(op_prec - calc_precision); 1416b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger current_sum = current_term; 1417b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger current_factor = current_term; 1418b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Current scaled Taylor series term 1419b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // before division by the exponent. 1420b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Accurate to 3 ulp at calc_precision. 1421b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm while (current_term.abs().compareTo(max_last_term) >= 0) { 14229666c57ab247046b716e9c1e56f0dcc7d4a1545fHans Boehm if (Thread.interrupted() || please_stop) throw new AbortedException(); 1423b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm exp += 2; 1424b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // current_factor = current_factor * op * op * (exp-1) * (exp-2) / 1425b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // (exp-1) * (exp-1), with the two exp-1 factors cancelling, 1426b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // giving 1427b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // current_factor = current_factor * op * op * (exp-2) / (exp-1) 1428b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Thus the error any in the previous term is multiplied by 1429b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // op^2, adding an error of < (1/2)^(2/3) < 2/3 the original 1430b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // error. 1431b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_factor = current_factor.multiply(BigInteger.valueOf(exp - 2)); 1432b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_factor = scale(current_factor.multiply(op_appr), op_prec + 2); 1433b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Carry 2 extra bits of precision forward; thus 1434b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // this effectively introduces 1/8 ulp error. 1435b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_factor = current_factor.multiply(op_appr); 1436b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm BigInteger divisor = BigInteger.valueOf(exp - 1); 1437b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_factor = current_factor.divide(divisor); 1438b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Another 1/4 ulp error here. 1439b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_factor = scale(current_factor, op_prec - 2); 1440b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Remove extra 2 bits. 1/2 ulp rounding error. 1441b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Current_factor has original 3 ulp rounding error, which we 1442b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // reduced by 1, plus < 1 ulp new rounding error. 1443b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_term = current_factor.divide(BigInteger.valueOf(exp)); 1444b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // Contributes 1 ulp error to sum plus at most 3 ulp 1445b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm // from current_factor. 1446b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm current_sum = current_sum.add(current_term); 1447b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 1448b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm return scale(current_sum, calc_precision - p); 1449b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 1450b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm } 1451b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 1452b849a8eb275670d2bd9b2dc535231cf333414f5eHans Boehm 1453349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehmclass sqrt_CR extends CR { 1454349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm CR op; 1455349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm sqrt_CR(CR x) { op = x; } 145636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm final int fp_prec = 50; // Conservative estimate of number of 145736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // significant bits in double precision 145836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // computation. 1459349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm final int fp_op_prec = 60; 1460349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm protected BigInteger approximate(int p) { 146136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int max_prec_needed = 2*p - 1; 146236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int msd = op.msd(max_prec_needed); 146336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (msd <= max_prec_needed) return big0; 146436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int result_msd = msd/2; // +- 1 146536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int result_digits = result_msd - p; // +- 2 146636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm if (result_digits > fp_prec) { 146736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Compute less precise approximation and use a Newton iter. 146836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int appr_digits = result_digits/2 + 6; 146936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // This should be conservative. Is fewer enough? 147036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int appr_prec = result_msd - appr_digits; 147136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger last_appr = get_appr(appr_prec); 147236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int prod_prec = 2*appr_prec; 147336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger op_appr = op.get_appr(prod_prec); 147436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Slightly fewer might be enough; 147536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Compute (last_appr * last_appr + op_appr)/(last_appr/2) 147636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // while adjusting the scaling to make everything work 147736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger prod_prec_scaled_numerator = 147836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm last_appr.multiply(last_appr).add(op_appr); 147936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_numerator = 148036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm scale(prod_prec_scaled_numerator, appr_prec - p); 148136dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger shifted_result = scaled_numerator.divide(last_appr); 148236dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return shifted_result.add(big1).shiftRight(1); 148336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } else { 148436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Use a double precision floating point approximation. 148536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm // Make sure all precisions are even 148636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int op_prec = (msd - fp_op_prec) & ~1; 148736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int working_prec = op_prec - fp_op_prec; 148836dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_bi_appr = op.get_appr(op_prec) 148936dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm .shiftLeft(fp_op_prec); 149036dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm double scaled_appr = scaled_bi_appr.doubleValue(); 14917b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm if (scaled_appr < 0.0) 14927b2383f4d4d94a6dd29a1f1ac581abf852875fcaHans Boehm throw new ArithmeticException("sqrt(negative)"); 149336dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm double scaled_fp_sqrt = Math.sqrt(scaled_appr); 149436dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm BigInteger scaled_sqrt = BigInteger.valueOf((long)scaled_fp_sqrt); 149536dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm int shift_count = working_prec/2 - p; 149636dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm return shift(scaled_sqrt, shift_count); 149736dc537535bd06a87d652c9d3c369e1f748aa45bHans Boehm } 1498349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm } 1499349dbd732e48cf93b11291291e056fc5a3829f31Hans Boehm} 1500