1 /* 2 * Copyright (C) 2011 The Guava Authors 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 17 package com.google.common.math; 18 19 import static com.google.common.base.Preconditions.checkArgument; 20 import static com.google.common.math.MathPreconditions.checkNonNegative; 21 import static com.google.common.math.MathPreconditions.checkPositive; 22 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary; 23 import static java.lang.Math.min; 24 import static java.math.RoundingMode.HALF_EVEN; 25 import static java.math.RoundingMode.HALF_UP; 26 27 import com.google.common.annotations.GwtCompatible; 28 import com.google.common.annotations.VisibleForTesting; 29 30 import java.math.RoundingMode; 31 32 /** 33 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and 34 * named analogously to their {@code BigInteger} counterparts. 35 * 36 * <p>The implementations of many methods in this class are based on material from Henry S. Warren, 37 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002). 38 * 39 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in 40 * {@link IntMath} and {@link BigIntegerMath} respectively. For other common operations on 41 * {@code long} values, see {@link com.google.common.primitives.Longs}. 42 * 43 * @author Louis Wasserman 44 * @since 11.0 45 */ 46 @GwtCompatible(emulated = true) 47 public final class LongMath { 48 // NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, || 49 50 /** 51 * Returns {@code true} if {@code x} represents a power of two. 52 * 53 * <p>This differs from {@code Long.bitCount(x) == 1}, because 54 * {@code Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two. 55 */ isPowerOfTwo(long x)56 public static boolean isPowerOfTwo(long x) { 57 return x > 0 & (x & (x - 1)) == 0; 58 } 59 60 /** 61 * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a 62 * signed long. The implementation is branch-free, and benchmarks suggest it is measurably 63 * faster than the straightforward ternary expression. 64 */ 65 @VisibleForTesting lessThanBranchFree(long x, long y)66 static int lessThanBranchFree(long x, long y) { 67 // Returns the sign bit of x - y. 68 return (int) (~~(x - y) >>> (Long.SIZE - 1)); 69 } 70 71 /** 72 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode. 73 * 74 * @throws IllegalArgumentException if {@code x <= 0} 75 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x} 76 * is not a power of two 77 */ 78 @SuppressWarnings("fallthrough") 79 // TODO(kevinb): remove after this warning is disabled globally log2(long x, RoundingMode mode)80 public static int log2(long x, RoundingMode mode) { 81 checkPositive("x", x); 82 switch (mode) { 83 case UNNECESSARY: 84 checkRoundingUnnecessary(isPowerOfTwo(x)); 85 // fall through 86 case DOWN: 87 case FLOOR: 88 return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x); 89 90 case UP: 91 case CEILING: 92 return Long.SIZE - Long.numberOfLeadingZeros(x - 1); 93 94 case HALF_DOWN: 95 case HALF_UP: 96 case HALF_EVEN: 97 // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5 98 int leadingZeros = Long.numberOfLeadingZeros(x); 99 long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros; 100 // floor(2^(logFloor + 0.5)) 101 int logFloor = (Long.SIZE - 1) - leadingZeros; 102 return logFloor + lessThanBranchFree(cmp, x); 103 104 default: 105 throw new AssertionError("impossible"); 106 } 107 } 108 109 /** The biggest half power of two that fits into an unsigned long */ 110 @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L; 111 112 // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i))) 113 @VisibleForTesting static final byte[] maxLog10ForLeadingZeros = { 114 19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 115 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, 116 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, 0 }; 117 118 // halfPowersOf10[i] = largest long less than 10^(i + 0.5) 119 120 /** 121 * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if 122 * {@code a == 0 && b == 0}. 123 * 124 * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0} 125 */ gcd(long a, long b)126 public static long gcd(long a, long b) { 127 /* 128 * The reason we require both arguments to be >= 0 is because otherwise, what do you return on 129 * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't 130 * an int. 131 */ 132 checkNonNegative("a", a); 133 checkNonNegative("b", b); 134 if (a == 0) { 135 // 0 % b == 0, so b divides a, but the converse doesn't hold. 136 // BigInteger.gcd is consistent with this decision. 137 return b; 138 } else if (b == 0) { 139 return a; // similar logic 140 } 141 /* 142 * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. 143 * This is >60% faster than the Euclidean algorithm in benchmarks. 144 */ 145 int aTwos = Long.numberOfTrailingZeros(a); 146 a >>= aTwos; // divide out all 2s 147 int bTwos = Long.numberOfTrailingZeros(b); 148 b >>= bTwos; // divide out all 2s 149 while (a != b) { // both a, b are odd 150 // The key to the binary GCD algorithm is as follows: 151 // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b). 152 // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two. 153 154 // We bend over backwards to avoid branching, adapting a technique from 155 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax 156 157 long delta = a - b; // can't overflow, since a and b are nonnegative 158 159 long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1)); 160 // equivalent to Math.min(delta, 0) 161 162 a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b) 163 // a is now nonnegative and even 164 165 b += minDeltaOrZero; // sets b to min(old a, b) 166 a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b 167 } 168 return a << min(aTwos, bTwos); 169 } 170 171 @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L; 172 173 static final long[] factorials = { 174 1L, 175 1L, 176 1L * 2, 177 1L * 2 * 3, 178 1L * 2 * 3 * 4, 179 1L * 2 * 3 * 4 * 5, 180 1L * 2 * 3 * 4 * 5 * 6, 181 1L * 2 * 3 * 4 * 5 * 6 * 7, 182 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8, 183 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9, 184 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10, 185 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11, 186 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12, 187 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13, 188 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14, 189 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15, 190 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16, 191 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17, 192 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18, 193 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19, 194 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20 195 }; 196 197 /** 198 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and 199 * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}. 200 * 201 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n} 202 */ binomial(int n, int k)203 public static long binomial(int n, int k) { 204 checkNonNegative("n", n); 205 checkNonNegative("k", k); 206 checkArgument(k <= n, "k (%s) > n (%s)", k, n); 207 if (k > (n >> 1)) { 208 k = n - k; 209 } 210 switch (k) { 211 case 0: 212 return 1; 213 case 1: 214 return n; 215 default: 216 if (n < factorials.length) { 217 return factorials[n] / (factorials[k] * factorials[n - k]); 218 } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) { 219 return Long.MAX_VALUE; 220 } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) { 221 // guaranteed not to overflow 222 long result = n--; 223 for (int i = 2; i <= k; n--, i++) { 224 result *= n; 225 result /= i; 226 } 227 return result; 228 } else { 229 int nBits = LongMath.log2(n, RoundingMode.CEILING); 230 231 long result = 1; 232 long numerator = n--; 233 long denominator = 1; 234 235 int numeratorBits = nBits; 236 // This is an upper bound on log2(numerator, ceiling). 237 238 /* 239 * We want to do this in long math for speed, but want to avoid overflow. We adapt the 240 * technique previously used by BigIntegerMath: maintain separate numerator and 241 * denominator accumulators, multiplying the fraction into result when near overflow. 242 */ 243 for (int i = 2; i <= k; i++, n--) { 244 if (numeratorBits + nBits < Long.SIZE - 1) { 245 // It's definitely safe to multiply into numerator and denominator. 246 numerator *= n; 247 denominator *= i; 248 numeratorBits += nBits; 249 } else { 250 // It might not be safe to multiply into numerator and denominator, 251 // so multiply (numerator / denominator) into result. 252 result = multiplyFraction(result, numerator, denominator); 253 numerator = n; 254 denominator = i; 255 numeratorBits = nBits; 256 } 257 } 258 return multiplyFraction(result, numerator, denominator); 259 } 260 } 261 } 262 263 /** 264 * Returns (x * numerator / denominator), which is assumed to come out to an integral value. 265 */ multiplyFraction(long x, long numerator, long denominator)266 static long multiplyFraction(long x, long numerator, long denominator) { 267 if (x == 1) { 268 return numerator / denominator; 269 } 270 long commonDivisor = gcd(x, denominator); 271 x /= commonDivisor; 272 denominator /= commonDivisor; 273 // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact, 274 // so denominator must be a divisor of numerator. 275 return x * (numerator / denominator); 276 } 277 278 /* 279 * binomial(biggestBinomials[k], k) fits in a long, but not 280 * binomial(biggestBinomials[k] + 1, k). 281 */ 282 static final int[] biggestBinomials = 283 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 3810779, 121977, 16175, 4337, 1733, 284 887, 534, 361, 265, 206, 169, 143, 125, 111, 101, 94, 88, 83, 79, 76, 74, 72, 70, 69, 68, 285 67, 67, 66, 66, 66, 66}; 286 287 /* 288 * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, 289 * but binomial(biggestSimpleBinomials[k] + 1, k) does. 290 */ 291 @VisibleForTesting static final int[] biggestSimpleBinomials = 292 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2642246, 86251, 11724, 3218, 1313, 293 684, 419, 287, 214, 169, 139, 119, 105, 95, 87, 81, 76, 73, 70, 68, 66, 64, 63, 62, 62, 294 61, 61, 61}; 295 // These values were generated by using checkedMultiply to see when the simple multiply/divide 296 // algorithm would lead to an overflow. 297 fitsInInt(long x)298 static boolean fitsInInt(long x) { 299 return (int) x == x; 300 } 301 302 /** 303 * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward 304 * negative infinity. This method is resilient to overflow. 305 * 306 * @since 14.0 307 */ mean(long x, long y)308 public static long mean(long x, long y) { 309 // Efficient method for computing the arithmetic mean. 310 // The alternative (x + y) / 2 fails for large values. 311 // The alternative (x + y) >>> 1 fails for negative values. 312 return (x & y) + ((x ^ y) >> 1); 313 } 314 LongMath()315 private LongMath() {} 316 } 317