001/* 002 * Copyright (C) 2011 The Guava Authors 003 * 004 * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except 005 * in compliance with the License. You may obtain a copy of the License at 006 * 007 * http://www.apache.org/licenses/LICENSE-2.0 008 * 009 * Unless required by applicable law or agreed to in writing, software distributed under the License 010 * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express 011 * or implied. See the License for the specific language governing permissions and limitations under 012 * the License. 013 */ 014 015package com.google.common.math; 016 017import static com.google.common.base.Preconditions.checkArgument; 018import static com.google.common.base.Preconditions.checkNotNull; 019import static com.google.common.math.MathPreconditions.checkNonNegative; 020import static com.google.common.math.MathPreconditions.checkPositive; 021import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary; 022import static java.math.RoundingMode.CEILING; 023import static java.math.RoundingMode.FLOOR; 024import static java.math.RoundingMode.HALF_EVEN; 025 026import com.google.common.annotations.GwtCompatible; 027import com.google.common.annotations.GwtIncompatible; 028import com.google.common.annotations.VisibleForTesting; 029import java.math.BigDecimal; 030import java.math.BigInteger; 031import java.math.RoundingMode; 032import java.util.ArrayList; 033import java.util.List; 034 035/** 036 * A class for arithmetic on values of type {@code BigInteger}. 037 * 038 * <p>The implementations of many methods in this class are based on material from Henry S. Warren, 039 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002). 040 * 041 * <p>Similar functionality for {@code int} and for {@code long} can be found in {@link IntMath} and 042 * {@link LongMath} respectively. 043 * 044 * @author Louis Wasserman 045 * @since 11.0 046 */ 047@GwtCompatible(emulated = true) 048@ElementTypesAreNonnullByDefault 049public final class BigIntegerMath { 050 /** 051 * Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to 052 * {@code BigInteger.valueOf(2).pow(log2(x, CEILING))}. 053 * 054 * @throws IllegalArgumentException if {@code x <= 0} 055 * @since 20.0 056 */ 057 public static BigInteger ceilingPowerOfTwo(BigInteger x) { 058 return BigInteger.ZERO.setBit(log2(x, CEILING)); 059 } 060 061 /** 062 * Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code 063 * BigInteger.valueOf(2).pow(log2(x, FLOOR))}. 064 * 065 * @throws IllegalArgumentException if {@code x <= 0} 066 * @since 20.0 067 */ 068 public static BigInteger floorPowerOfTwo(BigInteger x) { 069 return BigInteger.ZERO.setBit(log2(x, FLOOR)); 070 } 071 072 /** Returns {@code true} if {@code x} represents a power of two. */ 073 public static boolean isPowerOfTwo(BigInteger x) { 074 checkNotNull(x); 075 return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1; 076 } 077 078 /** 079 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode. 080 * 081 * @throws IllegalArgumentException if {@code x <= 0} 082 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 083 * is not a power of two 084 */ 085 @SuppressWarnings("fallthrough") 086 // TODO(kevinb): remove after this warning is disabled globally 087 public static int log2(BigInteger x, RoundingMode mode) { 088 checkPositive("x", checkNotNull(x)); 089 int logFloor = x.bitLength() - 1; 090 switch (mode) { 091 case UNNECESSARY: 092 checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through 093 case DOWN: 094 case FLOOR: 095 return logFloor; 096 097 case UP: 098 case CEILING: 099 return isPowerOfTwo(x) ? logFloor : logFloor + 1; 100 101 case HALF_DOWN: 102 case HALF_UP: 103 case HALF_EVEN: 104 if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) { 105 BigInteger halfPower = 106 SQRT2_PRECOMPUTED_BITS.shiftRight(SQRT2_PRECOMPUTE_THRESHOLD - logFloor); 107 if (x.compareTo(halfPower) <= 0) { 108 return logFloor; 109 } else { 110 return logFloor + 1; 111 } 112 } 113 // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5 114 // 115 // To determine which side of logFloor.5 the logarithm is, 116 // we compare x^2 to 2^(2 * logFloor + 1). 117 BigInteger x2 = x.pow(2); 118 int logX2Floor = x2.bitLength() - 1; 119 return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1; 120 121 default: 122 throw new AssertionError(); 123 } 124 } 125 126 /* 127 * The maximum number of bits in a square root for which we'll precompute an explicit half power 128 * of two. This can be any value, but higher values incur more class load time and linearly 129 * increasing memory consumption. 130 */ 131 @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256; 132 133 @VisibleForTesting 134 static final BigInteger SQRT2_PRECOMPUTED_BITS = 135 new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16); 136 137 /** 138 * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode. 139 * 140 * @throws IllegalArgumentException if {@code x <= 0} 141 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 142 * is not a power of ten 143 */ 144 @GwtIncompatible // TODO 145 @SuppressWarnings("fallthrough") 146 public static int log10(BigInteger x, RoundingMode mode) { 147 checkPositive("x", x); 148 if (fitsInLong(x)) { 149 return LongMath.log10(x.longValue(), mode); 150 } 151 152 int approxLog10 = (int) (log2(x, FLOOR) * LN_2 / LN_10); 153 BigInteger approxPow = BigInteger.TEN.pow(approxLog10); 154 int approxCmp = approxPow.compareTo(x); 155 156 /* 157 * We adjust approxLog10 and approxPow until they're equal to floor(log10(x)) and 158 * 10^floor(log10(x)). 159 */ 160 161 if (approxCmp > 0) { 162 /* 163 * The code is written so that even completely incorrect approximations will still yield the 164 * correct answer eventually, but in practice this branch should almost never be entered, and 165 * even then the loop should not run more than once. 166 */ 167 do { 168 approxLog10--; 169 approxPow = approxPow.divide(BigInteger.TEN); 170 approxCmp = approxPow.compareTo(x); 171 } while (approxCmp > 0); 172 } else { 173 BigInteger nextPow = BigInteger.TEN.multiply(approxPow); 174 int nextCmp = nextPow.compareTo(x); 175 while (nextCmp <= 0) { 176 approxLog10++; 177 approxPow = nextPow; 178 approxCmp = nextCmp; 179 nextPow = BigInteger.TEN.multiply(approxPow); 180 nextCmp = nextPow.compareTo(x); 181 } 182 } 183 184 int floorLog = approxLog10; 185 BigInteger floorPow = approxPow; 186 int floorCmp = approxCmp; 187 188 switch (mode) { 189 case UNNECESSARY: 190 checkRoundingUnnecessary(floorCmp == 0); 191 // fall through 192 case FLOOR: 193 case DOWN: 194 return floorLog; 195 196 case CEILING: 197 case UP: 198 return floorPow.equals(x) ? floorLog : floorLog + 1; 199 200 case HALF_DOWN: 201 case HALF_UP: 202 case HALF_EVEN: 203 // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5 204 BigInteger x2 = x.pow(2); 205 BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN); 206 return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1; 207 default: 208 throw new AssertionError(); 209 } 210 } 211 212 private static final double LN_10 = Math.log(10); 213 private static final double LN_2 = Math.log(2); 214 215 /** 216 * Returns the square root of {@code x}, rounded with the specified rounding mode. 217 * 218 * @throws IllegalArgumentException if {@code x < 0} 219 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code 220 * sqrt(x)} is not an integer 221 */ 222 @GwtIncompatible // TODO 223 @SuppressWarnings("fallthrough") 224 public static BigInteger sqrt(BigInteger x, RoundingMode mode) { 225 checkNonNegative("x", x); 226 if (fitsInLong(x)) { 227 return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode)); 228 } 229 BigInteger sqrtFloor = sqrtFloor(x); 230 switch (mode) { 231 case UNNECESSARY: 232 checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through 233 case FLOOR: 234 case DOWN: 235 return sqrtFloor; 236 case CEILING: 237 case UP: 238 int sqrtFloorInt = sqrtFloor.intValue(); 239 boolean sqrtFloorIsExact = 240 (sqrtFloorInt * sqrtFloorInt == x.intValue()) // fast check mod 2^32 241 && sqrtFloor.pow(2).equals(x); // slow exact check 242 return sqrtFloorIsExact ? sqrtFloor : sqrtFloor.add(BigInteger.ONE); 243 case HALF_DOWN: 244 case HALF_UP: 245 case HALF_EVEN: 246 BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor); 247 /* 248 * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x 249 * and halfSquare are integers, this is equivalent to testing whether or not x <= 250 * halfSquare. 251 */ 252 return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE); 253 default: 254 throw new AssertionError(); 255 } 256 } 257 258 @GwtIncompatible // TODO 259 private static BigInteger sqrtFloor(BigInteger x) { 260 /* 261 * Adapted from Hacker's Delight, Figure 11-1. 262 * 263 * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and 264 * then we can get a double approximation of the square root. Then, we iteratively improve this 265 * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2. 266 * This iteration has the following two properties: 267 * 268 * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is 269 * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean, 270 * and the arithmetic mean is always higher than the geometric mean. 271 * 272 * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles 273 * with each iteration, so this algorithm takes O(log(digits)) iterations. 274 * 275 * We start out with a double-precision approximation, which may be higher or lower than the 276 * true value. Therefore, we perform at least one Newton iteration to get a guess that's 277 * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point. 278 */ 279 BigInteger sqrt0; 280 int log2 = log2(x, FLOOR); 281 if (log2 < Double.MAX_EXPONENT) { 282 sqrt0 = sqrtApproxWithDoubles(x); 283 } else { 284 int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even! 285 /* 286 * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be 287 * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift). 288 */ 289 sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1); 290 } 291 BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1); 292 if (sqrt0.equals(sqrt1)) { 293 return sqrt0; 294 } 295 do { 296 sqrt0 = sqrt1; 297 sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1); 298 } while (sqrt1.compareTo(sqrt0) < 0); 299 return sqrt0; 300 } 301 302 @GwtIncompatible // TODO 303 private static BigInteger sqrtApproxWithDoubles(BigInteger x) { 304 return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN); 305 } 306 307 /** 308 * Returns {@code x}, rounded to a {@code double} with the specified rounding mode. If {@code x} 309 * is precisely representable as a {@code double}, its {@code double} value will be returned; 310 * otherwise, the rounding will choose between the two nearest representable values with {@code 311 * mode}. 312 * 313 * <p>For the case of {@link RoundingMode#HALF_DOWN}, {@code HALF_UP}, and {@code HALF_EVEN}, 314 * infinite {@code double} values are considered infinitely far away. For example, 2^2000 is not 315 * representable as a double, but {@code roundToDouble(BigInteger.valueOf(2).pow(2000), HALF_UP)} 316 * will return {@code Double.MAX_VALUE}, not {@code Double.POSITIVE_INFINITY}. 317 * 318 * <p>For the case of {@link RoundingMode#HALF_EVEN}, this implementation uses the IEEE 754 319 * default rounding mode: if the two nearest representable values are equally near, the one with 320 * the least significant bit zero is chosen. (In such cases, both of the nearest representable 321 * values are even integers; this method returns the one that is a multiple of a greater power of 322 * two.) 323 * 324 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 325 * is not precisely representable as a {@code double} 326 * @since 30.0 327 */ 328 @GwtIncompatible 329 public static double roundToDouble(BigInteger x, RoundingMode mode) { 330 return BigIntegerToDoubleRounder.INSTANCE.roundToDouble(x, mode); 331 } 332 333 @GwtIncompatible 334 private static class BigIntegerToDoubleRounder extends ToDoubleRounder<BigInteger> { 335 static final BigIntegerToDoubleRounder INSTANCE = new BigIntegerToDoubleRounder(); 336 337 private BigIntegerToDoubleRounder() {} 338 339 @Override 340 double roundToDoubleArbitrarily(BigInteger bigInteger) { 341 return DoubleUtils.bigToDouble(bigInteger); 342 } 343 344 @Override 345 int sign(BigInteger bigInteger) { 346 return bigInteger.signum(); 347 } 348 349 @Override 350 BigInteger toX(double d, RoundingMode mode) { 351 return DoubleMath.roundToBigInteger(d, mode); 352 } 353 354 @Override 355 BigInteger minus(BigInteger a, BigInteger b) { 356 return a.subtract(b); 357 } 358 } 359 360 /** 361 * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code 362 * RoundingMode}. 363 * 364 * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a} 365 * is not an integer multiple of {@code b} 366 */ 367 @GwtIncompatible // TODO 368 public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode) { 369 BigDecimal pDec = new BigDecimal(p); 370 BigDecimal qDec = new BigDecimal(q); 371 return pDec.divide(qDec, 0, mode).toBigIntegerExact(); 372 } 373 374 /** 375 * Returns {@code n!}, that is, the product of the first {@code n} positive integers, or {@code 1} 376 * if {@code n == 0}. 377 * 378 * <p><b>Warning:</b> the result takes <i>O(n log n)</i> space, so use cautiously. 379 * 380 * <p>This uses an efficient binary recursive algorithm to compute the factorial with balanced 381 * multiplies. It also removes all the 2s from the intermediate products (shifting them back in at 382 * the end). 383 * 384 * @throws IllegalArgumentException if {@code n < 0} 385 */ 386 public static BigInteger factorial(int n) { 387 checkNonNegative("n", n); 388 389 // If the factorial is small enough, just use LongMath to do it. 390 if (n < LongMath.factorials.length) { 391 return BigInteger.valueOf(LongMath.factorials[n]); 392 } 393 394 // Pre-allocate space for our list of intermediate BigIntegers. 395 int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING); 396 ArrayList<BigInteger> bignums = new ArrayList<>(approxSize); 397 398 // Start from the pre-computed maximum long factorial. 399 int startingNumber = LongMath.factorials.length; 400 long product = LongMath.factorials[startingNumber - 1]; 401 // Strip off 2s from this value. 402 int shift = Long.numberOfTrailingZeros(product); 403 product >>= shift; 404 405 // Use floor(log2(num)) + 1 to prevent overflow of multiplication. 406 int productBits = LongMath.log2(product, FLOOR) + 1; 407 int bits = LongMath.log2(startingNumber, FLOOR) + 1; 408 // Check for the next power of two boundary, to save us a CLZ operation. 409 int nextPowerOfTwo = 1 << (bits - 1); 410 411 // Iteratively multiply the longs as big as they can go. 412 for (long num = startingNumber; num <= n; num++) { 413 // Check to see if the floor(log2(num)) + 1 has changed. 414 if ((num & nextPowerOfTwo) != 0) { 415 nextPowerOfTwo <<= 1; 416 bits++; 417 } 418 // Get rid of the 2s in num. 419 int tz = Long.numberOfTrailingZeros(num); 420 long normalizedNum = num >> tz; 421 shift += tz; 422 // Adjust floor(log2(num)) + 1. 423 int normalizedBits = bits - tz; 424 // If it won't fit in a long, then we store off the intermediate product. 425 if (normalizedBits + productBits >= Long.SIZE) { 426 bignums.add(BigInteger.valueOf(product)); 427 product = 1; 428 productBits = 0; 429 } 430 product *= normalizedNum; 431 productBits = LongMath.log2(product, FLOOR) + 1; 432 } 433 // Check for leftovers. 434 if (product > 1) { 435 bignums.add(BigInteger.valueOf(product)); 436 } 437 // Efficiently multiply all the intermediate products together. 438 return listProduct(bignums).shiftLeft(shift); 439 } 440 441 static BigInteger listProduct(List<BigInteger> nums) { 442 return listProduct(nums, 0, nums.size()); 443 } 444 445 static BigInteger listProduct(List<BigInteger> nums, int start, int end) { 446 switch (end - start) { 447 case 0: 448 return BigInteger.ONE; 449 case 1: 450 return nums.get(start); 451 case 2: 452 return nums.get(start).multiply(nums.get(start + 1)); 453 case 3: 454 return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2)); 455 default: 456 // Otherwise, split the list in half and recursively do this. 457 int m = (end + start) >>> 1; 458 return listProduct(nums, start, m).multiply(listProduct(nums, m, end)); 459 } 460 } 461 462 /** 463 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and 464 * {@code k}, that is, {@code n! / (k! (n - k)!)}. 465 * 466 * <p><b>Warning:</b> the result can take as much as <i>O(k log n)</i> space. 467 * 468 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n} 469 */ 470 public static BigInteger binomial(int n, int k) { 471 checkNonNegative("n", n); 472 checkNonNegative("k", k); 473 checkArgument(k <= n, "k (%s) > n (%s)", k, n); 474 if (k > (n >> 1)) { 475 k = n - k; 476 } 477 if (k < LongMath.biggestBinomials.length && n <= LongMath.biggestBinomials[k]) { 478 return BigInteger.valueOf(LongMath.binomial(n, k)); 479 } 480 481 BigInteger accum = BigInteger.ONE; 482 483 long numeratorAccum = n; 484 long denominatorAccum = 1; 485 486 int bits = LongMath.log2(n, CEILING); 487 488 int numeratorBits = bits; 489 490 for (int i = 1; i < k; i++) { 491 int p = n - i; 492 int q = i + 1; 493 494 // log2(p) >= bits - 1, because p >= n/2 495 496 if (numeratorBits + bits >= Long.SIZE - 1) { 497 // The numerator is as big as it can get without risking overflow. 498 // Multiply numeratorAccum / denominatorAccum into accum. 499 accum = 500 accum 501 .multiply(BigInteger.valueOf(numeratorAccum)) 502 .divide(BigInteger.valueOf(denominatorAccum)); 503 numeratorAccum = p; 504 denominatorAccum = q; 505 numeratorBits = bits; 506 } else { 507 // We can definitely multiply into the long accumulators without overflowing them. 508 numeratorAccum *= p; 509 denominatorAccum *= q; 510 numeratorBits += bits; 511 } 512 } 513 return accum 514 .multiply(BigInteger.valueOf(numeratorAccum)) 515 .divide(BigInteger.valueOf(denominatorAccum)); 516 } 517 518 // Returns true if BigInteger.valueOf(x.longValue()).equals(x). 519 @GwtIncompatible // TODO 520 static boolean fitsInLong(BigInteger x) { 521 return x.bitLength() <= Long.SIZE - 1; 522 } 523 524 private BigIntegerMath() {} 525}