com.google.common.math.LongMath Maven / Gradle / Ivy
/*
* Copyright (C) 2011 The Guava Authors
*
* Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
* in compliance with the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software distributed under the License
* is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
* or implied. See the License for the specific language governing permissions and limitations under
* the License.
*/
package com.google.common.math;
import static com.google.common.base.Preconditions.checkArgument;
import static com.google.common.base.Preconditions.checkNotNull;
import static com.google.common.math.MathPreconditions.checkNoOverflow;
import static com.google.common.math.MathPreconditions.checkNonNegative;
import static com.google.common.math.MathPreconditions.checkPositive;
import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
import static java.lang.Math.abs;
import static java.lang.Math.min;
import static java.math.RoundingMode.HALF_EVEN;
import static java.math.RoundingMode.HALF_UP;
import com.google.common.annotations.Beta;
import com.google.common.annotations.GwtCompatible;
import com.google.common.annotations.GwtIncompatible;
import com.google.common.annotations.VisibleForTesting;
import com.google.common.primitives.UnsignedLongs;
import java.math.BigInteger;
import java.math.RoundingMode;
/**
* A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
* named analogously to their {@code BigInteger} counterparts.
*
* The implementations of many methods in this class are based on material from Henry S. Warren,
* Jr.'s Hacker's Delight, (Addison Wesley, 2002).
*
*
Similar functionality for {@code int} and for {@link BigInteger} can be found in {@link
* IntMath} and {@link BigIntegerMath} respectively. For other common operations on {@code long}
* values, see {@link com.google.common.primitives.Longs}.
*
* @author Louis Wasserman
* @since 11.0
*/
@GwtCompatible(emulated = true)
public final class LongMath {
// NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
@VisibleForTesting static final long MAX_SIGNED_POWER_OF_TWO = 1L << (Long.SIZE - 2);
/**
* Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to
* {@code checkedPow(2, log2(x, CEILING))}.
*
* @throws IllegalArgumentException if {@code x <= 0}
* @throws ArithmeticException of the next-higher power of two is not representable as a {@code
* long}, i.e. when {@code x > 2^62}
* @since 20.0
*/
@Beta
public static long ceilingPowerOfTwo(long x) {
checkPositive("x", x);
if (x > MAX_SIGNED_POWER_OF_TWO) {
throw new ArithmeticException("ceilingPowerOfTwo(" + x + ") is not representable as a long");
}
return 1L << -Long.numberOfLeadingZeros(x - 1);
}
/**
* Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code
* checkedPow(2, log2(x, FLOOR))}.
*
* @throws IllegalArgumentException if {@code x <= 0}
* @since 20.0
*/
@Beta
public static long floorPowerOfTwo(long x) {
checkPositive("x", x);
// Long.highestOneBit was buggy on GWT. We've fixed it, but I'm not certain when the fix will
// be released.
return 1L << ((Long.SIZE - 1) - Long.numberOfLeadingZeros(x));
}
/**
* Returns {@code true} if {@code x} represents a power of two.
*
*
This differs from {@code Long.bitCount(x) == 1}, because {@code
* Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
*/
public static boolean isPowerOfTwo(long x) {
return x > 0 & (x & (x - 1)) == 0;
}
/**
* Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
* signed long. The implementation is branch-free, and benchmarks suggest it is measurably faster
* than the straightforward ternary expression.
*/
@VisibleForTesting
static int lessThanBranchFree(long x, long y) {
// Returns the sign bit of x - y.
return (int) (~~(x - y) >>> (Long.SIZE - 1));
}
/**
* Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
*
* @throws IllegalArgumentException if {@code x <= 0}
* @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
* is not a power of two
*/
@SuppressWarnings("fallthrough")
// TODO(kevinb): remove after this warning is disabled globally
public static int log2(long x, RoundingMode mode) {
checkPositive("x", x);
switch (mode) {
case UNNECESSARY:
checkRoundingUnnecessary(isPowerOfTwo(x));
// fall through
case DOWN:
case FLOOR:
return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
case UP:
case CEILING:
return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
case HALF_DOWN:
case HALF_UP:
case HALF_EVEN:
// Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
int leadingZeros = Long.numberOfLeadingZeros(x);
long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
// floor(2^(logFloor + 0.5))
int logFloor = (Long.SIZE - 1) - leadingZeros;
return logFloor + lessThanBranchFree(cmp, x);
default:
throw new AssertionError("impossible");
}
}
/** The biggest half power of two that fits into an unsigned long */
@VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
/**
* Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
*
* @throws IllegalArgumentException if {@code x <= 0}
* @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
* is not a power of ten
*/
@GwtIncompatible // TODO
@SuppressWarnings("fallthrough")
// TODO(kevinb): remove after this warning is disabled globally
public static int log10(long x, RoundingMode mode) {
checkPositive("x", x);
int logFloor = log10Floor(x);
long floorPow = powersOf10[logFloor];
switch (mode) {
case UNNECESSARY:
checkRoundingUnnecessary(x == floorPow);
// fall through
case FLOOR:
case DOWN:
return logFloor;
case CEILING:
case UP:
return logFloor + lessThanBranchFree(floorPow, x);
case HALF_DOWN:
case HALF_UP:
case HALF_EVEN:
// sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x);
default:
throw new AssertionError();
}
}
@GwtIncompatible // TODO
static int log10Floor(long x) {
/*
* Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation.
*
* The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))), we
* can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x)) is 6,
* then 64 <= x < 128, so floor(log10(x)) is either 1 or 2.
*/
int y = maxLog10ForLeadingZeros[Long.numberOfLeadingZeros(x)];
/*
* y is the higher of the two possible values of floor(log10(x)). If x < 10^y, then we want the
* lower of the two possible values, or y - 1, otherwise, we want y.
*/
return y - lessThanBranchFree(x, powersOf10[y]);
}
// maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i)))
@VisibleForTesting
static final byte[] maxLog10ForLeadingZeros = {
19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 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, 3, 3, 3,
3, 2, 2, 2, 1, 1, 1, 0, 0, 0
};
@GwtIncompatible // TODO
@VisibleForTesting
static final long[] powersOf10 = {
1L,
10L,
100L,
1000L,
10000L,
100000L,
1000000L,
10000000L,
100000000L,
1000000000L,
10000000000L,
100000000000L,
1000000000000L,
10000000000000L,
100000000000000L,
1000000000000000L,
10000000000000000L,
100000000000000000L,
1000000000000000000L
};
// halfPowersOf10[i] = largest long less than 10^(i + 0.5)
@GwtIncompatible // TODO
@VisibleForTesting
static final long[] halfPowersOf10 = {
3L,
31L,
316L,
3162L,
31622L,
316227L,
3162277L,
31622776L,
316227766L,
3162277660L,
31622776601L,
316227766016L,
3162277660168L,
31622776601683L,
316227766016837L,
3162277660168379L,
31622776601683793L,
316227766016837933L,
3162277660168379331L
};
/**
* Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
* {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
* time.
*
* @throws IllegalArgumentException if {@code k < 0}
*/
@GwtIncompatible // TODO
public static long pow(long b, int k) {
checkNonNegative("exponent", k);
if (-2 <= b && b <= 2) {
switch ((int) b) {
case 0:
return (k == 0) ? 1 : 0;
case 1:
return 1;
case (-1):
return ((k & 1) == 0) ? 1 : -1;
case 2:
return (k < Long.SIZE) ? 1L << k : 0;
case (-2):
if (k < Long.SIZE) {
return ((k & 1) == 0) ? 1L << k : -(1L << k);
} else {
return 0;
}
default:
throw new AssertionError();
}
}
for (long accum = 1; ; k >>= 1) {
switch (k) {
case 0:
return accum;
case 1:
return accum * b;
default:
accum *= ((k & 1) == 0) ? 1 : b;
b *= b;
}
}
}
/**
* Returns the square root of {@code x}, rounded with the specified rounding mode.
*
* @throws IllegalArgumentException if {@code x < 0}
* @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
* sqrt(x)} is not an integer
*/
@GwtIncompatible // TODO
@SuppressWarnings("fallthrough")
public static long sqrt(long x, RoundingMode mode) {
checkNonNegative("x", x);
if (fitsInInt(x)) {
return IntMath.sqrt((int) x, mode);
}
/*
* Let k be the true value of floor(sqrt(x)), so that
*
* k * k <= x < (k + 1) * (k + 1)
* (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1))
* since casting to double is nondecreasing.
* Note that the right-hand inequality is no longer strict.
* Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1))
* since Math.sqrt is monotonic.
* (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1))
* since casting to long is monotonic
* k <= (long) Math.sqrt(x) <= k + 1
* since (long) Math.sqrt(k * k) == k, as checked exhaustively in
* {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect}
*/
long guess = (long) Math.sqrt(x);
// Note: guess is always <= FLOOR_SQRT_MAX_LONG.
long guessSquared = guess * guess;
// Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is
// faster here than using lessThanBranchFree.
switch (mode) {
case UNNECESSARY:
checkRoundingUnnecessary(guessSquared == x);
return guess;
case FLOOR:
case DOWN:
if (x < guessSquared) {
return guess - 1;
}
return guess;
case CEILING:
case UP:
if (x > guessSquared) {
return guess + 1;
}
return guess;
case HALF_DOWN:
case HALF_UP:
case HALF_EVEN:
long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0);
long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
/*
* We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
* and halfSquare are integers, this is equivalent to testing whether or not x <=
* halfSquare. (We have to deal with overflow, though.)
*
* If we treat halfSquare as an unsigned long, we know that
* sqrtFloor^2 <= x < (sqrtFloor + 1)^2
* halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1
* so |x - halfSquare| <= sqrtFloor. Therefore, it's safe to treat x - halfSquare as a
* signed long, so lessThanBranchFree is safe for use.
*/
return sqrtFloor + lessThanBranchFree(halfSquare, x);
default:
throw new AssertionError();
}
}
/**
* Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
* RoundingMode}.
*
* @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
* is not an integer multiple of {@code b}
*/
@GwtIncompatible // TODO
@SuppressWarnings("fallthrough")
public static long divide(long p, long q, RoundingMode mode) {
checkNotNull(mode);
long div = p / q; // throws if q == 0
long rem = p - q * div; // equals p % q
if (rem == 0) {
return div;
}
/*
* Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
* deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
* p / q.
*
* signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
*/
int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
boolean increment;
switch (mode) {
case UNNECESSARY:
checkRoundingUnnecessary(rem == 0);
// fall through
case DOWN:
increment = false;
break;
case UP:
increment = true;
break;
case CEILING:
increment = signum > 0;
break;
case FLOOR:
increment = signum < 0;
break;
case HALF_EVEN:
case HALF_DOWN:
case HALF_UP:
long absRem = abs(rem);
long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
// subtracting two nonnegative longs can't overflow
// cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
increment = (mode == HALF_UP | (mode == HALF_EVEN & (div & 1) != 0));
} else {
increment = cmpRemToHalfDivisor > 0; // closer to the UP value
}
break;
default:
throw new AssertionError();
}
return increment ? div + signum : div;
}
/**
* Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
* m}, which might be negative.
*
*
For example:
*
*
{@code
* mod(7, 4) == 3
* mod(-7, 4) == 1
* mod(-1, 4) == 3
* mod(-8, 4) == 0
* mod(8, 4) == 0
* }
*
* @throws ArithmeticException if {@code m <= 0}
* @see
* Remainder Operator
*/
@GwtIncompatible // TODO
public static int mod(long x, int m) {
// Cast is safe because the result is guaranteed in the range [0, m)
return (int) mod(x, (long) m);
}
/**
* Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
* m}, which might be negative.
*
* For example:
*
*
{@code
* mod(7, 4) == 3
* mod(-7, 4) == 1
* mod(-1, 4) == 3
* mod(-8, 4) == 0
* mod(8, 4) == 0
* }
*
* @throws ArithmeticException if {@code m <= 0}
* @see
* Remainder Operator
*/
@GwtIncompatible // TODO
public static long mod(long x, long m) {
if (m <= 0) {
throw new ArithmeticException("Modulus must be positive");
}
long result = x % m;
return (result >= 0) ? result : result + m;
}
/**
* Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if {@code a == 0 && b ==
* 0}.
*
* @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
*/
public static long gcd(long a, long b) {
/*
* The reason we require both arguments to be >= 0 is because otherwise, what do you return on
* gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't an
* int.
*/
checkNonNegative("a", a);
checkNonNegative("b", b);
if (a == 0) {
// 0 % b == 0, so b divides a, but the converse doesn't hold.
// BigInteger.gcd is consistent with this decision.
return b;
} else if (b == 0) {
return a; // similar logic
}
/*
* Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. This is
* >60% faster than the Euclidean algorithm in benchmarks.
*/
int aTwos = Long.numberOfTrailingZeros(a);
a >>= aTwos; // divide out all 2s
int bTwos = Long.numberOfTrailingZeros(b);
b >>= bTwos; // divide out all 2s
while (a != b) { // both a, b are odd
// The key to the binary GCD algorithm is as follows:
// Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
// But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
// We bend over backwards to avoid branching, adapting a technique from
// http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
long delta = a - b; // can't overflow, since a and b are nonnegative
long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
// equivalent to Math.min(delta, 0)
a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
// a is now nonnegative and even
b += minDeltaOrZero; // sets b to min(old a, b)
a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
}
return a << min(aTwos, bTwos);
}
/**
* Returns the sum of {@code a} and {@code b}, provided it does not overflow.
*
* @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
*/
@GwtIncompatible // TODO
public static long checkedAdd(long a, long b) {
long result = a + b;
checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0, "checkedAdd", a, b);
return result;
}
/**
* Returns the difference of {@code a} and {@code b}, provided it does not overflow.
*
* @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
*/
@GwtIncompatible // TODO
public static long checkedSubtract(long a, long b) {
long result = a - b;
checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0, "checkedSubtract", a, b);
return result;
}
/**
* Returns the product of {@code a} and {@code b}, provided it does not overflow.
*
* @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
*/
public static long checkedMultiply(long a, long b) {
// Hacker's Delight, Section 2-12
int leadingZeros =
Long.numberOfLeadingZeros(a)
+ Long.numberOfLeadingZeros(~a)
+ Long.numberOfLeadingZeros(b)
+ Long.numberOfLeadingZeros(~b);
/*
* If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
* bad. We do the leadingZeros check to avoid the division below if at all possible.
*
* Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
* care of all a < 0 with their own check, because in particular, the case a == -1 will
* incorrectly pass the division check below.
*
* In all other cases, we check that either a is 0 or the result is consistent with division.
*/
if (leadingZeros > Long.SIZE + 1) {
return a * b;
}
checkNoOverflow(leadingZeros >= Long.SIZE, "checkedMultiply", a, b);
checkNoOverflow(a >= 0 | b != Long.MIN_VALUE, "checkedMultiply", a, b);
long result = a * b;
checkNoOverflow(a == 0 || result / a == b, "checkedMultiply", a, b);
return result;
}
/**
* Returns the {@code b} to the {@code k}th power, provided it does not overflow.
*
* @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed {@code
* long} arithmetic
*/
@GwtIncompatible // TODO
public static long checkedPow(long b, int k) {
checkNonNegative("exponent", k);
if (b >= -2 & b <= 2) {
switch ((int) b) {
case 0:
return (k == 0) ? 1 : 0;
case 1:
return 1;
case (-1):
return ((k & 1) == 0) ? 1 : -1;
case 2:
checkNoOverflow(k < Long.SIZE - 1, "checkedPow", b, k);
return 1L << k;
case (-2):
checkNoOverflow(k < Long.SIZE, "checkedPow", b, k);
return ((k & 1) == 0) ? (1L << k) : (-1L << k);
default:
throw new AssertionError();
}
}
long accum = 1;
while (true) {
switch (k) {
case 0:
return accum;
case 1:
return checkedMultiply(accum, b);
default:
if ((k & 1) != 0) {
accum = checkedMultiply(accum, b);
}
k >>= 1;
if (k > 0) {
checkNoOverflow(
-FLOOR_SQRT_MAX_LONG <= b && b <= FLOOR_SQRT_MAX_LONG, "checkedPow", b, k);
b *= b;
}
}
}
}
/**
* Returns the sum of {@code a} and {@code b} unless it would overflow or underflow in which case
* {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
*
* @since 20.0
*/
@Beta
public static long saturatedAdd(long a, long b) {
long naiveSum = a + b;
if ((a ^ b) < 0 | (a ^ naiveSum) >= 0) {
// If a and b have different signs or a has the same sign as the result then there was no
// overflow, return.
return naiveSum;
}
// we did over/under flow, if the sign is negative we should return MAX otherwise MIN
return Long.MAX_VALUE + ((naiveSum >>> (Long.SIZE - 1)) ^ 1);
}
/**
* Returns the difference of {@code a} and {@code b} unless it would overflow or underflow in
* which case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
*
* @since 20.0
*/
@Beta
public static long saturatedSubtract(long a, long b) {
long naiveDifference = a - b;
if ((a ^ b) >= 0 | (a ^ naiveDifference) >= 0) {
// If a and b have the same signs or a has the same sign as the result then there was no
// overflow, return.
return naiveDifference;
}
// we did over/under flow
return Long.MAX_VALUE + ((naiveDifference >>> (Long.SIZE - 1)) ^ 1);
}
/**
* Returns the product of {@code a} and {@code b} unless it would overflow or underflow in which
* case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
*
* @since 20.0
*/
@Beta
public static long saturatedMultiply(long a, long b) {
// see checkedMultiply for explanation
int leadingZeros =
Long.numberOfLeadingZeros(a)
+ Long.numberOfLeadingZeros(~a)
+ Long.numberOfLeadingZeros(b)
+ Long.numberOfLeadingZeros(~b);
if (leadingZeros > Long.SIZE + 1) {
return a * b;
}
// the return value if we will overflow (which we calculate by overflowing a long :) )
long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
// overflow
return limit;
}
long result = a * b;
if (a == 0 || result / a == b) {
return result;
}
return limit;
}
/**
* Returns the {@code b} to the {@code k}th power, unless it would overflow or underflow in which
* case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
*
* @since 20.0
*/
@Beta
public static long saturatedPow(long b, int k) {
checkNonNegative("exponent", k);
if (b >= -2 & b <= 2) {
switch ((int) b) {
case 0:
return (k == 0) ? 1 : 0;
case 1:
return 1;
case (-1):
return ((k & 1) == 0) ? 1 : -1;
case 2:
if (k >= Long.SIZE - 1) {
return Long.MAX_VALUE;
}
return 1L << k;
case (-2):
if (k >= Long.SIZE) {
return Long.MAX_VALUE + (k & 1);
}
return ((k & 1) == 0) ? (1L << k) : (-1L << k);
default:
throw new AssertionError();
}
}
long accum = 1;
// if b is negative and k is odd then the limit is MIN otherwise the limit is MAX
long limit = Long.MAX_VALUE + ((b >>> Long.SIZE - 1) & (k & 1));
while (true) {
switch (k) {
case 0:
return accum;
case 1:
return saturatedMultiply(accum, b);
default:
if ((k & 1) != 0) {
accum = saturatedMultiply(accum, b);
}
k >>= 1;
if (k > 0) {
if (-FLOOR_SQRT_MAX_LONG > b | b > FLOOR_SQRT_MAX_LONG) {
return limit;
}
b *= b;
}
}
}
}
@VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
/**
* Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
* {@code n == 0}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
*
* @throws IllegalArgumentException if {@code n < 0}
*/
@GwtIncompatible // TODO
public static long factorial(int n) {
checkNonNegative("n", n);
return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
}
static final long[] factorials = {
1L,
1L,
1L * 2,
1L * 2 * 3,
1L * 2 * 3 * 4,
1L * 2 * 3 * 4 * 5,
1L * 2 * 3 * 4 * 5 * 6,
1L * 2 * 3 * 4 * 5 * 6 * 7,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
};
/**
* Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
* {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
*
* @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
*/
public static long binomial(int n, int k) {
checkNonNegative("n", n);
checkNonNegative("k", k);
checkArgument(k <= n, "k (%s) > n (%s)", k, n);
if (k > (n >> 1)) {
k = n - k;
}
switch (k) {
case 0:
return 1;
case 1:
return n;
default:
if (n < factorials.length) {
return factorials[n] / (factorials[k] * factorials[n - k]);
} else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
return Long.MAX_VALUE;
} else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
// guaranteed not to overflow
long result = n--;
for (int i = 2; i <= k; n--, i++) {
result *= n;
result /= i;
}
return result;
} else {
int nBits = LongMath.log2(n, RoundingMode.CEILING);
long result = 1;
long numerator = n--;
long denominator = 1;
int numeratorBits = nBits;
// This is an upper bound on log2(numerator, ceiling).
/*
* We want to do this in long math for speed, but want to avoid overflow. We adapt the
* technique previously used by BigIntegerMath: maintain separate numerator and
* denominator accumulators, multiplying the fraction into result when near overflow.
*/
for (int i = 2; i <= k; i++, n--) {
if (numeratorBits + nBits < Long.SIZE - 1) {
// It's definitely safe to multiply into numerator and denominator.
numerator *= n;
denominator *= i;
numeratorBits += nBits;
} else {
// It might not be safe to multiply into numerator and denominator,
// so multiply (numerator / denominator) into result.
result = multiplyFraction(result, numerator, denominator);
numerator = n;
denominator = i;
numeratorBits = nBits;
}
}
return multiplyFraction(result, numerator, denominator);
}
}
}
/** Returns (x * numerator / denominator), which is assumed to come out to an integral value. */
static long multiplyFraction(long x, long numerator, long denominator) {
if (x == 1) {
return numerator / denominator;
}
long commonDivisor = gcd(x, denominator);
x /= commonDivisor;
denominator /= commonDivisor;
// We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
// so denominator must be a divisor of numerator.
return x * (numerator / denominator);
}
/*
* binomial(biggestBinomials[k], k) fits in a long, but not binomial(biggestBinomials[k] + 1, k).
*/
static final int[] biggestBinomials = {
Integer.MAX_VALUE,
Integer.MAX_VALUE,
Integer.MAX_VALUE,
3810779,
121977,
16175,
4337,
1733,
887,
534,
361,
265,
206,
169,
143,
125,
111,
101,
94,
88,
83,
79,
76,
74,
72,
70,
69,
68,
67,
67,
66,
66,
66,
66
};
/*
* binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, but
* binomial(biggestSimpleBinomials[k] + 1, k) does.
*/
@VisibleForTesting
static final int[] biggestSimpleBinomials = {
Integer.MAX_VALUE,
Integer.MAX_VALUE,
Integer.MAX_VALUE,
2642246,
86251,
11724,
3218,
1313,
684,
419,
287,
214,
169,
139,
119,
105,
95,
87,
81,
76,
73,
70,
68,
66,
64,
63,
62,
62,
61,
61,
61
};
// These values were generated by using checkedMultiply to see when the simple multiply/divide
// algorithm would lead to an overflow.
static boolean fitsInInt(long x) {
return (int) x == x;
}
/**
* Returns the arithmetic mean of {@code x} and {@code y}, rounded toward negative infinity. This
* method is resilient to overflow.
*
* @since 14.0
*/
public static long mean(long x, long y) {
// Efficient method for computing the arithmetic mean.
// The alternative (x + y) / 2 fails for large values.
// The alternative (x + y) >>> 1 fails for negative values.
return (x & y) + ((x ^ y) >> 1);
}
/*
* This bitmask is used as an optimization for cheaply testing for divisiblity by 2, 3, or 5.
* Each bit is set to 1 for all remainders that indicate divisibility by 2, 3, or 5, so
* 1, 7, 11, 13, 17, 19, 23, 29 are set to 0. 30 and up don't matter because they won't be hit.
*/
private static final int SIEVE_30 =
~((1 << 1) | (1 << 7) | (1 << 11) | (1 << 13) | (1 << 17) | (1 << 19) | (1 << 23)
| (1 << 29));
/**
* Returns {@code true} if {@code n} is a prime number: an integer greater
* than one that cannot be factored into a product of smaller positive integers.
* Returns {@code false} if {@code n} is zero, one, or a composite number (one which can be
* factored into smaller positive integers).
*
* To test larger numbers, use {@link BigInteger#isProbablePrime}.
*
* @throws IllegalArgumentException if {@code n} is negative
* @since 20.0
*/
@GwtIncompatible // TODO
@Beta
public static boolean isPrime(long n) {
if (n < 2) {
checkNonNegative("n", n);
return false;
}
if (n == 2 || n == 3 || n == 5 || n == 7 || n == 11 || n == 13) {
return true;
}
if ((SIEVE_30 & (1 << (n % 30))) != 0) {
return false;
}
if (n % 7 == 0 || n % 11 == 0 || n % 13 == 0) {
return false;
}
if (n < 17 * 17) {
return true;
}
for (long[] baseSet : millerRabinBaseSets) {
if (n <= baseSet[0]) {
for (int i = 1; i < baseSet.length; i++) {
if (!MillerRabinTester.test(baseSet[i], n)) {
return false;
}
}
return true;
}
}
throw new AssertionError();
}
/*
* If n <= millerRabinBases[i][0], then testing n against bases millerRabinBases[i][1..] suffices
* to prove its primality. Values from miller-rabin.appspot.com.
*
* NOTE: We could get slightly better bases that would be treated as unsigned, but benchmarks
* showed negligible performance improvements.
*/
private static final long[][] millerRabinBaseSets = {
{291830, 126401071349994536L},
{885594168, 725270293939359937L, 3569819667048198375L},
{273919523040L, 15, 7363882082L, 992620450144556L},
{47636622961200L, 2, 2570940, 211991001, 3749873356L},
{
7999252175582850L,
2,
4130806001517L,
149795463772692060L,
186635894390467037L,
3967304179347715805L
},
{
585226005592931976L,
2,
123635709730000L,
9233062284813009L,
43835965440333360L,
761179012939631437L,
1263739024124850375L
},
{Long.MAX_VALUE, 2, 325, 9375, 28178, 450775, 9780504, 1795265022}
};
private enum MillerRabinTester {
/** Works for inputs ≤ FLOOR_SQRT_MAX_LONG. */
SMALL {
@Override
long mulMod(long a, long b, long m) {
/*
* NOTE(lowasser, 2015-Feb-12): Benchmarks suggest that changing this to
* UnsignedLongs.remainder and increasing the threshold to 2^32 doesn't pay for itself, and
* adding another enum constant hurts performance further -- I suspect because bimorphic
* implementation is a sweet spot for the JVM.
*/
return (a * b) % m;
}
@Override
long squareMod(long a, long m) {
return (a * a) % m;
}
},
/** Works for all nonnegative signed longs. */
LARGE {
/** Returns (a + b) mod m. Precondition: {@code 0 <= a}, {@code b < m < 2^63}. */
private long plusMod(long a, long b, long m) {
return (a >= m - b) ? (a + b - m) : (a + b);
}
/** Returns (a * 2^32) mod m. a may be any unsigned long. */
private long times2ToThe32Mod(long a, long m) {
int remainingPowersOf2 = 32;
do {
int shift = Math.min(remainingPowersOf2, Long.numberOfLeadingZeros(a));
// shift is either the number of powers of 2 left to multiply a by, or the biggest shift
// possible while keeping a in an unsigned long.
a = UnsignedLongs.remainder(a << shift, m);
remainingPowersOf2 -= shift;
} while (remainingPowersOf2 > 0);
return a;
}
@Override
long mulMod(long a, long b, long m) {
long aHi = a >>> 32; // < 2^31
long bHi = b >>> 32; // < 2^31
long aLo = a & 0xFFFFFFFFL; // < 2^32
long bLo = b & 0xFFFFFFFFL; // < 2^32
/*
* a * b == aHi * bHi * 2^64 + (aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo.
* == (aHi * bHi * 2^32 + aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo
*
* We carry out this computation in modular arithmetic. Since times2ToThe32Mod accepts any
* unsigned long, we don't have to do a mod on every operation, only when intermediate
* results can exceed 2^63.
*/
long result = times2ToThe32Mod(aHi * bHi /* < 2^62 */, m); // < m < 2^63
result += aHi * bLo; // aHi * bLo < 2^63, result < 2^64
if (result < 0) {
result = UnsignedLongs.remainder(result, m);
}
// result < 2^63 again
result += aLo * bHi; // aLo * bHi < 2^63, result < 2^64
result = times2ToThe32Mod(result, m); // result < m < 2^63
return plusMod(result, UnsignedLongs.remainder(aLo * bLo /* < 2^64 */, m), m);
}
@Override
long squareMod(long a, long m) {
long aHi = a >>> 32; // < 2^31
long aLo = a & 0xFFFFFFFFL; // < 2^32
/*
* a^2 == aHi^2 * 2^64 + aHi * aLo * 2^33 + aLo^2
* == (aHi^2 * 2^32 + aHi * aLo * 2) * 2^32 + aLo^2
* We carry out this computation in modular arithmetic. Since times2ToThe32Mod accepts any
* unsigned long, we don't have to do a mod on every operation, only when intermediate
* results can exceed 2^63.
*/
long result = times2ToThe32Mod(aHi * aHi /* < 2^62 */, m); // < m < 2^63
long hiLo = aHi * aLo * 2;
if (hiLo < 0) {
hiLo = UnsignedLongs.remainder(hiLo, m);
}
// hiLo < 2^63
result += hiLo; // result < 2^64
result = times2ToThe32Mod(result, m); // result < m < 2^63
return plusMod(result, UnsignedLongs.remainder(aLo * aLo /* < 2^64 */, m), m);
}
};
static boolean test(long base, long n) {
// Since base will be considered % n, it's okay if base > FLOOR_SQRT_MAX_LONG,
// so long as n <= FLOOR_SQRT_MAX_LONG.
return ((n <= FLOOR_SQRT_MAX_LONG) ? SMALL : LARGE).testWitness(base, n);
}
/** Returns a * b mod m. */
abstract long mulMod(long a, long b, long m);
/** Returns a^2 mod m. */
abstract long squareMod(long a, long m);
/** Returns a^p mod m. */
private long powMod(long a, long p, long m) {
long res = 1;
for (; p != 0; p >>= 1) {
if ((p & 1) != 0) {
res = mulMod(res, a, m);
}
a = squareMod(a, m);
}
return res;
}
/** Returns true if n is a strong probable prime relative to the specified base. */
private boolean testWitness(long base, long n) {
int r = Long.numberOfTrailingZeros(n - 1);
long d = (n - 1) >> r;
base %= n;
if (base == 0) {
return true;
}
// Calculate a := base^d mod n.
long a = powMod(base, d, n);
// n passes this test if
// base^d = 1 (mod n)
// or base^(2^j * d) = -1 (mod n) for some 0 <= j < r.
if (a == 1) {
return true;
}
int j = 0;
while (a != n - 1) {
if (++j == r) {
return false;
}
a = squareMod(a, n);
}
return true;
}
}
private LongMath() {}
}