001 /*
002 * Copyright (C) 2011 The Guava Authors
003 *
004 * Licensed under the Apache License, Version 2.0 (the "License");
005 * you may not use this file except in compliance with the License.
006 * You may obtain a copy of the License at
007 *
008 * http://www.apache.org/licenses/LICENSE-2.0
009 *
010 * Unless required by applicable law or agreed to in writing, software
011 * distributed under the License is distributed on an "AS IS" BASIS,
012 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013 * See the License for the specific language governing permissions and
014 * limitations under the License.
015 */
016
017 package com.google.common.math;
018
019 import static com.google.common.base.Preconditions.checkArgument;
020 import static com.google.common.base.Preconditions.checkNotNull;
021 import static com.google.common.math.MathPreconditions.checkNonNegative;
022 import static com.google.common.math.MathPreconditions.checkPositive;
023 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
024 import static java.math.RoundingMode.CEILING;
025 import static java.math.RoundingMode.FLOOR;
026 import static java.math.RoundingMode.HALF_EVEN;
027
028 import com.google.common.annotations.Beta;
029 import com.google.common.annotations.GwtCompatible;
030 import com.google.common.annotations.GwtIncompatible;
031 import com.google.common.annotations.VisibleForTesting;
032
033 import java.math.BigDecimal;
034 import java.math.BigInteger;
035 import java.math.RoundingMode;
036 import java.util.ArrayList;
037 import java.util.List;
038
039 /**
040 * A class for arithmetic on values of type {@code BigInteger}.
041 *
042 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
043 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
044 *
045 * <p>Similar functionality for {@code int} and for {@code long} can be found in
046 * {@link IntMath} and {@link LongMath} respectively.
047 *
048 * @author Louis Wasserman
049 * @since 11.0
050 */
051 @Beta
052 @GwtCompatible(emulated = true)
053 public final class BigIntegerMath {
054 /**
055 * Returns {@code true} if {@code x} represents a power of two.
056 */
057 public static boolean isPowerOfTwo(BigInteger x) {
058 checkNotNull(x);
059 return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1;
060 }
061
062 /**
063 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
064 *
065 * @throws IllegalArgumentException if {@code x <= 0}
066 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
067 * is not a power of two
068 */
069 @SuppressWarnings("fallthrough")
070 public static int log2(BigInteger x, RoundingMode mode) {
071 checkPositive("x", checkNotNull(x));
072 int logFloor = x.bitLength() - 1;
073 switch (mode) {
074 case UNNECESSARY:
075 checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through
076 case DOWN:
077 case FLOOR:
078 return logFloor;
079
080 case UP:
081 case CEILING:
082 return isPowerOfTwo(x) ? logFloor : logFloor + 1;
083
084 case HALF_DOWN:
085 case HALF_UP:
086 case HALF_EVEN:
087 if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) {
088 BigInteger halfPower = SQRT2_PRECOMPUTED_BITS.shiftRight(
089 SQRT2_PRECOMPUTE_THRESHOLD - logFloor);
090 if (x.compareTo(halfPower) <= 0) {
091 return logFloor;
092 } else {
093 return logFloor + 1;
094 }
095 }
096 /*
097 * Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
098 *
099 * To determine which side of logFloor.5 the logarithm is, we compare x^2 to 2^(2 *
100 * logFloor + 1).
101 */
102 BigInteger x2 = x.pow(2);
103 int logX2Floor = x2.bitLength() - 1;
104 return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1;
105
106 default:
107 throw new AssertionError();
108 }
109 }
110
111 /*
112 * The maximum number of bits in a square root for which we'll precompute an explicit half power
113 * of two. This can be any value, but higher values incur more class load time and linearly
114 * increasing memory consumption.
115 */
116 @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256;
117
118 @VisibleForTesting static final BigInteger SQRT2_PRECOMPUTED_BITS =
119 new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16);
120
121 /**
122 * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
123 *
124 * @throws IllegalArgumentException if {@code x <= 0}
125 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
126 * is not a power of ten
127 */
128 @GwtIncompatible("TODO")
129 @SuppressWarnings("fallthrough")
130 public static int log10(BigInteger x, RoundingMode mode) {
131 checkPositive("x", x);
132 if (fitsInLong(x)) {
133 return LongMath.log10(x.longValue(), mode);
134 }
135
136 // capacity of 10 suffices for all x <= 10^(2^10).
137 List<BigInteger> powersOf10 = new ArrayList<BigInteger>(10);
138 BigInteger powerOf10 = BigInteger.TEN;
139 while (x.compareTo(powerOf10) >= 0) {
140 powersOf10.add(powerOf10);
141 powerOf10 = powerOf10.pow(2);
142 }
143 BigInteger floorPow = BigInteger.ONE;
144 int floorLog = 0;
145 for (int i = powersOf10.size() - 1; i >= 0; i--) {
146 BigInteger powOf10 = powersOf10.get(i);
147 floorLog *= 2;
148 BigInteger tenPow = powOf10.multiply(floorPow);
149 if (x.compareTo(tenPow) >= 0) {
150 floorPow = tenPow;
151 floorLog++;
152 }
153 }
154 switch (mode) {
155 case UNNECESSARY:
156 checkRoundingUnnecessary(floorPow.equals(x));
157 // fall through
158 case FLOOR:
159 case DOWN:
160 return floorLog;
161
162 case CEILING:
163 case UP:
164 return floorPow.equals(x) ? floorLog : floorLog + 1;
165
166 case HALF_DOWN:
167 case HALF_UP:
168 case HALF_EVEN:
169 // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5
170 BigInteger x2 = x.pow(2);
171 BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN);
172 return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1;
173 default:
174 throw new AssertionError();
175 }
176 }
177
178 /**
179 * Returns the square root of {@code x}, rounded with the specified rounding mode.
180 *
181 * @throws IllegalArgumentException if {@code x < 0}
182 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and
183 * {@code sqrt(x)} is not an integer
184 */
185 @GwtIncompatible("TODO")
186 @SuppressWarnings("fallthrough")
187 public static BigInteger sqrt(BigInteger x, RoundingMode mode) {
188 checkNonNegative("x", x);
189 if (fitsInLong(x)) {
190 return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode));
191 }
192 BigInteger sqrtFloor = sqrtFloor(x);
193 switch (mode) {
194 case UNNECESSARY:
195 checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through
196 case FLOOR:
197 case DOWN:
198 return sqrtFloor;
199 case CEILING:
200 case UP:
201 return sqrtFloor.pow(2).equals(x) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
202 case HALF_DOWN:
203 case HALF_UP:
204 case HALF_EVEN:
205 BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor);
206 /*
207 * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both
208 * x and halfSquare are integers, this is equivalent to testing whether or not x <=
209 * halfSquare.
210 */
211 return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE);
212 default:
213 throw new AssertionError();
214 }
215 }
216
217 @GwtIncompatible("TODO")
218 private static BigInteger sqrtFloor(BigInteger x) {
219 /*
220 * Adapted from Hacker's Delight, Figure 11-1.
221 *
222 * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and
223 * then we can get a double approximation of the square root. Then, we iteratively improve this
224 * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2.
225 * This iteration has the following two properties:
226 *
227 * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is
228 * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean,
229 * and the arithmetic mean is always higher than the geometric mean.
230 *
231 * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles
232 * with each iteration, so this algorithm takes O(log(digits)) iterations.
233 *
234 * We start out with a double-precision approximation, which may be higher or lower than the
235 * true value. Therefore, we perform at least one Newton iteration to get a guess that's
236 * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point.
237 */
238 BigInteger sqrt0;
239 int log2 = log2(x, FLOOR);
240 if(log2 < Double.MAX_EXPONENT) {
241 sqrt0 = sqrtApproxWithDoubles(x);
242 } else {
243 int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even!
244 /*
245 * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be
246 * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift).
247 */
248 sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1);
249 }
250 BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
251 if (sqrt0.equals(sqrt1)) {
252 return sqrt0;
253 }
254 do {
255 sqrt0 = sqrt1;
256 sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1);
257 } while (sqrt1.compareTo(sqrt0) < 0);
258 return sqrt0;
259 }
260
261 @GwtIncompatible("TODO")
262 private static BigInteger sqrtApproxWithDoubles(BigInteger x) {
263 return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN);
264 }
265
266 /**
267 * Returns the result of dividing {@code p} by {@code q}, rounding using the specified
268 * {@code RoundingMode}.
269 *
270 * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
271 * is not an integer multiple of {@code b}
272 */
273 @GwtIncompatible("TODO")
274 public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode){
275 BigDecimal pDec = new BigDecimal(p);
276 BigDecimal qDec = new BigDecimal(q);
277 return pDec.divide(qDec, 0, mode).toBigIntegerExact();
278 }
279
280 /**
281 * Returns {@code n!}, that is, the product of the first {@code n} positive
282 * integers, or {@code 1} if {@code n == 0}.
283 *
284 * <p><b>Warning</b>: the result takes <i>O(n log n)</i> space, so use cautiously.
285 *
286 * <p>This uses an efficient binary recursive algorithm to compute the factorial
287 * with balanced multiplies. It also removes all the 2s from the intermediate
288 * products (shifting them back in at the end).
289 *
290 * @throws IllegalArgumentException if {@code n < 0}
291 */
292 public static BigInteger factorial(int n) {
293 checkNonNegative("n", n);
294
295 // If the factorial is small enough, just use LongMath to do it.
296 if (n < LongMath.FACTORIALS.length) {
297 return BigInteger.valueOf(LongMath.FACTORIALS[n]);
298 }
299
300 // Pre-allocate space for our list of intermediate BigIntegers.
301 int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING);
302 ArrayList<BigInteger> bignums = new ArrayList<BigInteger>(approxSize);
303
304 // Start from the pre-computed maximum long factorial.
305 int startingNumber = LongMath.FACTORIALS.length;
306 long product = LongMath.FACTORIALS[startingNumber - 1];
307 // Strip off 2s from this value.
308 int shift = Long.numberOfTrailingZeros(product);
309 product >>= shift;
310
311 // Use floor(log2(num)) + 1 to prevent overflow of multiplication.
312 int productBits = LongMath.log2(product, FLOOR) + 1;
313 int bits = LongMath.log2(startingNumber, FLOOR) + 1;
314 // Check for the next power of two boundary, to save us a CLZ operation.
315 int nextPowerOfTwo = 1 << (bits - 1);
316
317 // Iteratively multiply the longs as big as they can go.
318 for (long num = startingNumber; num <= n; num++) {
319 // Check to see if the floor(log2(num)) + 1 has changed.
320 if ((num & nextPowerOfTwo) != 0) {
321 nextPowerOfTwo <<= 1;
322 bits++;
323 }
324 // Get rid of the 2s in num.
325 int tz = Long.numberOfTrailingZeros(num);
326 long normalizedNum = num >> tz;
327 shift += tz;
328 // Adjust floor(log2(num)) + 1.
329 int normalizedBits = bits - tz;
330 // If it won't fit in a long, then we store off the intermediate product.
331 if (normalizedBits + productBits >= Long.SIZE) {
332 bignums.add(BigInteger.valueOf(product));
333 product = 1;
334 productBits = 0;
335 }
336 product *= normalizedNum;
337 productBits = LongMath.log2(product, FLOOR) + 1;
338 }
339 // Check for leftovers.
340 if (product > 1) {
341 bignums.add(BigInteger.valueOf(product));
342 }
343 // Efficiently multiply all the intermediate products together.
344 return listProduct(bignums).shiftLeft(shift);
345 }
346
347 static BigInteger listProduct(List<BigInteger> nums) {
348 return listProduct(nums, 0, nums.size());
349 }
350
351 static BigInteger listProduct(List<BigInteger> nums, int start, int end) {
352 switch (end - start) {
353 case 0:
354 return BigInteger.ONE;
355 case 1:
356 return nums.get(start);
357 case 2:
358 return nums.get(start).multiply(nums.get(start + 1));
359 case 3:
360 return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2));
361 default:
362 // Otherwise, split the list in half and recursively do this.
363 int m = (end + start) >>> 1;
364 return listProduct(nums, start, m).multiply(listProduct(nums, m, end));
365 }
366 }
367
368 /**
369 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
370 * {@code k}, that is, {@code n! / (k! (n - k)!)}.
371 *
372 * <p><b>Warning</b>: the result can take as much as <i>O(k log n)</i> space.
373 *
374 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
375 */
376 public static BigInteger binomial(int n, int k) {
377 checkNonNegative("n", n);
378 checkNonNegative("k", k);
379 checkArgument(k <= n, "k (%s) > n (%s)", k, n);
380 if (k > (n >> 1)) {
381 k = n - k;
382 }
383 if (k < LongMath.BIGGEST_BINOMIALS.length && n <= LongMath.BIGGEST_BINOMIALS[k]) {
384 return BigInteger.valueOf(LongMath.binomial(n, k));
385 }
386 BigInteger result = BigInteger.ONE;
387 for (int i = 0; i < k; i++) {
388 result = result.multiply(BigInteger.valueOf(n - i));
389 result = result.divide(BigInteger.valueOf(i + 1));
390 }
391 return result;
392 }
393
394 // Returns true if BigInteger.valueOf(x.longValue()).equals(x).
395 @GwtIncompatible("TODO")
396 static boolean fitsInLong(BigInteger x) {
397 return x.bitLength() <= Long.SIZE - 1;
398 }
399
400 private BigIntegerMath() {}
401 }