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
017package com.google.common.math;
018
019import static com.google.common.base.Preconditions.checkArgument;
020import static com.google.common.base.Preconditions.checkNotNull;
021import static com.google.common.math.MathPreconditions.checkNoOverflow;
022import static com.google.common.math.MathPreconditions.checkNonNegative;
023import static com.google.common.math.MathPreconditions.checkPositive;
024import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
025import static java.lang.Math.abs;
026import static java.lang.Math.min;
027import static java.math.RoundingMode.HALF_EVEN;
028import static java.math.RoundingMode.HALF_UP;
029
030import com.google.common.annotations.Beta;
031import com.google.common.annotations.GwtCompatible;
032import com.google.common.annotations.GwtIncompatible;
033import com.google.common.annotations.VisibleForTesting;
034
035import java.math.BigInteger;
036import 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)
054public 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}