001 /*
002 * Copyright (C) 2011 The Guava Authors
003 *
004 * Licensed under the Apache License, Version 2.0 (the "License");
005 * you may not use this file except in compliance with the License.
006 * You may obtain a copy of the License at
007 *
008 * http://www.apache.org/licenses/LICENSE-2.0
009 *
010 * Unless required by applicable law or agreed to in writing, software
011 * distributed under the License is distributed on an "AS IS" BASIS,
012 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013 * See the License for the specific language governing permissions and
014 * limitations under the License.
015 */
016
017 package com.google.common.math;
018
019 import static com.google.common.base.Preconditions.checkArgument;
020 import static com.google.common.base.Preconditions.checkNotNull;
021 import static com.google.common.math.MathPreconditions.checkNoOverflow;
022 import static com.google.common.math.MathPreconditions.checkNonNegative;
023 import static com.google.common.math.MathPreconditions.checkPositive;
024 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
025 import static java.lang.Math.abs;
026 import static java.lang.Math.min;
027 import static java.math.RoundingMode.HALF_EVEN;
028 import static java.math.RoundingMode.HALF_UP;
029
030 import com.google.common.annotations.Beta;
031 import com.google.common.annotations.GwtCompatible;
032 import com.google.common.annotations.GwtIncompatible;
033 import com.google.common.annotations.VisibleForTesting;
034
035 import java.math.BigInteger;
036 import java.math.RoundingMode;
037
038 /**
039 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
040 * named analogously to their {@code BigInteger} counterparts.
041 *
042 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
043 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
044 *
045 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in
046 * {@link IntMath} and {@link BigIntegerMath} respectively. For other common operations on
047 * {@code long} values, see {@link com.google.common.primitives.Longs}.
048 *
049 * @author Louis Wasserman
050 * @since 11.0
051 */
052 @Beta
053 @GwtCompatible(emulated = true)
054 public final class LongMath {
055 // NOTE: Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
056
057 /**
058 * Returns {@code true} if {@code x} represents a power of two.
059 *
060 * <p>This differs from {@code Long.bitCount(x) == 1}, because
061 * {@code Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
062 */
063 public static boolean isPowerOfTwo(long x) {
064 return x > 0 & (x & (x - 1)) == 0;
065 }
066
067 /**
068 * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
069 *
070 * @throws IllegalArgumentException if {@code x <= 0}
071 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
072 * is not a power of two
073 */
074 @SuppressWarnings("fallthrough")
075 public static int log2(long x, RoundingMode mode) {
076 checkPositive("x", x);
077 switch (mode) {
078 case UNNECESSARY:
079 checkRoundingUnnecessary(isPowerOfTwo(x));
080 // fall through
081 case DOWN:
082 case FLOOR:
083 return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
084
085 case UP:
086 case CEILING:
087 return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
088
089 case HALF_DOWN:
090 case HALF_UP:
091 case HALF_EVEN:
092 // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
093 int leadingZeros = Long.numberOfLeadingZeros(x);
094 long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
095 // floor(2^(logFloor + 0.5))
096 int logFloor = (Long.SIZE - 1) - leadingZeros;
097 return (x <= cmp) ? logFloor : logFloor + 1;
098
099 default:
100 throw new AssertionError("impossible");
101 }
102 }
103
104 /** The biggest half power of two that fits into an unsigned long */
105 @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
106
107 /**
108 * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
109 *
110 * @throws IllegalArgumentException if {@code x <= 0}
111 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
112 * is not a power of ten
113 */
114 @GwtIncompatible("TODO")
115 @SuppressWarnings("fallthrough")
116 public static int log10(long x, RoundingMode mode) {
117 checkPositive("x", x);
118 if (fitsInInt(x)) {
119 return IntMath.log10((int) x, mode);
120 }
121 int logFloor = log10Floor(x);
122 long floorPow = POWERS_OF_10[logFloor];
123 switch (mode) {
124 case UNNECESSARY:
125 checkRoundingUnnecessary(x == floorPow);
126 // fall through
127 case FLOOR:
128 case DOWN:
129 return logFloor;
130 case CEILING:
131 case UP:
132 return (x == floorPow) ? logFloor : logFloor + 1;
133 case HALF_DOWN:
134 case HALF_UP:
135 case HALF_EVEN:
136 // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
137 return (x <= HALF_POWERS_OF_10[logFloor]) ? logFloor : logFloor + 1;
138 default:
139 throw new AssertionError();
140 }
141 }
142
143 @GwtIncompatible("TODO")
144 static int log10Floor(long x) {
145 /*
146 * Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation.
147 *
148 * The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))),
149 * we can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x))
150 * is 6, then 64 <= x < 128, so floor(log10(x)) is either 1 or 2.
151 */
152 int y = MAX_LOG10_FOR_LEADING_ZEROS[Long.numberOfLeadingZeros(x)];
153 // y is the higher of the two possible values of floor(log10(x))
154
155 long sgn = (x - POWERS_OF_10[y]) >>> (Long.SIZE - 1);
156 /*
157 * sgn is the sign bit of x - 10^y; it is 1 if x < 10^y, and 0 otherwise. If x < 10^y, then we
158 * want the lower of the two possible values, or y - 1, otherwise, we want y.
159 */
160 return y - (int) sgn;
161 }
162
163 // MAX_LOG10_FOR_LEADING_ZEROS[i] == floor(log10(2^(Long.SIZE - i)))
164 @VisibleForTesting static final byte[] MAX_LOG10_FOR_LEADING_ZEROS = {
165 19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12,
166 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,
167 3, 3, 3, 3, 2, 2, 2, 1, 1, 1, 0, 0, 0 };
168
169 @GwtIncompatible("TODO")
170 @VisibleForTesting
171 static final long[] POWERS_OF_10 = {
172 1L,
173 10L,
174 100L,
175 1000L,
176 10000L,
177 100000L,
178 1000000L,
179 10000000L,
180 100000000L,
181 1000000000L,
182 10000000000L,
183 100000000000L,
184 1000000000000L,
185 10000000000000L,
186 100000000000000L,
187 1000000000000000L,
188 10000000000000000L,
189 100000000000000000L,
190 1000000000000000000L
191 };
192
193 // HALF_POWERS_OF_10[i] = largest long less than 10^(i + 0.5)
194 @GwtIncompatible("TODO")
195 @VisibleForTesting
196 static final long[] HALF_POWERS_OF_10 = {
197 3L,
198 31L,
199 316L,
200 3162L,
201 31622L,
202 316227L,
203 3162277L,
204 31622776L,
205 316227766L,
206 3162277660L,
207 31622776601L,
208 316227766016L,
209 3162277660168L,
210 31622776601683L,
211 316227766016837L,
212 3162277660168379L,
213 31622776601683793L,
214 316227766016837933L,
215 3162277660168379331L
216 };
217
218 /**
219 * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
220 * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
221 * time.
222 *
223 * @throws IllegalArgumentException if {@code k < 0}
224 */
225 @GwtIncompatible("TODO")
226 public static long pow(long b, int k) {
227 checkNonNegative("exponent", k);
228 if (-2 <= b && b <= 2) {
229 switch ((int) b) {
230 case 0:
231 return (k == 0) ? 1 : 0;
232 case 1:
233 return 1;
234 case (-1):
235 return ((k & 1) == 0) ? 1 : -1;
236 case 2:
237 return (k < Long.SIZE) ? 1L << k : 0;
238 case (-2):
239 if (k < Long.SIZE) {
240 return ((k & 1) == 0) ? 1L << k : -(1L << k);
241 } else {
242 return 0;
243 }
244 }
245 }
246 for (long accum = 1;; k >>= 1) {
247 switch (k) {
248 case 0:
249 return accum;
250 case 1:
251 return accum * b;
252 default:
253 accum *= ((k & 1) == 0) ? 1 : b;
254 b *= b;
255 }
256 }
257 }
258
259 /**
260 * Returns the square root of {@code x}, rounded with the specified rounding mode.
261 *
262 * @throws IllegalArgumentException if {@code x < 0}
263 * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and
264 * {@code sqrt(x)} is not an integer
265 */
266 @GwtIncompatible("TODO")
267 @SuppressWarnings("fallthrough")
268 public static long sqrt(long x, RoundingMode mode) {
269 checkNonNegative("x", x);
270 if (fitsInInt(x)) {
271 return IntMath.sqrt((int) x, mode);
272 }
273 long sqrtFloor = sqrtFloor(x);
274 switch (mode) {
275 case UNNECESSARY:
276 checkRoundingUnnecessary(sqrtFloor * sqrtFloor == x); // fall through
277 case FLOOR:
278 case DOWN:
279 return sqrtFloor;
280 case CEILING:
281 case UP:
282 return (sqrtFloor * sqrtFloor == x) ? sqrtFloor : sqrtFloor + 1;
283 case HALF_DOWN:
284 case HALF_UP:
285 case HALF_EVEN:
286 long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
287 /*
288 * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both
289 * x and halfSquare are integers, this is equivalent to testing whether or not x <=
290 * halfSquare. (We have to deal with overflow, though.)
291 */
292 return (halfSquare >= x | halfSquare < 0) ? sqrtFloor : sqrtFloor + 1;
293 default:
294 throw new AssertionError();
295 }
296 }
297
298 @GwtIncompatible("TODO")
299 private static long sqrtFloor(long x) {
300 // Hackers's Delight, Figure 11-1
301 long sqrt0 = (long) Math.sqrt(x);
302 // Precision can be lost in the cast to double, so we use this as a starting estimate.
303 long sqrt1 = (sqrt0 + (x / sqrt0)) >> 1;
304 if (sqrt1 == sqrt0) {
305 return sqrt0;
306 }
307 do {
308 sqrt0 = sqrt1;
309 sqrt1 = (sqrt0 + (x / sqrt0)) >> 1;
310 } while (sqrt1 < sqrt0);
311 return sqrt0;
312 }
313
314 /**
315 * Returns the result of dividing {@code p} by {@code q}, rounding using the specified
316 * {@code RoundingMode}.
317 *
318 * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
319 * is not an integer multiple of {@code b}
320 */
321 @GwtIncompatible("TODO")
322 @SuppressWarnings("fallthrough")
323 public static long divide(long p, long q, RoundingMode mode) {
324 checkNotNull(mode);
325 long div = p / q; // throws if q == 0
326 long rem = p - q * div; // equals p % q
327
328 if (rem == 0) {
329 return div;
330 }
331
332 /*
333 * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
334 * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
335 * p / q.
336 *
337 * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
338 */
339 int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
340 boolean increment;
341 switch (mode) {
342 case UNNECESSARY:
343 checkRoundingUnnecessary(rem == 0);
344 // fall through
345 case DOWN:
346 increment = false;
347 break;
348 case UP:
349 increment = true;
350 break;
351 case CEILING:
352 increment = signum > 0;
353 break;
354 case FLOOR:
355 increment = signum < 0;
356 break;
357 case HALF_EVEN:
358 case HALF_DOWN:
359 case HALF_UP:
360 long absRem = abs(rem);
361 long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
362 // subtracting two nonnegative longs can't overflow
363 // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
364 if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
365 increment = (mode == HALF_UP | (mode == HALF_EVEN & (div & 1) != 0));
366 } else {
367 increment = cmpRemToHalfDivisor > 0; // closer to the UP value
368 }
369 break;
370 default:
371 throw new AssertionError();
372 }
373 return increment ? div + signum : div;
374 }
375
376 /**
377 * Returns {@code x mod m}. This differs from {@code x % m} in that it always returns a
378 * non-negative result.
379 *
380 * <p>For example:
381 *
382 * <pre> {@code
383 *
384 * mod(7, 4) == 3
385 * mod(-7, 4) == 1
386 * mod(-1, 4) == 3
387 * mod(-8, 4) == 0
388 * mod(8, 4) == 0}</pre>
389 *
390 * @throws ArithmeticException if {@code m <= 0}
391 */
392 @GwtIncompatible("TODO")
393 public static int mod(long x, int m) {
394 // Cast is safe because the result is guaranteed in the range [0, m)
395 return (int) mod(x, (long) m);
396 }
397
398 /**
399 * Returns {@code x mod m}. This differs from {@code x % m} in that it always returns a
400 * non-negative result.
401 *
402 * <p>For example:
403 *
404 * <pre> {@code
405 *
406 * mod(7, 4) == 3
407 * mod(-7, 4) == 1
408 * mod(-1, 4) == 3
409 * mod(-8, 4) == 0
410 * mod(8, 4) == 0}</pre>
411 *
412 * @throws ArithmeticException if {@code m <= 0}
413 */
414 @GwtIncompatible("TODO")
415 public static long mod(long x, long m) {
416 if (m <= 0) {
417 throw new ArithmeticException("Modulus " + m + " must be > 0");
418 }
419 long result = x % m;
420 return (result >= 0) ? result : result + m;
421 }
422
423 /**
424 * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if
425 * {@code a == 0 && b == 0}.
426 *
427 * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
428 */
429 @GwtIncompatible("TODO")
430 public static long gcd(long a, long b) {
431 /*
432 * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
433 * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't
434 * an int.
435 */
436 checkNonNegative("a", a);
437 checkNonNegative("b", b);
438 if (a == 0) {
439 // 0 % b == 0, so b divides a, but the converse doesn't hold.
440 // BigInteger.gcd is consistent with this decision.
441 return b;
442 } else if (b == 0) {
443 return a; // similar logic
444 }
445 /*
446 * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm.
447 * This is >60% faster than the Euclidean algorithm in benchmarks.
448 */
449 int aTwos = Long.numberOfTrailingZeros(a);
450 a >>= aTwos; // divide out all 2s
451 int bTwos = Long.numberOfTrailingZeros(b);
452 b >>= bTwos; // divide out all 2s
453 while (a != b) { // both a, b are odd
454 // The key to the binary GCD algorithm is as follows:
455 // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
456 // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
457
458 // We bend over backwards to avoid branching, adapting a technique from
459 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
460
461 long delta = a - b; // can't overflow, since a and b are nonnegative
462
463 long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
464 // equivalent to Math.min(delta, 0)
465
466 a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
467 // a is now nonnegative and even
468
469 b += minDeltaOrZero; // sets b to min(old a, b)
470 a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
471 }
472 return a << min(aTwos, bTwos);
473 }
474
475 /**
476 * Returns the sum of {@code a} and {@code b}, provided it does not overflow.
477 *
478 * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
479 */
480 @GwtIncompatible("TODO")
481 public static long checkedAdd(long a, long b) {
482 long result = a + b;
483 checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0);
484 return result;
485 }
486
487 /**
488 * Returns the difference of {@code a} and {@code b}, provided it does not overflow.
489 *
490 * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
491 */
492 @GwtIncompatible("TODO")
493 public static long checkedSubtract(long a, long b) {
494 long result = a - b;
495 checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0);
496 return result;
497 }
498
499 /**
500 * Returns the product of {@code a} and {@code b}, provided it does not overflow.
501 *
502 * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
503 */
504 @GwtIncompatible("TODO")
505 public static long checkedMultiply(long a, long b) {
506 // Hacker's Delight, Section 2-12
507 int leadingZeros = Long.numberOfLeadingZeros(a) + Long.numberOfLeadingZeros(~a)
508 + Long.numberOfLeadingZeros(b) + Long.numberOfLeadingZeros(~b);
509 /*
510 * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
511 * bad. We do the leadingZeros check to avoid the division below if at all possible.
512 *
513 * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
514 * care of all a < 0 with their own check, because in particular, the case a == -1 will
515 * incorrectly pass the division check below.
516 *
517 * In all other cases, we check that either a is 0 or the result is consistent with division.
518 */
519 if (leadingZeros > Long.SIZE + 1) {
520 return a * b;
521 }
522 checkNoOverflow(leadingZeros >= Long.SIZE);
523 checkNoOverflow(a >= 0 | b != Long.MIN_VALUE);
524 long result = a * b;
525 checkNoOverflow(a == 0 || result / a == b);
526 return result;
527 }
528
529 /**
530 * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
531 *
532 * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed
533 * {@code long} arithmetic
534 */
535 @GwtIncompatible("TODO")
536 public static long checkedPow(long b, int k) {
537 checkNonNegative("exponent", k);
538 if (b >= -2 & b <= 2) {
539 switch ((int) b) {
540 case 0:
541 return (k == 0) ? 1 : 0;
542 case 1:
543 return 1;
544 case (-1):
545 return ((k & 1) == 0) ? 1 : -1;
546 case 2:
547 checkNoOverflow(k < Long.SIZE - 1);
548 return 1L << k;
549 case (-2):
550 checkNoOverflow(k < Long.SIZE);
551 return ((k & 1) == 0) ? (1L << k) : (-1L << k);
552 }
553 }
554 long accum = 1;
555 while (true) {
556 switch (k) {
557 case 0:
558 return accum;
559 case 1:
560 return checkedMultiply(accum, b);
561 default:
562 if ((k & 1) != 0) {
563 accum = checkedMultiply(accum, b);
564 }
565 k >>= 1;
566 if (k > 0) {
567 checkNoOverflow(b <= FLOOR_SQRT_MAX_LONG);
568 b *= b;
569 }
570 }
571 }
572 }
573
574 @GwtIncompatible("TODO")
575 @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
576
577 /**
578 * Returns {@code n!}, that is, the product of the first {@code n} positive
579 * integers, {@code 1} if {@code n == 0}, or {@link Long#MAX_VALUE} if the
580 * result does not fit in a {@code long}.
581 *
582 * @throws IllegalArgumentException if {@code n < 0}
583 */
584 @GwtIncompatible("TODO")
585 public static long factorial(int n) {
586 checkNonNegative("n", n);
587 return (n < FACTORIALS.length) ? FACTORIALS[n] : Long.MAX_VALUE;
588 }
589
590 static final long[] FACTORIALS = {
591 1L,
592 1L,
593 1L * 2,
594 1L * 2 * 3,
595 1L * 2 * 3 * 4,
596 1L * 2 * 3 * 4 * 5,
597 1L * 2 * 3 * 4 * 5 * 6,
598 1L * 2 * 3 * 4 * 5 * 6 * 7,
599 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
600 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
601 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
602 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
603 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
604 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
605 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
606 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
607 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
608 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
609 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
610 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
611 1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
612 };
613
614 /**
615 * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
616 * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
617 *
618 * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
619 */
620 public static long binomial(int n, int k) {
621 checkNonNegative("n", n);
622 checkNonNegative("k", k);
623 checkArgument(k <= n, "k (%s) > n (%s)", k, n);
624 if (k > (n >> 1)) {
625 k = n - k;
626 }
627 if (k >= BIGGEST_BINOMIALS.length || n > BIGGEST_BINOMIALS[k]) {
628 return Long.MAX_VALUE;
629 }
630 long result = 1;
631 if (k < BIGGEST_SIMPLE_BINOMIALS.length && n <= BIGGEST_SIMPLE_BINOMIALS[k]) {
632 // guaranteed not to overflow
633 for (int i = 0; i < k; i++) {
634 result *= n - i;
635 result /= i + 1;
636 }
637 } else {
638 // We want to do this in long math for speed, but want to avoid overflow.
639 // Dividing by the GCD suffices to avoid overflow in all the remaining cases.
640 for (int i = 1; i <= k; i++, n--) {
641 int d = IntMath.gcd(n, i);
642 result /= i / d; // (i/d) is guaranteed to divide result
643 result *= n / d;
644 }
645 }
646 return result;
647 }
648
649 /*
650 * binomial(BIGGEST_BINOMIALS[k], k) fits in a long, but not
651 * binomial(BIGGEST_BINOMIALS[k] + 1, k).
652 */
653 static final int[] BIGGEST_BINOMIALS =
654 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 3810779, 121977, 16175, 4337, 1733,
655 887, 534, 361, 265, 206, 169, 143, 125, 111, 101, 94, 88, 83, 79, 76, 74, 72, 70, 69, 68,
656 67, 67, 66, 66, 66, 66};
657
658 /*
659 * binomial(BIGGEST_SIMPLE_BINOMIALS[k], k) doesn't need to use the slower GCD-based impl,
660 * but binomial(BIGGEST_SIMPLE_BINOMIALS[k] + 1, k) does.
661 */
662 @VisibleForTesting static final int[] BIGGEST_SIMPLE_BINOMIALS =
663 {Integer.MAX_VALUE, Integer.MAX_VALUE, Integer.MAX_VALUE, 2642246, 86251, 11724, 3218, 1313,
664 684, 419, 287, 214, 169, 139, 119, 105, 95, 87, 81, 76, 73, 70, 68, 66, 64, 63, 62, 62,
665 61, 61, 61};
666 // These values were generated by using checkedMultiply to see when the simple multiply/divide
667 // algorithm would lead to an overflow.
668
669 @GwtIncompatible("TODO")
670 static boolean fitsInInt(long x) {
671 return (int) x == x;
672 }
673
674 private LongMath() {}
675 }