11d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert/* 21d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Copyright (C) 2011 The Guava Authors 31d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 41d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Licensed under the Apache License, Version 2.0 (the "License"); 51d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * you may not use this file except in compliance with the License. 61d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * You may obtain a copy of the License at 71d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 81d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * http://www.apache.org/licenses/LICENSE-2.0 91d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Unless required by applicable law or agreed to in writing, software 111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * distributed under the License is distributed on an "AS IS" BASIS, 121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * See the License for the specific language governing permissions and 141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * limitations under the License. 151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertpackage com.google.common.math; 181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.base.Preconditions.checkArgument; 201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.base.Preconditions.checkNotNull; 211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkNonNegative; 221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkPositive; 231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkRoundingUnnecessary; 241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.math.RoundingMode.CEILING; 251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.math.RoundingMode.FLOOR; 261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.math.RoundingMode.HALF_EVEN; 271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport com.google.common.annotations.Beta; 291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport com.google.common.annotations.VisibleForTesting; 301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.BigDecimal; 321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.BigInteger; 331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.RoundingMode; 341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.util.ArrayList; 351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.util.List; 361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert/** 381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * A class for arithmetic on values of type {@code BigInteger}. 391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>The implementations of many methods in this class are based on material from Henry S. Warren, 411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002). 421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>Similar functionality for {@code int} and for {@code long} can be found in 441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@link IntMath} and {@link LongMath} respectively. 451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @author Louis Wasserman 471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @since 11.0 481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert@Beta 501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertpublic final class BigIntegerMath { 511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code true} if {@code x} represents a power of two. 531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static boolean isPowerOfTwo(BigInteger x) { 551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNotNull(x); 561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return x.signum() > 0 && x.getLowestSetBit() == x.bitLength() - 1; 571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode. 611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x <= 0} 631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not a power of two 651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static int log2(BigInteger x, RoundingMode mode) { 681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkPositive("x", checkNotNull(x)); 691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int logFloor = x.bitLength() - 1; 701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(isPowerOfTwo(x)); // fall through 731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return logFloor; 761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return isPowerOfTwo(x) ? logFloor : logFloor + 1; 801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (logFloor < SQRT2_PRECOMPUTE_THRESHOLD) { 851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger halfPower = SQRT2_PRECOMPUTED_BITS.shiftRight( 861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert SQRT2_PRECOMPUTE_THRESHOLD - logFloor); 871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (x.compareTo(halfPower) <= 0) { 881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return logFloor; 891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } else { 901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return logFloor + 1; 911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5 951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * To determine which side of logFloor.5 the logarithm is, we compare x^2 to 2^(2 * 971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * logFloor + 1). 981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger x2 = x.pow(2); 1001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int logX2Floor = x2.bitLength() - 1; 1011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (logX2Floor < 2 * logFloor + 1) ? logFloor : logFloor + 1; 1021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 1041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 1051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 1091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * The maximum number of bits in a square root for which we'll precompute an explicit half power 1101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * of two. This can be any value, but higher values incur more class load time and linearly 1111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * increasing memory consumption. 1121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 1131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @VisibleForTesting static final int SQRT2_PRECOMPUTE_THRESHOLD = 256; 1141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @VisibleForTesting static final BigInteger SQRT2_PRECOMPUTED_BITS = 1161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert new BigInteger("16a09e667f3bcc908b2fb1366ea957d3e3adec17512775099da2f590b0667322a", 16); 1171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 1191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode. 1201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 1211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x <= 0} 1221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 1231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not a power of ten 1241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 1251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 1261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static int log10(BigInteger x, RoundingMode mode) { 1271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkPositive("x", x); 1281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (fitsInLong(x)) { 1291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return LongMath.log10(x.longValue(), mode); 1301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // capacity of 10 suffices for all x <= 10^(2^10). 1331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert List<BigInteger> powersOf10 = new ArrayList<BigInteger>(10); 1341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger powerOf10 = BigInteger.TEN; 1351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert while (x.compareTo(powerOf10) >= 0) { 1361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert powersOf10.add(powerOf10); 1371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert powerOf10 = powerOf10.pow(2); 1381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger floorPow = BigInteger.ONE; 1401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int floorLog = 0; 1411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert for (int i = powersOf10.size() - 1; i >= 0; i--) { 1421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger powOf10 = powersOf10.get(i); 1431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert floorLog *= 2; 1441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger tenPow = powOf10.multiply(floorPow); 1451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (x.compareTo(tenPow) >= 0) { 1461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert floorPow = tenPow; 1471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert floorLog++; 1481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 1511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 1521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(floorPow.equals(x)); 1531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // fall through 1541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 1551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 1561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return floorLog; 1571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 1591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 1601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return floorPow.equals(x) ? floorLog : floorLog + 1; 1611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 1631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 1641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 1651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Since sqrt(10) is irrational, log10(x) - floorLog can never be exactly 0.5 1661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger x2 = x.pow(2); 1671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger halfPowerSquared = floorPow.pow(2).multiply(BigInteger.TEN); 1681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (x2.compareTo(halfPowerSquared) <= 0) ? floorLog : floorLog + 1; 1691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 1701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 1711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 1751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the square root of {@code x}, rounded with the specified rounding mode. 1761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 1771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x < 0} 1781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and 1791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code sqrt(x)} is not an integer 1801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 1811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 1821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static BigInteger sqrt(BigInteger x, RoundingMode mode) { 1831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("x", x); 1841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (fitsInLong(x)) { 1851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return BigInteger.valueOf(LongMath.sqrt(x.longValue(), mode)); 1861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger sqrtFloor = sqrtFloor(x); 1881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 1891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 1901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(sqrtFloor.pow(2).equals(x)); // fall through 1911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 1921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 1931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return sqrtFloor; 1941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 1951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 1961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return sqrtFloor.pow(2).equals(x) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE); 1971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 1981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 1991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 2001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger halfSquare = sqrtFloor.pow(2).add(sqrtFloor); 2011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 2021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both 2031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * x and halfSquare are integers, this is equivalent to testing whether or not x <= 2041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * halfSquare. 2051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (halfSquare.compareTo(x) >= 0) ? sqrtFloor : sqrtFloor.add(BigInteger.ONE); 2071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 2081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 2091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert private static BigInteger sqrtFloor(BigInteger x) { 2131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 2141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Adapted from Hacker's Delight, Figure 11-1. 2151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Using DoubleUtils.bigToDouble, getting a double approximation of x is extremely fast, and 2171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * then we can get a double approximation of the square root. Then, we iteratively improve this 2181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * guess with an application of Newton's method, which sets guess := (guess + (x / guess)) / 2. 2191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * This iteration has the following two properties: 2201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * a) every iteration (except potentially the first) has guess >= floor(sqrt(x)). This is 2221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * because guess' is the arithmetic mean of guess and x / guess, sqrt(x) is the geometric mean, 2231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * and the arithmetic mean is always higher than the geometric mean. 2241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * b) this iteration converges to floor(sqrt(x)). In fact, the number of correct digits doubles 2261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * with each iteration, so this algorithm takes O(log(digits)) iterations. 2271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * We start out with a double-precision approximation, which may be higher or lower than the 2291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * true value. Therefore, we perform at least one Newton iteration to get a guess that's 2301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * definitely >= floor(sqrt(x)), and then continue the iteration until we reach a fixed point. 2311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger sqrt0; 2331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int log2 = log2(x, FLOOR); 2341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if(log2 < DoubleUtils.MAX_DOUBLE_EXPONENT) { 2351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert sqrt0 = sqrtApproxWithDoubles(x); 2361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } else { 2371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int shift = (log2 - DoubleUtils.SIGNIFICAND_BITS) & ~1; // even! 2381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 2391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * We have that x / 2^shift < 2^54. Our initial approximation to sqrtFloor(x) will be 2401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2^(shift/2) * sqrtApproxWithDoubles(x / 2^shift). 2411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert sqrt0 = sqrtApproxWithDoubles(x.shiftRight(shift)).shiftLeft(shift >> 1); 2431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1); 2451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (sqrt0.equals(sqrt1)) { 2461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return sqrt0; 2471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert do { 2491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert sqrt0 = sqrt1; 2501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert sqrt1 = sqrt0.add(x.divide(sqrt0)).shiftRight(1); 2511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } while (sqrt1.compareTo(sqrt0) < 0); 2521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return sqrt0; 2531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert private static BigInteger sqrtApproxWithDoubles(BigInteger x) { 2561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return DoubleMath.roundToBigInteger(Math.sqrt(DoubleUtils.bigToDouble(x)), HALF_EVEN); 2571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 2601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the result of dividing {@code p} by {@code q}, rounding using the specified 2611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code RoundingMode}. 2621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a} 2641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not an integer multiple of {@code b} 2651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static BigInteger divide(BigInteger p, BigInteger q, RoundingMode mode){ 2671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigDecimal pDec = new BigDecimal(p); 2681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigDecimal qDec = new BigDecimal(q); 2691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return pDec.divide(qDec, 0, mode).toBigIntegerExact(); 2701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 2731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code n!}, that is, the product of the first {@code n} positive 2741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * integers, or {@code 1} if {@code n == 0}. 2751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p><b>Warning</b>: the result takes <i>O(n log n)</i> space, so use cautiously. 2771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>This uses an efficient binary recursive algorithm to compute the factorial 2791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * with balanced multiplies. It also removes all the 2s from the intermediate 2801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * products (shifting them back in at the end). 2811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code n < 0} 2831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static BigInteger factorial(int n) { 2851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("n", n); 2861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // If the factorial is small enough, just use LongMath to do it. 2881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (n < LongMath.FACTORIALS.length) { 2891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return BigInteger.valueOf(LongMath.FACTORIALS[n]); 2901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Pre-allocate space for our list of intermediate BigIntegers. 2931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int approxSize = IntMath.divide(n * IntMath.log2(n, CEILING), Long.SIZE, CEILING); 2941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert ArrayList<BigInteger> bignums = new ArrayList<BigInteger>(approxSize); 2951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Start from the pre-computed maximum long factorial. 2971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int startingNumber = LongMath.FACTORIALS.length; 2981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long product = LongMath.FACTORIALS[startingNumber - 1]; 2991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Strip off 2s from this value. 3001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int shift = Long.numberOfTrailingZeros(product); 3011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert product >>= shift; 3021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Use floor(log2(num)) + 1 to prevent overflow of multiplication. 3041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int productBits = LongMath.log2(product, FLOOR) + 1; 3051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int bits = LongMath.log2(startingNumber, FLOOR) + 1; 3061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Check for the next power of two boundary, to save us a CLZ operation. 3071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int nextPowerOfTwo = 1 << (bits - 1); 3081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Iteratively multiply the longs as big as they can go. 3101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert for (long num = startingNumber; num <= n; num++) { 3111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Check to see if the floor(log2(num)) + 1 has changed. 3121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if ((num & nextPowerOfTwo) != 0) { 3131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert nextPowerOfTwo <<= 1; 3141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert bits++; 3151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Get rid of the 2s in num. 3171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int tz = Long.numberOfTrailingZeros(num); 3181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long normalizedNum = num >> tz; 3191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert shift += tz; 3201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Adjust floor(log2(num)) + 1. 3211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int normalizedBits = bits - tz; 3221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // If it won't fit in a long, then we store off the intermediate product. 3231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (normalizedBits + productBits >= Long.SIZE) { 3241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert bignums.add(BigInteger.valueOf(product)); 3251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert product = 1; 3261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert productBits = 0; 3271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert product *= normalizedNum; 3291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert productBits = LongMath.log2(product, FLOOR) + 1; 3301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Check for leftovers. 3321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (product > 1) { 3331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert bignums.add(BigInteger.valueOf(product)); 3341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Efficiently multiply all the intermediate products together. 3361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return listProduct(bignums).shiftLeft(shift); 3371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert static BigInteger listProduct(List<BigInteger> nums) { 3401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return listProduct(nums, 0, nums.size()); 3411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert static BigInteger listProduct(List<BigInteger> nums, int start, int end) { 3441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (end - start) { 3451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 0: 3461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return BigInteger.ONE; 3471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 1: 3481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return nums.get(start); 3491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 2: 3501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return nums.get(start).multiply(nums.get(start + 1)); 3511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 3: 3521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return nums.get(start).multiply(nums.get(start + 1)).multiply(nums.get(start + 2)); 3531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 3541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Otherwise, split the list in half and recursively do this. 3551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int m = (end + start) >>> 1; 3561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return listProduct(nums, start, m).multiply(listProduct(nums, m, end)); 3571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 3611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and 3621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code k}, that is, {@code n! / (k! (n - k)!)}. 3631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 3641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p><b>Warning</b>: the result can take as much as <i>O(k log n)</i> space. 3651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 3661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n} 3671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 3681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static BigInteger binomial(int n, int k) { 3691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("n", n); 3701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("k", k); 3711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkArgument(k <= n, "k (%s) > n (%s)", k, n); 3721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (k > (n >> 1)) { 3731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert k = n - k; 3741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (k < LongMath.BIGGEST_BINOMIALS.length && n <= LongMath.BIGGEST_BINOMIALS[k]) { 3761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return BigInteger.valueOf(LongMath.binomial(n, k)); 3771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert BigInteger result = BigInteger.ONE; 3791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert for (int i = 0; i < k; i++) { 3801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert result = result.multiply(BigInteger.valueOf(n - i)); 3811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert result = result.divide(BigInteger.valueOf(i + 1)); 3821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return result; 3841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Returns true if BigInteger.valueOf(x.longValue()).equals(x). 3871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert static boolean fitsInLong(BigInteger x) { 3881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return x.bitLength() <= Long.SIZE - 1; 3891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert private BigIntegerMath() {} 3921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert} 393