1995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm/* 2995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Copyright (C) 2016 The Android Open Source Project 3995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * 4995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Licensed under the Apache License, Version 2.0 (the "License"); 5995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * you may not use this file except in compliance with the License. 6995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * You may obtain a copy of the License at 7995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * 8995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * http://www.apache.org/licenses/LICENSE-2.0 9995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * 10995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Unless required by applicable law or agreed to in writing, software 11995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * distributed under the License is distributed on an "AS IS" BASIS, 12995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * See the License for the specific language governing permissions and 14995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * limitations under the License. 15995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 16995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 17995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehmpackage com.android.calculator2; 18995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 19995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehmimport java.math.BigInteger; 20995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehmimport com.hp.creals.CR; 21995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehmimport com.hp.creals.UnaryCRFunction; 22995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 23995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm/** 24995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Computable real numbers, represented so that we can get exact decidable comparisons 25995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * for a number of interesting special cases, including rational computations. 26995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * 27995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * A real number is represented as the product of two numbers with different representations: 28995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * A) A BoundedRational that can only represent a subset of the rationals, but supports 29995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * exact computable comparisons. 30995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * B) A lazily evaluated "constructive real number" that provides operations to evaluate 31995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * itself to any requested number of digits. 32995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Whenever possible, we choose (B) to be one of a small set of known constants about which we 33995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * know more. For example, whenever we can, we represent rationals such that (B) is 1. 34995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * This scheme allows us to do some very limited symbolic computation on numbers when both 35995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * have the same (B) value, as well as in some other situations. We try to maximize that 36995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * possibility. 37995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * 38995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Arithmetic operations and operations that produce finite approximations may throw unchecked 39995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * exceptions produced by the underlying CR and BoundedRational packages, including 40995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * CR.PrecisionOverflowException and CR.AbortedException. 41995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 42995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehmpublic class UnifiedReal { 43995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 44995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final BoundedRational mRatFactor; 45995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final CR mCrFactor; 46995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // TODO: It would be helpful to add flags to indicate whether the result is known 47995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // irrational, etc. This sometimes happens even if mCrFactor is not one of the known ones. 48995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // And exact comparisons between rationals and known irrationals are decidable. 49995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 50995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 51995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Perform some nontrivial consistency checks. 52995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * @hide 53995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 54995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static boolean enableChecks = true; 55995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 56995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static void check(boolean b) { 57995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (!b) { 58995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new AssertionError(); 59995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 60995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 61995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 62995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private UnifiedReal(BoundedRational rat, CR cr) { 63995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (rat == null) { 64995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("Building UnifiedReal from null"); 65995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 66995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // We don't normally traffic in null CRs, and hence don't test explicitly. 67995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm mCrFactor = cr; 68995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm mRatFactor = rat; 69995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 70995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 71995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal(CR cr) { 72995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm this(BoundedRational.ONE, cr); 73995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 74995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 75995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal(BoundedRational rat) { 76995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm this(rat, CR_ONE); 77995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 78995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 79995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal(BigInteger n) { 80995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm this(new BoundedRational(n)); 81995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 82995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 83995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal(long n) { 84995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm this(new BoundedRational(n)); 85995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 86995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 87995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Various helpful constants 88995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static BigInteger BIG_24 = BigInteger.valueOf(24); 89995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static int DEFAULT_COMPARE_TOLERANCE = -1000; 90995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 91995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Well-known CR constants we try to use in the mCrFactor position: 92995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_ONE = CR.ONE; 93995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_PI = CR.PI; 94995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_E = CR.ONE.exp(); 95995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_SQRT2 = CR.valueOf(2).sqrt(); 96995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_SQRT3 = CR.valueOf(3).sqrt(); 97995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN2 = CR.valueOf(2).ln(); 98995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN3 = CR.valueOf(3).ln(); 99995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN5 = CR.valueOf(5).ln(); 100995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN6 = CR.valueOf(6).ln(); 101995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN7 = CR.valueOf(7).ln(); 102995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR CR_LN10 = CR.valueOf(10).ln(); 103995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 104995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Square roots that we try to recognize. 105995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // We currently recognize only a small fixed collection, since the sqrt() function needs to 106995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // identify numbers of the form <SQRT[i]>*n^2, and we don't otherwise know of a good 107995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // algorithm for that. 108995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR sSqrts[] = { 109995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 110995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR.ONE, 111995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_SQRT2, 112995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_SQRT3, 113995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 114995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR.valueOf(5).sqrt(), 115995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR.valueOf(6).sqrt(), 116995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR.valueOf(7).sqrt(), 117995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 118995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 119995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR.valueOf(10).sqrt() }; 120995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 121995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Natural logs of small integers that we try to recognize. 122995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static CR sLogs[] = { 123995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 124995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 125995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN2, 126995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN3, 127995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 128995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN5, 129995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN6, 130995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN7, 131995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 132995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm null, 133995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR_LN10 }; 134995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 135995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 136995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Some convenient UnifiedReal constants. 137995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal PI = new UnifiedReal(CR_PI); 138995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal E = new UnifiedReal(CR_E); 139995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal ZERO = new UnifiedReal(BoundedRational.ZERO); 140995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal ONE = new UnifiedReal(BoundedRational.ONE); 141995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal MINUS_ONE = new UnifiedReal(BoundedRational.MINUS_ONE); 142995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal TWO = new UnifiedReal(BoundedRational.TWO); 143995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal MINUS_TWO = new UnifiedReal(BoundedRational.MINUS_TWO); 144995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal HALF = new UnifiedReal(BoundedRational.HALF); 145995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal MINUS_HALF = new UnifiedReal(BoundedRational.MINUS_HALF); 146995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal TEN = new UnifiedReal(BoundedRational.TEN); 147995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static final UnifiedReal RADIANS_PER_DEGREE 148995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm = new UnifiedReal(new BoundedRational(1, 180), CR_PI); 149995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal SIX = new UnifiedReal(6); 150995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal HALF_SQRT2 = new UnifiedReal(BoundedRational.HALF, CR_SQRT2); 151995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal SQRT3 = new UnifiedReal(CR_SQRT3); 152995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal HALF_SQRT3 = new UnifiedReal(BoundedRational.HALF, CR_SQRT3); 153995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal THIRD_SQRT3 = new UnifiedReal(BoundedRational.THIRD, CR_SQRT3); 154995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal PI_OVER_2 = new UnifiedReal(BoundedRational.HALF, CR_PI); 155995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal PI_OVER_3 = new UnifiedReal(BoundedRational.THIRD, CR_PI); 156995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal PI_OVER_4 = new UnifiedReal(BoundedRational.QUARTER, CR_PI); 157995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final UnifiedReal PI_OVER_6 = new UnifiedReal(BoundedRational.SIXTH, CR_PI); 158995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 159995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 160995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 161995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Given a constructive real cr, try to determine whether cr is the square root of 162995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * a small integer. If so, return its square as a BoundedRational. Otherwise return null. 163995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * We make this determination by simple table lookup, so spurious null returns are 164995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * entirely possible, or even likely. 165995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 166995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static BoundedRational getSquare(CR cr) { 167995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int i = 0; i < sSqrts.length; ++i) { 168995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sSqrts[i] == cr) { 169995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new BoundedRational(i); 170995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 171995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 172995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 173995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 174995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 175995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 176995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Given a constructive real cr, try to determine whether cr is the square root of 177995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * a small integer. If so, return its square as a BoundedRational. Otherwise return null. 178995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * We make this determination by simple table lookup, so spurious null returns are 179995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * entirely possible, or even likely. 180995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 181995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private BoundedRational getExp(CR cr) { 182995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int i = 0; i < sLogs.length; ++i) { 183995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sLogs[i] == cr) { 184995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new BoundedRational(i); 185995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 186995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 187995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 188995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 189995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 190995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 191995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * If the argument is a well-known constructive real, return its name. 192995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * The name of "CR_ONE" is the empty string. 193995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * No named constructive reals are rational multiples of each other. 194995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Thus two UnifiedReals with different named mCrFactors can be equal only if both 195995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * mRatFactors are zero or possibly if one is CR_PI and the other is CR_E. 196995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * (The latter is apparently an open problem.) 197995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 198995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static String crName(CR cr) { 199995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == CR_ONE) { 200995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ""; 201995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 202995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == CR_PI) { 203995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return "\u03C0"; // GREEK SMALL LETTER PI 204995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 205995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == CR_E) { 206995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return "e"; 207995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 208995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int i = 0; i < sSqrts.length; ++i) { 209995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == sSqrts[i]) { 210995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return "\u221A" /* SQUARE ROOT */ + i; 211995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 212995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 213995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int i = 0; i < sLogs.length; ++i) { 214995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == sLogs[i]) { 215995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return "ln(" + i + ")"; 216995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 217995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 218995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 219995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 220995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 221995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 222995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Would crName() return non-Null? 223995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 224995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static boolean isNamed(CR cr) { 225995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == CR_ONE || cr == CR_PI || cr == CR_E) { 226995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return true; 227995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 228995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (CR r: sSqrts) { 229995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == r) { 230995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return true; 231995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 232995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 233995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (CR r: sLogs) { 234995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (cr == r) { 235995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return true; 236995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 237995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 238995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return false; 239995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 240995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 241995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 242995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is cr known to be algebraic (as opposed to transcendental)? 243995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Currently only produces meaningful results for the above known special 244995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * constructive reals. 245995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 246995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static boolean definitelyAlgebraic(CR cr) { 247995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return cr == CR_ONE || getSquare(cr) != null; 248995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 249995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 250995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 251995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is this number known to be rational? 252995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 253995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyRational() { 254995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mCrFactor == CR_ONE || mRatFactor.signum() == 0; 255995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 256995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 257995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 258995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is this number known to be irrational? 259995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * TODO: We could track the fact that something is irrational with an explicit flag, which 260995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * could cover many more cases. Whether that matters in practice is TBD. 261995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 262995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyIrrational() { 263995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return !definitelyRational() && isNamed(mCrFactor); 264995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 265995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 266995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 267995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is this number known to be algebraic? 268995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 269995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyAlgebraic() { 270995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return definitelyAlgebraic(mCrFactor) || mRatFactor.signum() == 0; 271995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 272995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 273995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 274995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is this number known to be transcendental? 275995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 276995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyTranscendental() { 277995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return !definitelyAlgebraic() && isNamed(mCrFactor); 278995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 279995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 280995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 281995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 282995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is it known that the two constructive reals differ by something other than a 283995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * a rational factor, i.e. is it known that two UnifiedReals 284995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * with those mCrFactors will compare unequal unless both mRatFactors are zero? 285995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * If this returns true, then a comparison of two UnifiedReals using those two 286995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * mCrFactors cannot diverge, though we don't know of a good runtime bound. 287995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 288995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static boolean definitelyIndependent(CR r1, CR r2) { 289995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // The question here is whether r1 = x*r2, where x is rational, where r1 and r2 290995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // are in our set of special known CRs, can have a solution. 291995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // This cannot happen if one is CR_ONE and the other is not. 292995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // (Since all others are irrational.) 293995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // This cannot happen for two named square roots, which have no repeated factors. 294995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // (To see this, square both sides of the equation and factor. Each prime 295995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // factor in the numerator and denominator occurs twice.) 296995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // This cannot happen for e or pi on one side, and a square root on the other. 297995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // (One is transcendental, the other is algebraic.) 298995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // This cannot happen for two of our special natural logs. 299995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // (Otherwise ln(m) = (a/b)ln(n) ==> m = n^(a/b) ==> m^b = n^a, which is impossible 300995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // because either m or n includes a prime factor not shared by the other.) 301995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // This cannot happen for a log and a square root. 302995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // (The Lindemann-Weierstrass theorem tells us, among other things, that if 303995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // a is algebraic, then exp(a) is transcendental. Thus if l in our finite 304995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // set of logs where algebraic, expl(l), must be transacendental. 305995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // But exp(l) is an integer. Thus the logs are transcendental. But of course the 306995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // square roots are algebraic. Thus they can't be rational multiples.) 307995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Unfortunately, we do not know whether e/pi is rational. 308995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (r1 == r2) { 309995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return false; 310995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 311995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm CR other; 312995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (r1 == CR_E || r1 == CR_PI) { 313995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return definitelyAlgebraic(r2); 314995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 315995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (r2 == CR_E || r2 == CR_PI) { 316995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return definitelyAlgebraic(r1); 317995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 318995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return isNamed(r1) && isNamed(r2); 319995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 320995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 321995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 322995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Convert to String reflecting raw representation. 323995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Debug or log messages only, not pretty. 324995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 325995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public String toString() { 326995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.toString() + "*" + mCrFactor.toString(); 327995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 328995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 329995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 330995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Convert to readable String. 331995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Intended for user output. Produces exact expression when possible. 332995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 333995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public String toNiceString() { 334995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) { 335995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.toNiceString(); 336995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 337995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm String name = crName(mCrFactor); 338995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (name != null) { 339995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger bi = BoundedRational.asBigInteger(mRatFactor); 340995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (bi != null) { 341995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (bi.equals(BigInteger.ONE)) { 342995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return name; 343995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 344995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.toNiceString() + name; 345995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 346995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return "(" + mRatFactor.toNiceString() + ")" + name; 347995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 348995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mRatFactor.equals(BoundedRational.ONE)) { 349995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mCrFactor.toString(); 350995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 351995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().toString(); 352995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 353995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 354995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 355995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Will toNiceString() produce an exact representation? 356995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 357995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean exactlyDisplayable() { 358995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crName(mCrFactor) != null; 359995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 360995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 361995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Number of extra bits used in evaluation below to prefer truncation to rounding. 362995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Must be <= 30. 363995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private final static int EXTRA_PREC = 10; 364995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 365995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /* 366995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Returns a truncated representation of the result. 367995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * If exactlyTruncatable(), we round correctly towards zero. Otherwise the resulting digit 368995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * string may occasionally be rounded up instead. 3693465c782410b9e19f9a4e92cb8c74647e38a5bdcHans Boehm * Always includes a decimal point in the result. 370995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * The result includes n digits to the right of the decimal point. 371995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * @param n result precision, >= 0 372995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 373995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public String toStringTruncated(int n) { 374995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO) { 375995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.toStringTruncated(n); 376995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 377995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final CR scaled = CR.valueOf(BigInteger.TEN.pow(n)).multiply(crValue()); 378995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm boolean negative = false; 379995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger intScaled; 380995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (exactlyTruncatable()) { 381995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = scaled.get_appr(0); 382995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (intScaled.signum() < 0) { 383995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm negative = true; 384995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = intScaled.negate(); 385995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 386995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (CR.valueOf(intScaled).compareTo(scaled.abs()) > 0) { 387995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = intScaled.subtract(BigInteger.ONE); 388995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 389995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm check(CR.valueOf(intScaled).compareTo(scaled.abs()) < 0); 390995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 391995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Approximate case. Exact comparisons are impossible. 392995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = scaled.get_appr(-EXTRA_PREC); 393995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (intScaled.signum() < 0) { 394995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm negative = true; 395995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = intScaled.negate(); 396995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 397995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm intScaled = intScaled.shiftRight(EXTRA_PREC); 398995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 399995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm String digits = intScaled.toString(); 400995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int len = digits.length(); 401995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (len < n + 1) { 40224c91edb89b3bfecb5167d3ba4f76cf0203cb547Hans Boehm digits = StringUtils.repeat('0', n + 1 - len) + digits; 403995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm len = n + 1; 404995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 405995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return (negative ? "-" : "") + digits.substring(0, len - n) + "." 406995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm + digits.substring(len - n); 407995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 408995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 409995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /* 410995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Can we compute correctly truncated approximations of this number? 411995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 412995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean exactlyTruncatable() { 413995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // If the value is known rational, we can do exact comparisons. 414995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // If the value is known irrational, then we can safely compare to rational approximations; 415995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // equality is impossible; hence the comparison must converge. 416995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // The only problem cases are the ones in which we don't know. 417995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mCrFactor == CR_ONE || mRatFactor == BoundedRational.ZERO || definitelyIrrational(); 418995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 419995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 420995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 421995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return a double approximation. 422995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * TODO: Result is correctly rounded if known to be rational. 423995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 424995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public double doubleValue() { 425995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE) { 426995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.doubleValue(); // Hopefully correctly rounded 427995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 428995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().doubleValue(); // Approximately correctly rounded 429995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 430995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 431995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 432995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public CR crValue() { 433995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.crValue().multiply(mCrFactor); 434995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 435995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 436995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 437995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Are this and r exactly comparable? 438995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 439995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean isComparable(UnifiedReal u) { 440995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // We check for ONE only to speed up the common case. 441995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // The use of a tolerance here means we can spuriously return false, not true. 442995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mCrFactor == u.mCrFactor 443995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm && (isNamed(mCrFactor) || mCrFactor.signum(DEFAULT_COMPARE_TOLERANCE) != 0) 444995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm || mRatFactor.signum() == 0 && u.mRatFactor.signum() == 0 445995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm || definitelyIndependent(mCrFactor, u.mCrFactor) 446995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm || crValue().compareTo(u.crValue(), DEFAULT_COMPARE_TOLERANCE) != 0; 447995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 448995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 449995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 450995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return +1 if this is greater than r, -1 if this is less than r, or 0 of the two are 451995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * known to be equal. 452995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * May diverge if the two are equal and !isComparable(r). 453995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 454995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int compareTo(UnifiedReal u) { 455995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyZero() && u.definitelyZero()) return 0; 456995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == u.mCrFactor) { 457995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int signum = mCrFactor.signum(); // Can diverge if mCRFactor == 0. 458995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return signum * mRatFactor.compareTo(u.mRatFactor); 459995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 460995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().compareTo(u.crValue()); // Can also diverge. 461995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 462995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 463995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 464995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return +1 if this is greater than r, -1 if this is less than r, or possibly 0 of the two are 465995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * within 2^a of each other. 466995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 467995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int compareTo(UnifiedReal u, int a) { 468995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isComparable(u)) { 469995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return compareTo(u); 470995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 471995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().compareTo(u.crValue(), a); 472995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 473995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 474995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 475995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 476995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return compareTo(ZERO, a). 477995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 478995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int signum(int a) { 479995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return compareTo(ZERO, a); 480995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 481995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 482995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 483995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return compareTo(ZERO). 484995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * May diverge for ZERO argument if !isComparable(ZERO). 485995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 486995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int signum() { 487995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return compareTo(ZERO); 488995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 489995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 490995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 491995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Equality comparison. May erroneously return true if values differ by less than 2^a, 492995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * and !isComparable(u). 493995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 494995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean approxEquals(UnifiedReal u, int a) { 495995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isComparable(u)) { 496995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyIndependent(mCrFactor, u.mCrFactor) 497995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm && (mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0)) { 498995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // No need to actually evaluate, though we don't know which is larger. 499995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return false; 500995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 501995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return compareTo(u) == 0; 502995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 503995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 504995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().compareTo(u.crValue(), a) == 0; 505995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 506995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 507995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 508995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Returns true if values are definitely known to be equal, false in all other cases. 509995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 510995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyEquals(UnifiedReal u) { 511995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return isComparable(u) && compareTo(u) == 0; 512995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 513995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 514995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 515995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Returns true if values are definitely known not to be equal, false in all other cases. 5164452c78e94b522c74d63056b0814a00281be296aHans Boehm * Performs no approximate evaluation. 517995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 518995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyNotEquals(UnifiedReal u) { 519995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm boolean isNamed = isNamed(mCrFactor); 520995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm boolean uIsNamed = isNamed(u.mCrFactor); 521995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isNamed && uIsNamed) { 522995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyIndependent(mCrFactor, u.mCrFactor)) { 523995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.signum() != 0 || u.mRatFactor.signum() != 0; 524995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else if (mCrFactor == u.mCrFactor) { 525995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return !mRatFactor.equals(u.mRatFactor); 526995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 527995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return !mRatFactor.equals(u.mRatFactor); 528995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 529995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mRatFactor.signum() == 0) { 530995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return uIsNamed && u.mRatFactor.signum() != 0; 531995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 532995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (u.mRatFactor.signum() == 0) { 533995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return isNamed && mRatFactor.signum() != 0; 534995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 535995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return false; 536995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 537995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 538995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // And some slightly faster convenience functions for special cases: 539995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 540995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyZero() { 541995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.signum() == 0; 542995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 543995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 5444452c78e94b522c74d63056b0814a00281be296aHans Boehm /** 5454452c78e94b522c74d63056b0814a00281be296aHans Boehm * Can this number be determined to be definitely nonzero without performing approximate 5464452c78e94b522c74d63056b0814a00281be296aHans Boehm * evaluation? 5474452c78e94b522c74d63056b0814a00281be296aHans Boehm */ 548995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyNonZero() { 549995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return isNamed(mCrFactor) && mRatFactor.signum() != 0; 550995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 551995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 552995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean definitelyOne() { 553995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mCrFactor == CR_ONE && mRatFactor.equals(BoundedRational.ONE); 554995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 555995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 556995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 557995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return equivalent BoundedRational, if known to exist, null otherwise 558995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 559995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public BoundedRational boundedRationalValue() { 560995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) { 561995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor; 562995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 563995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 564995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 565995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 566995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 567995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Returns equivalent BigInteger result if it exists, null if not. 568995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 569995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public BigInteger bigIntegerValue() { 570995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BoundedRational r = boundedRationalValue(); 571995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return BoundedRational.asBigInteger(r); 572995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 573995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 574995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal add(UnifiedReal u) { 575995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == u.mCrFactor) { 576995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.add(mRatFactor, u.mRatFactor); 577995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 578995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, mCrFactor); 579995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 580995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 581995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyZero()) { 582995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Avoid creating new mCrFactor, even if they don't currently match. 583995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return u; 584995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 585995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (u.definitelyZero()) { 586995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return this; 587995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 588995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().add(u.crValue())); 589995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 590995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 591995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal negate() { 592995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.negate(mRatFactor), mCrFactor); 593995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 594995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 595995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal subtract(UnifiedReal u) { 596995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return add(u.negate()); 597995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 598995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 599995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal multiply(UnifiedReal u) { 600995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Preserve a preexisting mCrFactor when we can. 601995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE) { 602995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor); 603995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 604995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, u.mCrFactor); 605995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 606995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 607995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (u.mCrFactor == CR_ONE) { 608995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor); 609995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 610995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, mCrFactor); 611995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 612995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 613995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyZero() || u.definitelyZero()) { 614995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ZERO; 615995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 616995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == u.mCrFactor) { 617995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational square = getSquare(mCrFactor); 618995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (square != null) { 619995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.multiply( 620995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.multiply(square, mRatFactor), u.mRatFactor); 621995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 622995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor); 623995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 624995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 625995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 626995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Probably a bit cheaper to multiply component-wise. 627995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.multiply(mRatFactor, u.mRatFactor); 628995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 629995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, mCrFactor.multiply(u.mCrFactor)); 630995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 631995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().multiply(u.crValue())); 632995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 633995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 634995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static class ZeroDivisionException extends ArithmeticException { 635995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public ZeroDivisionException() { 636995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm super("Division by zero"); 637995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 638995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 639995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 640995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 641995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return the reciprocal. 642995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 643995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal inverse() { 644995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyZero()) { 645995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ZeroDivisionException(); 646995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 647995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational square = getSquare(mCrFactor); 648995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (square != null) { 649995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // 1/sqrt(n) = sqrt(n)/n 650995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.inverse( 651995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.multiply(mRatFactor, square)); 652995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 653995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, mCrFactor); 654995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 655995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 656995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.inverse(mRatFactor), mCrFactor.inverse()); 657995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 658995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 659995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal divide(UnifiedReal u) { 660995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == u.mCrFactor) { 661995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (u.definitelyZero()) { 662995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ZeroDivisionException(); 663995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 664995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.divide(mRatFactor, u.mRatFactor); 665995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 666995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, CR_ONE); 667995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 668995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 669995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return multiply(u.inverse()); 670995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 671995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 672995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal sqrt() { 673995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE) { 674995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational ratSqrt; 675995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Check for all arguments of the form <perfect rational square> * small_int, 676995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // where small_int has a known sqrt. This includes the small_int = 1 case. 677995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int divisor = 1; divisor < sSqrts.length; ++divisor) { 678995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sSqrts[divisor] != null) { 679995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm ratSqrt = BoundedRational.sqrt( 680995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.divide(mRatFactor, new BoundedRational(divisor))); 681995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (ratSqrt != null) { 682995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(ratSqrt, sSqrts[divisor]); 683995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 684995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 685995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 686995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 687995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().sqrt()); 688995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 689995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 690995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 691995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return (this mod 2pi)/(pi/6) as a BigInteger, or null if that isn't easily possible. 692995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 693995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private BigInteger getPiTwelfths() { 694995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyZero()) return BigInteger.ZERO; 695995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_PI) { 696995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger quotient = BoundedRational.asBigInteger( 697995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.multiply(mRatFactor, BoundedRational.TWELVE)); 698995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (quotient == null) { 699995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 700995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 701995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return quotient.mod(BIG_24); 702995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 703995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 704995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 705995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 706995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 707995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Computer the sin() for an integer multiple n of pi/12, if easily representable. 708995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * @param n value between 0 and 23 inclusive. 709995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 710995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static UnifiedReal sinPiTwelfths(int n) { 711995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n >= 12) { 712995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal negResult = sinPiTwelfths(n - 12); 713995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return negResult == null ? null : negResult.negate(); 714995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 715995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm switch (n) { 716995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 0: 717995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ZERO; 718995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 2: // 30 degrees 719995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF; 720995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 3: // 45 degrees 721995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF_SQRT2; 722995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 4: // 60 degrees 723995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF_SQRT3; 724995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 6: 725995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ONE; 726995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 8: 727995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF_SQRT3; 728995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 9: 729995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF_SQRT2; 730995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 10: 731995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return HALF; 732995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm default: 733995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return null; 734995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 735995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 736995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 737995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal sin() { 738995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger piTwelfths = getPiTwelfths(); 739995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (piTwelfths != null) { 740995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal result = sinPiTwelfths(piTwelfths.intValue()); 741995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (result != null) { 742995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return result; 743995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 744995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 745995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().sin()); 746995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 747995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 748995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static UnifiedReal cosPiTwelfths(int n) { 749995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int sinArg = n + 6; 750995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sinArg >= 24) { 751995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm sinArg -= 24; 752995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 753995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return sinPiTwelfths(sinArg); 754995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 755995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 756995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal cos() { 757995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger piTwelfths = getPiTwelfths(); 758995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (piTwelfths != null) { 759995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal result = cosPiTwelfths(piTwelfths.intValue()); 760995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (result != null) { 761995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return result; 762995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 763995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 764995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().cos()); 765995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 766995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 767995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal tan() { 768995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger piTwelfths = getPiTwelfths(); 769995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (piTwelfths != null) { 770995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int i = piTwelfths.intValue(); 771995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (i == 6 || i == 18) { 772995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("Tangent undefined"); 773995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 774995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal top = sinPiTwelfths(i); 775995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal bottom = cosPiTwelfths(i); 776995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (top != null && bottom != null) { 777995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return top.divide(bottom); 778995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 779995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 780995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return sin().divide(cos()); 781995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 782995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 783995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Throw an exception if the argument is definitely out of bounds for asin or acos. 784995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private void checkAsinDomain() { 785995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isComparable(ONE) && (compareTo(ONE) > 0 || compareTo(MINUS_ONE) < 0)) { 786995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("inverse trig argument out of range"); 787995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 788995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 789995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 790995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 791995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return asin(n/2). n is between -2 and 2. 792995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 793995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public static UnifiedReal asinHalves(int n){ 794995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n < 0) { 795995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return (asinHalves(-n).negate()); 796995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 797995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm switch (n) { 798995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 0: 799995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ZERO; 800995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 1: 801995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.SIXTH, CR.PI); 802995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 2: 803995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.HALF, CR.PI); 804995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 805995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new AssertionError("asinHalves: Bad argument"); 806995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 807995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 808995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 809995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return asin of this, assuming this is not an integral multiple of a half. 810995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 811995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal asinNonHalves() 812995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm { 813995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (compareTo(ZERO, -10) < 0) { 814995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return negate().asinNonHalves().negate(); 815995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 816995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(HALF_SQRT2)) { 817995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.QUARTER, CR_PI); 818995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 819995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(HALF_SQRT3)) { 820995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(BoundedRational.THIRD, CR_PI); 821995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 822995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().asin()); 823995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 824995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 825995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal asin() { 826995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm checkAsinDomain(); 827995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BigInteger halves = multiply(TWO).bigIntegerValue(); 828995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (halves != null) { 829995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return asinHalves(halves.intValue()); 830995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 831995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR.ONE || mCrFactor != CR_SQRT2 ||mCrFactor != CR_SQRT3) { 832995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return asinNonHalves(); 833995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 834995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().asin()); 835995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 836995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 837995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal acos() { 838995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return PI_OVER_2.subtract(asin()); 839995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 840995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 841995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal atan() { 842995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (compareTo(ZERO, -10) < 0) { 843995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return negate().atan().negate(); 844995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 845995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BigInteger asBI = bigIntegerValue(); 846995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (asBI != null && asBI.compareTo(BigInteger.ONE) <= 0) { 847995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final int asInt = asBI.intValue(); 848995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // These seem to be all rational cases: 849995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm switch (asInt) { 850995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 0: 851995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ZERO; 852995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm case 1: 853995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return PI_OVER_4; 854995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm default: 855995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new AssertionError("Impossible r_int"); 856995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 857995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 858995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(THIRD_SQRT3)) { 859995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return PI_OVER_6; 860995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 861995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(SQRT3)) { 862995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return PI_OVER_3; 863995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 864995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(UnaryCRFunction.atanFunction.execute(crValue())); 865995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 866995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 867995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static final BigInteger BIG_TWO = BigInteger.valueOf(2); 868995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 869995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 8704452c78e94b522c74d63056b0814a00281be296aHans Boehm * Compute an integral power of a constrive real, using the standard recursive algorithm. 8714452c78e94b522c74d63056b0814a00281be296aHans Boehm * exp is known to be positive. 8724452c78e94b522c74d63056b0814a00281be296aHans Boehm */ 8734452c78e94b522c74d63056b0814a00281be296aHans Boehm private static CR recursivePow(CR base, BigInteger exp) { 8744452c78e94b522c74d63056b0814a00281be296aHans Boehm if (exp.equals(BigInteger.ONE)) { 8754452c78e94b522c74d63056b0814a00281be296aHans Boehm return base; 8764452c78e94b522c74d63056b0814a00281be296aHans Boehm } 8774452c78e94b522c74d63056b0814a00281be296aHans Boehm if (exp.and(BigInteger.ONE).intValue() == 1) { 8784452c78e94b522c74d63056b0814a00281be296aHans Boehm return base.multiply(recursivePow(base, exp.subtract(BigInteger.ONE))); 8794452c78e94b522c74d63056b0814a00281be296aHans Boehm } 8804452c78e94b522c74d63056b0814a00281be296aHans Boehm CR tmp = recursivePow(base, exp.shiftRight(1)); 8814452c78e94b522c74d63056b0814a00281be296aHans Boehm if (Thread.interrupted()) { 8824452c78e94b522c74d63056b0814a00281be296aHans Boehm throw new CR.AbortedException(); 8834452c78e94b522c74d63056b0814a00281be296aHans Boehm } 8844452c78e94b522c74d63056b0814a00281be296aHans Boehm return tmp.multiply(tmp); 8854452c78e94b522c74d63056b0814a00281be296aHans Boehm } 8864452c78e94b522c74d63056b0814a00281be296aHans Boehm 8874452c78e94b522c74d63056b0814a00281be296aHans Boehm /** 888995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Compute an integral power of this. 8894452c78e94b522c74d63056b0814a00281be296aHans Boehm * This recurses roughly as deeply as the number of bits in the exponent, and can, in 8904452c78e94b522c74d63056b0814a00281be296aHans Boehm * ridiculous cases, result in a stack overflow. 891995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 892995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private UnifiedReal pow(BigInteger exp) { 893995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (exp.signum() < 0) { 894995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return pow(exp.negate()).inverse(); 895995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 896995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (exp.equals(BigInteger.ONE)) { 897995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return this; 898995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 899995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (exp.signum() == 0) { 900995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Questionable if base has undefined value. Java.lang.Math.pow() returns 1 anyway, 901995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // so we do the same. 902995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ONE; 903995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 904995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE) { 905995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BoundedRational ratPow = mRatFactor.pow(exp); 906995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (ratPow != null) { 907995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(mRatFactor.pow(exp)); 908995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 909995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 910995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational square = getSquare(mCrFactor); 911995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (square != null) { 912995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BoundedRational nRatFactor = 913995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.multiply(mRatFactor.pow(exp), square.pow(exp.shiftRight(1))); 914995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 915995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (exp.and(BigInteger.ONE).intValue() == 1) { 916995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Odd power: Multiply by remaining square root. 917995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, mCrFactor); 918995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 919995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor); 920995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 921995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 922995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 9234452c78e94b522c74d63056b0814a00281be296aHans Boehm if (signum(DEFAULT_COMPARE_TOLERANCE) > 0) { 9244452c78e94b522c74d63056b0814a00281be296aHans Boehm // Safe to take the log. This avoids deep recursion for huge exponents, which 9254452c78e94b522c74d63056b0814a00281be296aHans Boehm // may actually make sense here. 9264452c78e94b522c74d63056b0814a00281be296aHans Boehm return new UnifiedReal(crValue().ln().multiply(CR.valueOf(exp)).exp()); 9274452c78e94b522c74d63056b0814a00281be296aHans Boehm } else { 9284452c78e94b522c74d63056b0814a00281be296aHans Boehm // Possibly negative base with integer exponent. Use a recursive computation. 9294452c78e94b522c74d63056b0814a00281be296aHans Boehm // (Another possible option would be to use the absolute value of the base, and then 9304452c78e94b522c74d63056b0814a00281be296aHans Boehm // adjust the sign at the end. But that would have to be done in the CR 9314452c78e94b522c74d63056b0814a00281be296aHans Boehm // implementation.) 9324452c78e94b522c74d63056b0814a00281be296aHans Boehm return new UnifiedReal(recursivePow(crValue(), exp)); 9334452c78e94b522c74d63056b0814a00281be296aHans Boehm } 934995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 935995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 936995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal pow(UnifiedReal expon) { 937995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_E) { 938995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mRatFactor.equals(BoundedRational.ONE)) { 939995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return expon.exp(); 940995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 941995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal ratPart = new UnifiedReal(mRatFactor).pow(expon); 942995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return expon.exp().multiply(ratPart); 943995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 944995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 945995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BoundedRational expAsBR = expon.boundedRationalValue(); 946995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (expAsBR != null) { 947995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger expAsBI = BoundedRational.asBigInteger(expAsBR); 948995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (expAsBI != null) { 949995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return pow(expAsBI); 950995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 951995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Check for exponent that is a multiple of a half. 952995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm expAsBI = BoundedRational.asBigInteger( 953995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.multiply(BoundedRational.TWO, expAsBR)); 954995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (expAsBI != null) { 955995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return pow(expAsBI).sqrt(); 956995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 957995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 958995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 959995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().ln().multiply(expon.crValue()).exp()); 960995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 961995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 962995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 963995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Raise the argument to the 16th power. 964995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 965995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static long pow16(int n) { 966995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n > 10) { 967995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new AssertionError("Unexpexted pow16 argument"); 968995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 969995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm long result = n*n; 970995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm result *= result; 971995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm result *= result; 972995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm result *= result; 973995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return result; 974995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 975995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 976995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 977995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return the integral log with respect to the given base if it exists, 0 otherwise. 978995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * n is presumed positive. 979995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 980995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static long getIntLog(BigInteger n, int base) { 981995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm double nAsDouble = n.doubleValue(); 982995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm double approx = Math.log(nAsDouble)/Math.log(base); 983995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // A relatively quick test first. 984995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Unfortunately, this doesn't help for values to big to fit in a Double. 985995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (!Double.isInfinite(nAsDouble) && Math.abs(approx - Math.rint(approx)) > 1.0e-6) { 986995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return 0; 987995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 988995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm long result = 0; 989995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger remaining = n; 990995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger bigBase = BigInteger.valueOf(base); 991995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger base16th = null; // base^16, computed lazily 992995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm while (n.mod(bigBase).signum() == 0) { 993995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (Thread.interrupted()) { 994995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new CR.AbortedException(); 995995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 996995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm n = n.divide(bigBase); 997995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm ++result; 998995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // And try a slightly faster computation for large n: 999995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (base16th == null) { 1000995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm base16th = BigInteger.valueOf(pow16(base)); 1001995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1002995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm while (n.mod(base16th).signum() == 0) { 1003995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm n = n.divide(base16th); 1004995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm result += 16; 1005995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1006995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1007995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n.equals(BigInteger.ONE)) { 1008995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return result; 1009995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1010995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return 0; 1011995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1012995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1013995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal ln() { 1014995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isComparable(ZERO)) { 1015995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (signum() <= 0) { 1016995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("log(non-positive)"); 1017995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1018995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int compare1 = compareTo(ONE, DEFAULT_COMPARE_TOLERANCE); 1019995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (compare1 == 0) { 1020995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(ONE)) { 1021995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ZERO; 1022995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1023995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else if (compare1 < 0) { 1024995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return inverse().ln().negate(); 1025995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1026995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BigInteger bi = BoundedRational.asBigInteger(mRatFactor); 1027995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (bi != null) { 1028995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE) { 1029995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Check for a power of a small integer. We can use sLogs[] to return 1030995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // a more useful answer for those. 1031995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (int i = 0; i < sLogs.length; ++i) { 1032995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sLogs[i] != null) { 1033995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm long intLog = getIntLog(bi, i); 1034995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (intLog != 0) { 1035995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(new BoundedRational(intLog), sLogs[i]); 1036995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1037995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1038995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1039995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 1040995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Check for n^k * sqrt(n), for which we can also return a more useful answer. 1041995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational square = getSquare(mCrFactor); 1042995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (square != null) { 1043995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm int intSquare = square.intValue(); 1044995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (sLogs[intSquare] != null) { 1045995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm long intLog = getIntLog(bi, intSquare); 1046995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (intLog != 0) { 1047995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = 1048995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.add(new BoundedRational(intLog), 1049995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational.HALF); 1050995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 1051995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor, sLogs[intSquare]); 1052995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1053995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1054995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1055995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1056995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1057995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1058995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1059995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().ln()); 1060995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1061995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1062995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal exp() { 1063995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (definitelyEquals(ZERO)) { 1064995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return ONE; 1065995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1066e7111cfb93a073227a360f219ef3e395cb99f55dHans Boehm if (definitelyEquals(ONE)) { 1067e7111cfb93a073227a360f219ef3e395cb99f55dHans Boehm // Avoid redundant computations, and ensure we recognize all instances as equal. 1068e7111cfb93a073227a360f219ef3e395cb99f55dHans Boehm return E; 1069e7111cfb93a073227a360f219ef3e395cb99f55dHans Boehm } 1070995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final BoundedRational crExp = getExp(mCrFactor); 1071995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (crExp != null) { 1072995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mRatFactor.signum() < 0) { 1073995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return negate().exp().inverse(); 1074995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1075995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm boolean needSqrt = false; 1076995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational ratExponent = mRatFactor; 1077995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger asBI = BoundedRational.asBigInteger(ratExponent); 1078995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (asBI == null) { 1079995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // check for multiple of one half. 1080995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm needSqrt = true; 1081995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm ratExponent = BoundedRational.multiply(ratExponent, BoundedRational.TWO); 1082995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1083995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = BoundedRational.pow(crExp, ratExponent); 1084995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (nRatFactor != null) { 1085995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm UnifiedReal result = new UnifiedReal(nRatFactor); 1086995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (needSqrt) { 1087995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm result = result.sqrt(); 1088995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1089995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return result; 1090995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1091995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1092995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(crValue().exp()); 1093995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1094995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1095995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1096995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 1097995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Generalized factorial. 1098995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Compute n * (n - step) * (n - 2 * step) * etc. This can be used to compute factorial a bit 1099995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * faster, especially if BigInteger uses sub-quadratic multiplication. 1100995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 1101995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm private static BigInteger genFactorial(long n, long step) { 1102995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n > 4 * step) { 1103995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger prod1 = genFactorial(n, 2 * step); 1104995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (Thread.interrupted()) { 1105995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new CR.AbortedException(); 1106995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1107995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger prod2 = genFactorial(n - step, 2 * step); 1108995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (Thread.interrupted()) { 1109995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new CR.AbortedException(); 1110995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1111995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return prod1.multiply(prod2); 1112995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 1113995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (n == 0) { 1114995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return BigInteger.ONE; 1115995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1116995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger res = BigInteger.valueOf(n); 1117995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm for (long i = n - step; i > 1; i -= step) { 1118995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm res = res.multiply(BigInteger.valueOf(i)); 1119995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1120995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return res; 1121995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1122995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1123995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1124995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1125995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 1126995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Factorial function. 1127995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Fails if argument is clearly not an integer. 1128995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * May round to nearest integer if value is close. 1129995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 1130995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public UnifiedReal fact() { 1131995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger asBI = bigIntegerValue(); 1132995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (asBI == null) { 1133995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm asBI = crValue().get_appr(0); // Correct if it was an integer. 1134995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (!approxEquals(new UnifiedReal(asBI), DEFAULT_COMPARE_TOLERANCE)) { 1135995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("Non-integral factorial argument"); 1136995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1137995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1138995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (asBI.signum() < 0) { 1139995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("Negative factorial argument"); 1140995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1141995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (asBI.bitLength() > 20) { 1142995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Will fail. LongValue() may not work. Punt now. 1143995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm throw new ArithmeticException("Factorial argument too big"); 1144995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1145995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BigInteger biResult = genFactorial(asBI.longValue(), 1); 1146995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm BoundedRational nRatFactor = new BoundedRational(biResult); 1147995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return new UnifiedReal(nRatFactor); 1148995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1149995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1150995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 1151995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return the number of decimal digits to the right of the decimal point required to represent 1152995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * the argument exactly. 1153995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return Integer.MAX_VALUE if that's not possible. Never returns a value less than zero, even 1154995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * if r is a power of ten. 1155995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 1156995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int digitsRequired() { 1157995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (mCrFactor == CR_ONE || mRatFactor.signum() == 0) { 1158995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return BoundedRational.digitsRequired(mRatFactor); 1159995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 1160995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return Integer.MAX_VALUE; 1161995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1162995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1163995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1164995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 1165995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return an upper bound on the number of leading zero bits. 1166995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * These are the number of 0 bits 1167995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * to the right of the binary point and to the left of the most significant digit. 1168995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Return Integer.MAX_VALUE if we cannot bound it. 1169995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 1170995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public int leadingBinaryZeroes() { 1171995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isNamed(mCrFactor)) { 1172995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Only ln(2) is smaller than one, and could possibly add one zero bit. 1173995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm // Adding 3 gives us a somewhat sloppy upper bound. 1174995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm final int wholeBits = mRatFactor.wholeNumberBits(); 1175995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (wholeBits == Integer.MIN_VALUE) { 1176995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return Integer.MAX_VALUE; 1177995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1178995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (wholeBits >= 3) { 1179995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return 0; 1180995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 1181995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return -wholeBits + 3; 1182995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1183995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1184995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return Integer.MAX_VALUE; 1185995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1186995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm 1187995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm /** 1188995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * Is the number of bits to the left of the decimal point greater than bound? 1189995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm * The result is inexact: We roughly approximate the whole number bits. 1190995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm */ 1191995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm public boolean approxWholeNumberBitsGreaterThan(int bound) { 1192995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm if (isNamed(mCrFactor)) { 1193995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return mRatFactor.wholeNumberBits() > bound; 1194995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } else { 1195995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm return crValue().get_appr(bound - 2).bitLength() > 2; 1196995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1197995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm } 1198995e5eb3b58a5cbef2686feb34f53493c1825866Hans Boehm} 1199