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.checkNoOverflow; 221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkNonNegative; 231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkPositive; 241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.math.MathPreconditions.checkRoundingUnnecessary; 251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.lang.Math.abs; 261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.lang.Math.min; 271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.math.RoundingMode.HALF_EVEN; 281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static java.math.RoundingMode.HALF_UP; 291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 307dd252788645e940eada959bdde927426e2531c9Paul Duffinimport com.google.common.annotations.GwtCompatible; 317dd252788645e940eada959bdde927426e2531c9Paul Duffinimport com.google.common.annotations.GwtIncompatible; 321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport com.google.common.annotations.VisibleForTesting; 331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.BigInteger; 351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.RoundingMode; 361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert/** 381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and 391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * named analogously to their {@code BigInteger} counterparts. 401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>The implementations of many methods in this class are based on material from Henry S. Warren, 421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002). 431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in 451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@link IntMath} and {@link BigIntegerMath} respectively. For other common operations on 461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code long} values, see {@link com.google.common.primitives.Longs}. 471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @author Louis Wasserman 491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @since 11.0 501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 517dd252788645e940eada959bdde927426e2531c9Paul Duffin@GwtCompatible(emulated = true) 521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertpublic final class LongMath { 531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, || 541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code true} if {@code x} represents a power of two. 571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>This differs from {@code Long.bitCount(x) == 1}, because 591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two. 601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static boolean isPowerOfTwo(long x) { 621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return x > 0 & (x & (x - 1)) == 0; 631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 640888a09821a98ac0680fad765217302858e70fa4Paul Duffin 650888a09821a98ac0680fad765217302858e70fa4Paul Duffin /** 660888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a 670888a09821a98ac0680fad765217302858e70fa4Paul Duffin * signed long. The implementation is branch-free, and benchmarks suggest it is measurably 680888a09821a98ac0680fad765217302858e70fa4Paul Duffin * faster than the straightforward ternary expression. 690888a09821a98ac0680fad765217302858e70fa4Paul Duffin */ 700888a09821a98ac0680fad765217302858e70fa4Paul Duffin @VisibleForTesting 710888a09821a98ac0680fad765217302858e70fa4Paul Duffin static int lessThanBranchFree(long x, long y) { 720888a09821a98ac0680fad765217302858e70fa4Paul Duffin // Returns the sign bit of x - y. 730888a09821a98ac0680fad765217302858e70fa4Paul Duffin return (int) (~~(x - y) >>> (Long.SIZE - 1)); 740888a09821a98ac0680fad765217302858e70fa4Paul Duffin } 751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode. 781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x <= 0} 801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not a power of two 821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 847dd252788645e940eada959bdde927426e2531c9Paul Duffin // TODO(kevinb): remove after this warning is disabled globally 851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static int log2(long x, RoundingMode mode) { 861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkPositive("x", x); 871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(isPowerOfTwo(x)); 901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // fall through 911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x); 941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return Long.SIZE - Long.numberOfLeadingZeros(x - 1); 981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 1001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 1011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 1021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5 1031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int leadingZeros = Long.numberOfLeadingZeros(x); 1041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros; 1051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // floor(2^(logFloor + 0.5)) 1061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int logFloor = (Long.SIZE - 1) - leadingZeros; 1070888a09821a98ac0680fad765217302858e70fa4Paul Duffin return logFloor + lessThanBranchFree(cmp, x); 1081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 1101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError("impossible"); 1111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** The biggest half power of two that fits into an unsigned long */ 1150888a09821a98ac0680fad765217302858e70fa4Paul Duffin @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L; 1161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 1181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode. 1191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 1201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x <= 0} 1211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 1221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not a power of ten 1231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 1247dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 1251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 1267dd252788645e940eada959bdde927426e2531c9Paul Duffin // TODO(kevinb): remove after this warning is disabled globally 1271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static int log10(long x, RoundingMode mode) { 1281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkPositive("x", x); 1291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int logFloor = log10Floor(x); 1307dd252788645e940eada959bdde927426e2531c9Paul Duffin long floorPow = powersOf10[logFloor]; 1311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 1321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 1331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(x == floorPow); 1341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // fall through 1351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 1361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 1371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return logFloor; 1381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 1391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 1400888a09821a98ac0680fad765217302858e70fa4Paul Duffin return logFloor + lessThanBranchFree(floorPow, x); 1411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 1421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 1431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 1441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5 1450888a09821a98ac0680fad765217302858e70fa4Paul Duffin return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x); 1461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 1471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 1481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1517dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 1521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert static int log10Floor(long x) { 1537dd252788645e940eada959bdde927426e2531c9Paul Duffin /* 1547dd252788645e940eada959bdde927426e2531c9Paul Duffin * Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation. 1557dd252788645e940eada959bdde927426e2531c9Paul Duffin * 1567dd252788645e940eada959bdde927426e2531c9Paul Duffin * The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))), 1577dd252788645e940eada959bdde927426e2531c9Paul Duffin * we can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x)) 1587dd252788645e940eada959bdde927426e2531c9Paul Duffin * is 6, then 64 <= x < 128, so floor(log10(x)) is either 1 or 2. 1597dd252788645e940eada959bdde927426e2531c9Paul Duffin */ 1607dd252788645e940eada959bdde927426e2531c9Paul Duffin int y = maxLog10ForLeadingZeros[Long.numberOfLeadingZeros(x)]; 1617dd252788645e940eada959bdde927426e2531c9Paul Duffin /* 1620888a09821a98ac0680fad765217302858e70fa4Paul Duffin * y is the higher of the two possible values of floor(log10(x)). If x < 10^y, then we want the 1630888a09821a98ac0680fad765217302858e70fa4Paul Duffin * lower of the two possible values, or y - 1, otherwise, we want y. 1647dd252788645e940eada959bdde927426e2531c9Paul Duffin */ 1650888a09821a98ac0680fad765217302858e70fa4Paul Duffin return y - lessThanBranchFree(x, powersOf10[y]); 1661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 1671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1687dd252788645e940eada959bdde927426e2531c9Paul Duffin // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i))) 1690888a09821a98ac0680fad765217302858e70fa4Paul Duffin @VisibleForTesting static final byte[] maxLog10ForLeadingZeros = { 1700888a09821a98ac0680fad765217302858e70fa4Paul Duffin 19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 1710888a09821a98ac0680fad765217302858e70fa4Paul Duffin 12, 12, 11, 11, 11, 10, 10, 10, 9, 9, 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4, 1720888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, 0 }; 1737dd252788645e940eada959bdde927426e2531c9Paul Duffin 1747dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 1751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @VisibleForTesting 1760888a09821a98ac0680fad765217302858e70fa4Paul Duffin static final long[] powersOf10 = { 1770888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L, 1780888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10L, 1790888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100L, 1800888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000L, 1810888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10000L, 1820888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100000L, 1830888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000000L, 1840888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10000000L, 1850888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100000000L, 1860888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000000000L, 1870888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10000000000L, 1880888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100000000000L, 1890888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000000000000L, 1900888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10000000000000L, 1910888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100000000000000L, 1920888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000000000000000L, 1930888a09821a98ac0680fad765217302858e70fa4Paul Duffin 10000000000000000L, 1940888a09821a98ac0680fad765217302858e70fa4Paul Duffin 100000000000000000L, 1950888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1000000000000000000L 1960888a09821a98ac0680fad765217302858e70fa4Paul Duffin }; 1977dd252788645e940eada959bdde927426e2531c9Paul Duffin 1987dd252788645e940eada959bdde927426e2531c9Paul Duffin // halfPowersOf10[i] = largest long less than 10^(i + 0.5) 1997dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 2007dd252788645e940eada959bdde927426e2531c9Paul Duffin @VisibleForTesting 2010888a09821a98ac0680fad765217302858e70fa4Paul Duffin static final long[] halfPowersOf10 = { 2020888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3L, 2030888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31L, 2040888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316L, 2050888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162L, 2060888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31622L, 2070888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316227L, 2080888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162277L, 2090888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31622776L, 2100888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316227766L, 2110888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162277660L, 2120888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31622776601L, 2130888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316227766016L, 2140888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162277660168L, 2150888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31622776601683L, 2160888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316227766016837L, 2170888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162277660168379L, 2180888a09821a98ac0680fad765217302858e70fa4Paul Duffin 31622776601683793L, 2190888a09821a98ac0680fad765217302858e70fa4Paul Duffin 316227766016837933L, 2200888a09821a98ac0680fad765217302858e70fa4Paul Duffin 3162277660168379331L 2210888a09821a98ac0680fad765217302858e70fa4Paul Duffin }; 2221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 2241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to 2251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)} 2261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * time. 2271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code k < 0} 2291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2307dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 2311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long pow(long b, int k) { 2321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("exponent", k); 2331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (-2 <= b && b <= 2) { 2341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch ((int) b) { 2351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 0: 2361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (k == 0) ? 1 : 0; 2371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 1: 2381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return 1; 2391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case (-1): 2401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return ((k & 1) == 0) ? 1 : -1; 2411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 2: 2421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (k < Long.SIZE) ? 1L << k : 0; 2431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case (-2): 2441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (k < Long.SIZE) { 2451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return ((k & 1) == 0) ? 1L << k : -(1L << k); 2461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } else { 2471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return 0; 2481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2497dd252788645e940eada959bdde927426e2531c9Paul Duffin default: 2507dd252788645e940eada959bdde927426e2531c9Paul Duffin throw new AssertionError(); 2511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert for (long accum = 1;; k >>= 1) { 2541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (k) { 2551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 0: 2561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return accum; 2571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 1: 2581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return accum * b; 2591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 2601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert accum *= ((k & 1) == 0) ? 1 : b; 2611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert b *= b; 2621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 2661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 2671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the square root of {@code x}, rounded with the specified rounding mode. 2681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 2691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code x < 0} 2701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and 2711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code sqrt(x)} is not an integer 2721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 2737dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 2741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 2751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long sqrt(long x, RoundingMode mode) { 2761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("x", x); 2771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (fitsInInt(x)) { 2781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return IntMath.sqrt((int) x, mode); 2791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 2800888a09821a98ac0680fad765217302858e70fa4Paul Duffin /* 2810888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Let k be the true value of floor(sqrt(x)), so that 2820888a09821a98ac0680fad765217302858e70fa4Paul Duffin * 2830888a09821a98ac0680fad765217302858e70fa4Paul Duffin * k * k <= x < (k + 1) * (k + 1) 2840888a09821a98ac0680fad765217302858e70fa4Paul Duffin * (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1)) 2850888a09821a98ac0680fad765217302858e70fa4Paul Duffin * since casting to double is nondecreasing. 2860888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Note that the right-hand inequality is no longer strict. 2870888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1)) 2880888a09821a98ac0680fad765217302858e70fa4Paul Duffin * since Math.sqrt is monotonic. 2890888a09821a98ac0680fad765217302858e70fa4Paul Duffin * (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1)) 2900888a09821a98ac0680fad765217302858e70fa4Paul Duffin * since casting to long is monotonic 2910888a09821a98ac0680fad765217302858e70fa4Paul Duffin * k <= (long) Math.sqrt(x) <= k + 1 2920888a09821a98ac0680fad765217302858e70fa4Paul Duffin * since (long) Math.sqrt(k * k) == k, as checked exhaustively in 2930888a09821a98ac0680fad765217302858e70fa4Paul Duffin * {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect} 2940888a09821a98ac0680fad765217302858e70fa4Paul Duffin */ 2950888a09821a98ac0680fad765217302858e70fa4Paul Duffin long guess = (long) Math.sqrt(x); 2960888a09821a98ac0680fad765217302858e70fa4Paul Duffin // Note: guess is always <= FLOOR_SQRT_MAX_LONG. 2970888a09821a98ac0680fad765217302858e70fa4Paul Duffin long guessSquared = guess * guess; 2980888a09821a98ac0680fad765217302858e70fa4Paul Duffin // Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is 2990888a09821a98ac0680fad765217302858e70fa4Paul Duffin // faster here than using lessThanBranchFree. 3001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 3011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 3020888a09821a98ac0680fad765217302858e70fa4Paul Duffin checkRoundingUnnecessary(guessSquared == x); 3030888a09821a98ac0680fad765217302858e70fa4Paul Duffin return guess; 3041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 3051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 3060888a09821a98ac0680fad765217302858e70fa4Paul Duffin if (x < guessSquared) { 3070888a09821a98ac0680fad765217302858e70fa4Paul Duffin return guess - 1; 3080888a09821a98ac0680fad765217302858e70fa4Paul Duffin } 3090888a09821a98ac0680fad765217302858e70fa4Paul Duffin return guess; 3101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 3111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 3120888a09821a98ac0680fad765217302858e70fa4Paul Duffin if (x > guessSquared) { 3130888a09821a98ac0680fad765217302858e70fa4Paul Duffin return guess + 1; 3140888a09821a98ac0680fad765217302858e70fa4Paul Duffin } 3150888a09821a98ac0680fad765217302858e70fa4Paul Duffin return guess; 3161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 3171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 3180888a09821a98ac0680fad765217302858e70fa4Paul Duffin case HALF_EVEN: 3190888a09821a98ac0680fad765217302858e70fa4Paul Duffin long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0); 3201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor; 3211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 3221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both 3231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * x and halfSquare are integers, this is equivalent to testing whether or not x <= 3241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * halfSquare. (We have to deal with overflow, though.) 3250888a09821a98ac0680fad765217302858e70fa4Paul Duffin * 3260888a09821a98ac0680fad765217302858e70fa4Paul Duffin * If we treat halfSquare as an unsigned long, we know that 3270888a09821a98ac0680fad765217302858e70fa4Paul Duffin * sqrtFloor^2 <= x < (sqrtFloor + 1)^2 3280888a09821a98ac0680fad765217302858e70fa4Paul Duffin * halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1 3290888a09821a98ac0680fad765217302858e70fa4Paul Duffin * so |x - halfSquare| <= sqrtFloor. Therefore, it's safe to treat x - halfSquare as a 3300888a09821a98ac0680fad765217302858e70fa4Paul Duffin * signed long, so lessThanBranchFree is safe for use. 3311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 3320888a09821a98ac0680fad765217302858e70fa4Paul Duffin return sqrtFloor + lessThanBranchFree(halfSquare, x); 3331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 3341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 3351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 3391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the result of dividing {@code p} by {@code q}, rounding using the specified 3401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code RoundingMode}. 3411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 3421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a} 3431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * is not an integer multiple of {@code b} 3441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 3457dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 3461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert @SuppressWarnings("fallthrough") 3471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long divide(long p, long q, RoundingMode mode) { 3481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNotNull(mode); 3491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long div = p / q; // throws if q == 0 3501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long rem = p - q * div; // equals p % q 3511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (rem == 0) { 3531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return div; 3541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 3561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 3571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to 3581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of 3591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * p / q. 3601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 3611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise. 3621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 3631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1)); 3641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert boolean increment; 3651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (mode) { 3661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UNNECESSARY: 3671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkRoundingUnnecessary(rem == 0); 3681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // fall through 3691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case DOWN: 3701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = false; 3711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert break; 3721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case UP: 3731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = true; 3741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert break; 3751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case CEILING: 3761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = signum > 0; 3771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert break; 3781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case FLOOR: 3791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = signum < 0; 3801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert break; 3811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_EVEN: 3821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_DOWN: 3831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case HALF_UP: 3841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long absRem = abs(rem); 3851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long cmpRemToHalfDivisor = absRem - (abs(q) - absRem); 3861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // subtracting two nonnegative longs can't overflow 3871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2). 3881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (cmpRemToHalfDivisor == 0) { // exactly on the half mark 3891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = (mode == HALF_UP | (mode == HALF_EVEN & (div & 1) != 0)); 3901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } else { 3911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert increment = cmpRemToHalfDivisor > 0; // closer to the UP value 3921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert break; 3941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 3951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert throw new AssertionError(); 3961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return increment ? div + signum : div; 3981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 3991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 4001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 4010888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Returns {@code x mod m}, a non-negative value less than {@code m}. 4020888a09821a98ac0680fad765217302858e70fa4Paul Duffin * This differs from {@code x % m}, which might be negative. 4031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>For example: 4051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <pre> {@code 4071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(7, 4) == 3 4091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-7, 4) == 1 4101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-1, 4) == 3 4111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-8, 4) == 0 4121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(8, 4) == 0}</pre> 4131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code m <= 0} 4150888a09821a98ac0680fad765217302858e70fa4Paul Duffin * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3"> 4160888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Remainder Operator</a> 4171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 4187dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 4191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static int mod(long x, int m) { 4201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Cast is safe because the result is guaranteed in the range [0, m) 4211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (int) mod(x, (long) m); 4221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 4231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 4241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 4250888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Returns {@code x mod m}, a non-negative value less than {@code m}. 4260888a09821a98ac0680fad765217302858e70fa4Paul Duffin * This differs from {@code x % m}, which might be negative. 4271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <p>For example: 4291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * <pre> {@code 4311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(7, 4) == 3 4331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-7, 4) == 1 4341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-1, 4) == 3 4351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(-8, 4) == 0 4361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * mod(8, 4) == 0}</pre> 4371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code m <= 0} 4390888a09821a98ac0680fad765217302858e70fa4Paul Duffin * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3"> 4400888a09821a98ac0680fad765217302858e70fa4Paul Duffin * Remainder Operator</a> 4411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 4427dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 4431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long mod(long x, long m) { 4441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (m <= 0) { 4450888a09821a98ac0680fad765217302858e70fa4Paul Duffin throw new ArithmeticException("Modulus must be positive"); 4461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 4471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long result = x % m; 4481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (result >= 0) ? result : result + m; 4491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 4501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 4511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 4521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if 4531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code a == 0 && b == 0}. 4541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 4551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0} 4561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 4571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long gcd(long a, long b) { 4581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 4591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * The reason we require both arguments to be >= 0 is because otherwise, what do you return on 4601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't 4611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * an int. 4621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 4631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("a", a); 4641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("b", b); 4657dd252788645e940eada959bdde927426e2531c9Paul Duffin if (a == 0) { 4667dd252788645e940eada959bdde927426e2531c9Paul Duffin // 0 % b == 0, so b divides a, but the converse doesn't hold. 4677dd252788645e940eada959bdde927426e2531c9Paul Duffin // BigInteger.gcd is consistent with this decision. 4687dd252788645e940eada959bdde927426e2531c9Paul Duffin return b; 4697dd252788645e940eada959bdde927426e2531c9Paul Duffin } else if (b == 0) { 4707dd252788645e940eada959bdde927426e2531c9Paul Duffin return a; // similar logic 4711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 4721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 4731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. 4747dd252788645e940eada959bdde927426e2531c9Paul Duffin * This is >60% faster than the Euclidean algorithm in benchmarks. 4751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 4761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int aTwos = Long.numberOfTrailingZeros(a); 4771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert a >>= aTwos; // divide out all 2s 4781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int bTwos = Long.numberOfTrailingZeros(b); 4791d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert b >>= bTwos; // divide out all 2s 4801d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert while (a != b) { // both a, b are odd 4817dd252788645e940eada959bdde927426e2531c9Paul Duffin // The key to the binary GCD algorithm is as follows: 4827dd252788645e940eada959bdde927426e2531c9Paul Duffin // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b). 4837dd252788645e940eada959bdde927426e2531c9Paul Duffin // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two. 4847dd252788645e940eada959bdde927426e2531c9Paul Duffin 4857dd252788645e940eada959bdde927426e2531c9Paul Duffin // We bend over backwards to avoid branching, adapting a technique from 4867dd252788645e940eada959bdde927426e2531c9Paul Duffin // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax 4877dd252788645e940eada959bdde927426e2531c9Paul Duffin 4887dd252788645e940eada959bdde927426e2531c9Paul Duffin long delta = a - b; // can't overflow, since a and b are nonnegative 4897dd252788645e940eada959bdde927426e2531c9Paul Duffin 4907dd252788645e940eada959bdde927426e2531c9Paul Duffin long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1)); 4917dd252788645e940eada959bdde927426e2531c9Paul Duffin // equivalent to Math.min(delta, 0) 4927dd252788645e940eada959bdde927426e2531c9Paul Duffin 4937dd252788645e940eada959bdde927426e2531c9Paul Duffin a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b) 4947dd252788645e940eada959bdde927426e2531c9Paul Duffin // a is now nonnegative and even 4957dd252788645e940eada959bdde927426e2531c9Paul Duffin 4967dd252788645e940eada959bdde927426e2531c9Paul Duffin b += minDeltaOrZero; // sets b to min(old a, b) 4971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b 4981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 4991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return a << min(aTwos, bTwos); 5001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 5021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 5031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the sum of {@code a} and {@code b}, provided it does not overflow. 5041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic 5061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 5077dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 5081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long checkedAdd(long a, long b) { 5091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long result = a + b; 5101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0); 5111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return result; 5121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 5141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 5151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the difference of {@code a} and {@code b}, provided it does not overflow. 5161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic 5181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 5197dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 5201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long checkedSubtract(long a, long b) { 5211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long result = a - b; 5221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0); 5231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return result; 5241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 5261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 5271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the product of {@code a} and {@code b}, provided it does not overflow. 5281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic 5301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 5317dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 5321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long checkedMultiply(long a, long b) { 5331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // Hacker's Delight, Section 2-12 5341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert int leadingZeros = Long.numberOfLeadingZeros(a) + Long.numberOfLeadingZeros(~a) 5351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert + Long.numberOfLeadingZeros(b) + Long.numberOfLeadingZeros(~b); 5361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 5371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely 5381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * bad. We do the leadingZeros check to avoid the division below if at all possible. 5391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take 5411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * care of all a < 0 with their own check, because in particular, the case a == -1 will 5421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * incorrectly pass the division check below. 5431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * In all other cases, we check that either a is 0 or the result is consistent with division. 5451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 5461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (leadingZeros > Long.SIZE + 1) { 5471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return a * b; 5481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(leadingZeros >= Long.SIZE); 5501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(a >= 0 | b != Long.MIN_VALUE); 5511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long result = a * b; 5521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(a == 0 || result / a == b); 5531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return result; 5541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 5561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 5571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns the {@code b} to the {@code k}th power, provided it does not overflow. 5581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 5591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed 5601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code long} arithmetic 5611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 5627dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 5631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long checkedPow(long b, int k) { 5641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("exponent", k); 5651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (b >= -2 & b <= 2) { 5661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch ((int) b) { 5671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 0: 5681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (k == 0) ? 1 : 0; 5691d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 1: 5701d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return 1; 5711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case (-1): 5721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return ((k & 1) == 0) ? 1 : -1; 5731d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 2: 5741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(k < Long.SIZE - 1); 5751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return 1L << k; 5761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case (-2): 5771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(k < Long.SIZE); 5781d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return ((k & 1) == 0) ? (1L << k) : (-1L << k); 5797dd252788645e940eada959bdde927426e2531c9Paul Duffin default: 5807dd252788645e940eada959bdde927426e2531c9Paul Duffin throw new AssertionError(); 5811d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5821d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5831d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert long accum = 1; 5841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert while (true) { 5851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert switch (k) { 5861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 0: 5871d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return accum; 5881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert case 1: 5891d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return checkedMultiply(accum, b); 5901d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert default: 5911d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if ((k & 1) != 0) { 5921d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert accum = checkedMultiply(accum, b); 5931d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5941d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert k >>= 1; 5951d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (k > 0) { 5961d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNoOverflow(b <= FLOOR_SQRT_MAX_LONG); 5971d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert b *= b; 5981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 5991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 6001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 6011d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 6021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 6030888a09821a98ac0680fad765217302858e70fa4Paul Duffin @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L; 6041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 6051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 6061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code n!}, that is, the product of the first {@code n} positive 6071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * integers, {@code 1} if {@code n == 0}, or {@link Long#MAX_VALUE} if the 6081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * result does not fit in a {@code long}. 6091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 6101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code n < 0} 6111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 6127dd252788645e940eada959bdde927426e2531c9Paul Duffin @GwtIncompatible("TODO") 6131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long factorial(int n) { 6141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("n", n); 6157dd252788645e940eada959bdde927426e2531c9Paul Duffin return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE; 6161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 6171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 6180888a09821a98ac0680fad765217302858e70fa4Paul Duffin static final long[] factorials = { 6190888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L, 6200888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L, 6210888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2, 6220888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3, 6230888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4, 6240888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5, 6250888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6, 6260888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7, 6270888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8, 6280888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9, 6290888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10, 6300888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11, 6311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12, 6321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13, 6331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14, 6341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15, 6351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16, 6361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17, 6371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18, 6381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19, 6390888a09821a98ac0680fad765217302858e70fa4Paul Duffin 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 6400888a09821a98ac0680fad765217302858e70fa4Paul Duffin }; 6411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 6421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /** 6431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and 6441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}. 6451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * 6461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n} 6471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 6481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert public static long binomial(int n, int k) { 6491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("n", n); 6501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkNonNegative("k", k); 6511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert checkArgument(k <= n, "k (%s) > n (%s)", k, n); 6521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert if (k > (n >> 1)) { 6531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert k = n - k; 6541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 6557dd252788645e940eada959bdde927426e2531c9Paul Duffin switch (k) { 6567dd252788645e940eada959bdde927426e2531c9Paul Duffin case 0: 6577dd252788645e940eada959bdde927426e2531c9Paul Duffin return 1; 6587dd252788645e940eada959bdde927426e2531c9Paul Duffin case 1: 6597dd252788645e940eada959bdde927426e2531c9Paul Duffin return n; 6607dd252788645e940eada959bdde927426e2531c9Paul Duffin default: 6617dd252788645e940eada959bdde927426e2531c9Paul Duffin if (n < factorials.length) { 6627dd252788645e940eada959bdde927426e2531c9Paul Duffin return factorials[n] / (factorials[k] * factorials[n - k]); 6637dd252788645e940eada959bdde927426e2531c9Paul Duffin } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) { 6647dd252788645e940eada959bdde927426e2531c9Paul Duffin return Long.MAX_VALUE; 6657dd252788645e940eada959bdde927426e2531c9Paul Duffin } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) { 6667dd252788645e940eada959bdde927426e2531c9Paul Duffin // guaranteed not to overflow 6677dd252788645e940eada959bdde927426e2531c9Paul Duffin long result = n--; 6687dd252788645e940eada959bdde927426e2531c9Paul Duffin for (int i = 2; i <= k; n--, i++) { 6697dd252788645e940eada959bdde927426e2531c9Paul Duffin result *= n; 6707dd252788645e940eada959bdde927426e2531c9Paul Duffin result /= i; 6717dd252788645e940eada959bdde927426e2531c9Paul Duffin } 6727dd252788645e940eada959bdde927426e2531c9Paul Duffin return result; 6737dd252788645e940eada959bdde927426e2531c9Paul Duffin } else { 6747dd252788645e940eada959bdde927426e2531c9Paul Duffin int nBits = LongMath.log2(n, RoundingMode.CEILING); 6757dd252788645e940eada959bdde927426e2531c9Paul Duffin 6767dd252788645e940eada959bdde927426e2531c9Paul Duffin long result = 1; 6777dd252788645e940eada959bdde927426e2531c9Paul Duffin long numerator = n--; 6787dd252788645e940eada959bdde927426e2531c9Paul Duffin long denominator = 1; 6797dd252788645e940eada959bdde927426e2531c9Paul Duffin 6807dd252788645e940eada959bdde927426e2531c9Paul Duffin int numeratorBits = nBits; 6817dd252788645e940eada959bdde927426e2531c9Paul Duffin // This is an upper bound on log2(numerator, ceiling). 6827dd252788645e940eada959bdde927426e2531c9Paul Duffin 6837dd252788645e940eada959bdde927426e2531c9Paul Duffin /* 6847dd252788645e940eada959bdde927426e2531c9Paul Duffin * We want to do this in long math for speed, but want to avoid overflow. We adapt the 6857dd252788645e940eada959bdde927426e2531c9Paul Duffin * technique previously used by BigIntegerMath: maintain separate numerator and 6867dd252788645e940eada959bdde927426e2531c9Paul Duffin * denominator accumulators, multiplying the fraction into result when near overflow. 6877dd252788645e940eada959bdde927426e2531c9Paul Duffin */ 6887dd252788645e940eada959bdde927426e2531c9Paul Duffin for (int i = 2; i <= k; i++, n--) { 6897dd252788645e940eada959bdde927426e2531c9Paul Duffin if (numeratorBits + nBits < Long.SIZE - 1) { 6907dd252788645e940eada959bdde927426e2531c9Paul Duffin // It's definitely safe to multiply into numerator and denominator. 6917dd252788645e940eada959bdde927426e2531c9Paul Duffin numerator *= n; 6927dd252788645e940eada959bdde927426e2531c9Paul Duffin denominator *= i; 6937dd252788645e940eada959bdde927426e2531c9Paul Duffin numeratorBits += nBits; 6947dd252788645e940eada959bdde927426e2531c9Paul Duffin } else { 6957dd252788645e940eada959bdde927426e2531c9Paul Duffin // It might not be safe to multiply into numerator and denominator, 6967dd252788645e940eada959bdde927426e2531c9Paul Duffin // so multiply (numerator / denominator) into result. 6977dd252788645e940eada959bdde927426e2531c9Paul Duffin result = multiplyFraction(result, numerator, denominator); 6987dd252788645e940eada959bdde927426e2531c9Paul Duffin numerator = n; 6997dd252788645e940eada959bdde927426e2531c9Paul Duffin denominator = i; 7007dd252788645e940eada959bdde927426e2531c9Paul Duffin numeratorBits = nBits; 7017dd252788645e940eada959bdde927426e2531c9Paul Duffin } 7027dd252788645e940eada959bdde927426e2531c9Paul Duffin } 7037dd252788645e940eada959bdde927426e2531c9Paul Duffin return multiplyFraction(result, numerator, denominator); 7047dd252788645e940eada959bdde927426e2531c9Paul Duffin } 7051d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 7067dd252788645e940eada959bdde927426e2531c9Paul Duffin } 7077dd252788645e940eada959bdde927426e2531c9Paul Duffin 7087dd252788645e940eada959bdde927426e2531c9Paul Duffin /** 7097dd252788645e940eada959bdde927426e2531c9Paul Duffin * Returns (x * numerator / denominator), which is assumed to come out to an integral value. 7107dd252788645e940eada959bdde927426e2531c9Paul Duffin */ 7117dd252788645e940eada959bdde927426e2531c9Paul Duffin static long multiplyFraction(long x, long numerator, long denominator) { 7127dd252788645e940eada959bdde927426e2531c9Paul Duffin if (x == 1) { 7137dd252788645e940eada959bdde927426e2531c9Paul Duffin return numerator / denominator; 7141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 7157dd252788645e940eada959bdde927426e2531c9Paul Duffin long commonDivisor = gcd(x, denominator); 7167dd252788645e940eada959bdde927426e2531c9Paul Duffin x /= commonDivisor; 7177dd252788645e940eada959bdde927426e2531c9Paul Duffin denominator /= commonDivisor; 7187dd252788645e940eada959bdde927426e2531c9Paul Duffin // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact, 7197dd252788645e940eada959bdde927426e2531c9Paul Duffin // so denominator must be a divisor of numerator. 7207dd252788645e940eada959bdde927426e2531c9Paul Duffin return x * (numerator / denominator); 7211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 7221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 7231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 7247dd252788645e940eada959bdde927426e2531c9Paul Duffin * binomial(biggestBinomials[k], k) fits in a long, but not 7257dd252788645e940eada959bdde927426e2531c9Paul Duffin * binomial(biggestBinomials[k] + 1, k). 7261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 7270888a09821a98ac0680fad765217302858e70fa4Paul Duffin static final int[] biggestBinomials = 7280888a09821a98ac0680fad765217302858e70fa4Paul Duffin {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 3810779, 121977, 16175, 4337, 1733, 7290888a09821a98ac0680fad765217302858e70fa4Paul Duffin 887, 534, 361, 265, 206, 169, 143, 125, 111, 101, 94, 88, 83, 79, 76, 74, 72, 70, 69, 68, 7300888a09821a98ac0680fad765217302858e70fa4Paul Duffin 67, 67, 66, 66, 66, 66}; 7311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 7321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert /* 7337dd252788645e940eada959bdde927426e2531c9Paul Duffin * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, 7347dd252788645e940eada959bdde927426e2531c9Paul Duffin * but binomial(biggestSimpleBinomials[k] + 1, k) does. 7351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */ 7360888a09821a98ac0680fad765217302858e70fa4Paul Duffin @VisibleForTesting static final int[] biggestSimpleBinomials = 7370888a09821a98ac0680fad765217302858e70fa4Paul Duffin {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2642246, 86251, 11724, 3218, 1313, 7380888a09821a98ac0680fad765217302858e70fa4Paul Duffin 684, 419, 287, 214, 169, 139, 119, 105, 95, 87, 81, 76, 73, 70, 68, 66, 64, 63, 62, 62, 7390888a09821a98ac0680fad765217302858e70fa4Paul Duffin 61, 61, 61}; 7401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // These values were generated by using checkedMultiply to see when the simple multiply/divide 7411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert // algorithm would lead to an overflow. 7421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 7431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert static boolean fitsInInt(long x) { 7441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert return (int) x == x; 7451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert } 7461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert 7477dd252788645e940eada959bdde927426e2531c9Paul Duffin /** 7487dd252788645e940eada959bdde927426e2531c9Paul Duffin * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward 7497dd252788645e940eada959bdde927426e2531c9Paul Duffin * negative infinity. This method is resilient to overflow. 7507dd252788645e940eada959bdde927426e2531c9Paul Duffin * 7517dd252788645e940eada959bdde927426e2531c9Paul Duffin * @since 14.0 7527dd252788645e940eada959bdde927426e2531c9Paul Duffin */ 7537dd252788645e940eada959bdde927426e2531c9Paul Duffin public static long mean(long x, long y) { 7547dd252788645e940eada959bdde927426e2531c9Paul Duffin // Efficient method for computing the arithmetic mean. 7557dd252788645e940eada959bdde927426e2531c9Paul Duffin // The alternative (x + y) / 2 fails for large values. 7567dd252788645e940eada959bdde927426e2531c9Paul Duffin // The alternative (x + y) >>> 1 fails for negative values. 7577dd252788645e940eada959bdde927426e2531c9Paul Duffin return (x & y) + ((x ^ y) >> 1); 7587dd252788645e940eada959bdde927426e2531c9Paul Duffin } 7597dd252788645e940eada959bdde927426e2531c9Paul Duffin 7601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert private LongMath() {} 7611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert} 762