001/*
002 * Copyright (C) 2011 The Guava Authors
003 *
004 * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
005 * in compliance with the License. You may obtain a copy of the License at
006 *
007 * http://www.apache.org/licenses/LICENSE-2.0
008 *
009 * Unless required by applicable law or agreed to in writing, software distributed under the License
010 * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
011 * or implied. See the License for the specific language governing permissions and limitations under
012 * the License.
013 */
014
015package com.google.common.math;
016
017import static com.google.common.base.Preconditions.checkArgument;
018import static com.google.common.base.Preconditions.checkNotNull;
019import static com.google.common.math.MathPreconditions.checkNoOverflow;
020import static com.google.common.math.MathPreconditions.checkNonNegative;
021import static com.google.common.math.MathPreconditions.checkPositive;
022import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
023import static java.lang.Math.abs;
024import static java.lang.Math.min;
025import static java.math.RoundingMode.HALF_EVEN;
026import static java.math.RoundingMode.HALF_UP;
027
028import com.google.common.annotations.GwtCompatible;
029import com.google.common.annotations.GwtIncompatible;
030import com.google.common.annotations.VisibleForTesting;
031import com.google.common.primitives.UnsignedLongs;
032import java.math.BigInteger;
033import java.math.RoundingMode;
034
035/**
036 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
037 * named analogously to their {@code BigInteger} counterparts.
038 *
039 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
040 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
041 *
042 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in {@link
043 * IntMath} and {@link BigIntegerMath} respectively. For other common operations on {@code long}
044 * values, see {@link com.google.common.primitives.Longs}.
045 *
046 * @author Louis Wasserman
047 * @since 11.0
048 */
049@GwtCompatible(emulated = true)
050@ElementTypesAreNonnullByDefault
051public final class LongMath {
052  @VisibleForTesting static final long MAX_SIGNED_POWER_OF_TWO = 1L << (Long.SIZE - 2);
053
054  /**
055   * Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to
056   * {@code checkedPow(2, log2(x, CEILING))}.
057   *
058   * @throws IllegalArgumentException if {@code x <= 0}
059   * @throws ArithmeticException of the next-higher power of two is not representable as a {@code
060   *     long}, i.e. when {@code x > 2^62}
061   * @since 20.0
062   */
063  public static long ceilingPowerOfTwo(long x) {
064    checkPositive("x", x);
065    if (x > MAX_SIGNED_POWER_OF_TWO) {
066      throw new ArithmeticException("ceilingPowerOfTwo(" + x + ") is not representable as a long");
067    }
068    return 1L << -Long.numberOfLeadingZeros(x - 1);
069  }
070
071  /**
072   * Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code
073   * checkedPow(2, log2(x, FLOOR))}.
074   *
075   * @throws IllegalArgumentException if {@code x <= 0}
076   * @since 20.0
077   */
078  public static long floorPowerOfTwo(long x) {
079    checkPositive("x", x);
080
081    // Long.highestOneBit was buggy on GWT.  We've fixed it, but I'm not certain when the fix will
082    // be released.
083    return 1L << ((Long.SIZE - 1) - Long.numberOfLeadingZeros(x));
084  }
085
086  /**
087   * Returns {@code true} if {@code x} represents a power of two.
088   *
089   * <p>This differs from {@code Long.bitCount(x) == 1}, because {@code
090   * Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
091   */
092  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
093  @SuppressWarnings("ShortCircuitBoolean")
094  public static boolean isPowerOfTwo(long x) {
095    return x > 0 & (x & (x - 1)) == 0;
096  }
097
098  /**
099   * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
100   * signed long. The implementation is branch-free, and benchmarks suggest it is measurably faster
101   * than the straightforward ternary expression.
102   */
103  @VisibleForTesting
104  static int lessThanBranchFree(long x, long y) {
105    // Returns the sign bit of x - y.
106    return (int) (~~(x - y) >>> (Long.SIZE - 1));
107  }
108
109  /**
110   * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
111   *
112   * @throws IllegalArgumentException if {@code x <= 0}
113   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
114   *     is not a power of two
115   */
116  @SuppressWarnings("fallthrough")
117  // TODO(kevinb): remove after this warning is disabled globally
118  public static int log2(long x, RoundingMode mode) {
119    checkPositive("x", x);
120    switch (mode) {
121      case UNNECESSARY:
122        checkRoundingUnnecessary(isPowerOfTwo(x));
123        // fall through
124      case DOWN:
125      case FLOOR:
126        return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
127
128      case UP:
129      case CEILING:
130        return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
131
132      case HALF_DOWN:
133      case HALF_UP:
134      case HALF_EVEN:
135        // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
136        int leadingZeros = Long.numberOfLeadingZeros(x);
137        long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
138        // floor(2^(logFloor + 0.5))
139        int logFloor = (Long.SIZE - 1) - leadingZeros;
140        return logFloor + lessThanBranchFree(cmp, x);
141    }
142    throw new AssertionError("impossible");
143  }
144
145  /** The biggest half power of two that fits into an unsigned long */
146  @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
147
148  /**
149   * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
150   *
151   * @throws IllegalArgumentException if {@code x <= 0}
152   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
153   *     is not a power of ten
154   */
155  @GwtIncompatible // TODO
156  @SuppressWarnings("fallthrough")
157  // TODO(kevinb): remove after this warning is disabled globally
158  public static int log10(long x, RoundingMode mode) {
159    checkPositive("x", x);
160    int logFloor = log10Floor(x);
161    long floorPow = powersOf10[logFloor];
162    switch (mode) {
163      case UNNECESSARY:
164        checkRoundingUnnecessary(x == floorPow);
165        // fall through
166      case FLOOR:
167      case DOWN:
168        return logFloor;
169      case CEILING:
170      case UP:
171        return logFloor + lessThanBranchFree(floorPow, x);
172      case HALF_DOWN:
173      case HALF_UP:
174      case HALF_EVEN:
175        // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
176        return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x);
177    }
178    throw new AssertionError();
179  }
180
181  @GwtIncompatible // TODO
182  static int log10Floor(long x) {
183    /*
184     * Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation.
185     *
186     * The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))), we
187     * can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x)) is 6,
188     * then 64 <= x < 128, so floor(log10(x)) is either 1 or 2.
189     */
190    int y = maxLog10ForLeadingZeros[Long.numberOfLeadingZeros(x)];
191    /*
192     * y is the higher of the two possible values of floor(log10(x)). If x < 10^y, then we want the
193     * lower of the two possible values, or y - 1, otherwise, we want y.
194     */
195    return y - lessThanBranchFree(x, powersOf10[y]);
196  }
197
198  // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i)))
199  @VisibleForTesting
200  static final byte[] maxLog10ForLeadingZeros = {
201    19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 12,
202    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,
203    3, 2, 2, 2, 1, 1, 1, 0, 0, 0
204  };
205
206  @GwtIncompatible // TODO
207  @VisibleForTesting
208  static final long[] powersOf10 = {
209    1L,
210    10L,
211    100L,
212    1000L,
213    10000L,
214    100000L,
215    1000000L,
216    10000000L,
217    100000000L,
218    1000000000L,
219    10000000000L,
220    100000000000L,
221    1000000000000L,
222    10000000000000L,
223    100000000000000L,
224    1000000000000000L,
225    10000000000000000L,
226    100000000000000000L,
227    1000000000000000000L
228  };
229
230  // halfPowersOf10[i] = largest long less than 10^(i + 0.5)
231  @GwtIncompatible // TODO
232  @VisibleForTesting
233  static final long[] halfPowersOf10 = {
234    3L,
235    31L,
236    316L,
237    3162L,
238    31622L,
239    316227L,
240    3162277L,
241    31622776L,
242    316227766L,
243    3162277660L,
244    31622776601L,
245    316227766016L,
246    3162277660168L,
247    31622776601683L,
248    316227766016837L,
249    3162277660168379L,
250    31622776601683793L,
251    316227766016837933L,
252    3162277660168379331L
253  };
254
255  /**
256   * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
257   * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
258   * time.
259   *
260   * @throws IllegalArgumentException if {@code k < 0}
261   */
262  @GwtIncompatible // TODO
263  public static long pow(long b, int k) {
264    checkNonNegative("exponent", k);
265    if (-2 <= b && b <= 2) {
266      switch ((int) b) {
267        case 0:
268          return (k == 0) ? 1 : 0;
269        case 1:
270          return 1;
271        case (-1):
272          return ((k & 1) == 0) ? 1 : -1;
273        case 2:
274          return (k < Long.SIZE) ? 1L << k : 0;
275        case (-2):
276          if (k < Long.SIZE) {
277            return ((k & 1) == 0) ? 1L << k : -(1L << k);
278          } else {
279            return 0;
280          }
281        default:
282          throw new AssertionError();
283      }
284    }
285    for (long accum = 1; ; k >>= 1) {
286      switch (k) {
287        case 0:
288          return accum;
289        case 1:
290          return accum * b;
291        default:
292          accum *= ((k & 1) == 0) ? 1 : b;
293          b *= b;
294      }
295    }
296  }
297
298  /**
299   * Returns the square root of {@code x}, rounded with the specified rounding mode.
300   *
301   * @throws IllegalArgumentException if {@code x < 0}
302   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
303   *     sqrt(x)} is not an integer
304   */
305  @GwtIncompatible // TODO
306  public static long sqrt(long x, RoundingMode mode) {
307    checkNonNegative("x", x);
308    if (fitsInInt(x)) {
309      return IntMath.sqrt((int) x, mode);
310    }
311    /*
312     * Let k be the true value of floor(sqrt(x)), so that
313     *
314     *            k * k <= x          <  (k + 1) * (k + 1)
315     * (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1))
316     *          since casting to double is nondecreasing.
317     *          Note that the right-hand inequality is no longer strict.
318     * Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1))
319     *          since Math.sqrt is monotonic.
320     * (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1))
321     *          since casting to long is monotonic
322     * k <= (long) Math.sqrt(x) <= k + 1
323     *          since (long) Math.sqrt(k * k) == k, as checked exhaustively in
324     *          {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect}
325     */
326    long guess = (long) Math.sqrt((double) x);
327    // Note: guess is always <= FLOOR_SQRT_MAX_LONG.
328    long guessSquared = guess * guess;
329    // Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is
330    // faster here than using lessThanBranchFree.
331    switch (mode) {
332      case UNNECESSARY:
333        checkRoundingUnnecessary(guessSquared == x);
334        return guess;
335      case FLOOR:
336      case DOWN:
337        if (x < guessSquared) {
338          return guess - 1;
339        }
340        return guess;
341      case CEILING:
342      case UP:
343        if (x > guessSquared) {
344          return guess + 1;
345        }
346        return guess;
347      case HALF_DOWN:
348      case HALF_UP:
349      case HALF_EVEN:
350        long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0);
351        long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
352        /*
353         * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
354         * and halfSquare are integers, this is equivalent to testing whether or not x <=
355         * halfSquare. (We have to deal with overflow, though.)
356         *
357         * If we treat halfSquare as an unsigned long, we know that
358         *            sqrtFloor^2 <= x < (sqrtFloor + 1)^2
359         * halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1
360         * so |x - halfSquare| <= sqrtFloor.  Therefore, it's safe to treat x - halfSquare as a
361         * signed long, so lessThanBranchFree is safe for use.
362         */
363        return sqrtFloor + lessThanBranchFree(halfSquare, x);
364    }
365    throw new AssertionError();
366  }
367
368  /**
369   * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
370   * RoundingMode}.
371   *
372   * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
373   *     is not an integer multiple of {@code b}
374   */
375  @GwtIncompatible // TODO
376  @SuppressWarnings("fallthrough")
377  public static long divide(long p, long q, RoundingMode mode) {
378    checkNotNull(mode);
379    long div = p / q; // throws if q == 0
380    long rem = p - q * div; // equals p % q
381
382    if (rem == 0) {
383      return div;
384    }
385
386    /*
387     * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
388     * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
389     * p / q.
390     *
391     * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
392     */
393    int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
394    boolean increment;
395    switch (mode) {
396      case UNNECESSARY:
397        checkRoundingUnnecessary(rem == 0);
398        // fall through
399      case DOWN:
400        increment = false;
401        break;
402      case UP:
403        increment = true;
404        break;
405      case CEILING:
406        increment = signum > 0;
407        break;
408      case FLOOR:
409        increment = signum < 0;
410        break;
411      case HALF_EVEN:
412      case HALF_DOWN:
413      case HALF_UP:
414        long absRem = abs(rem);
415        long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
416        // subtracting two nonnegative longs can't overflow
417        // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
418        if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
419          increment = (mode == HALF_UP || (mode == HALF_EVEN && (div & 1) != 0));
420        } else {
421          increment = cmpRemToHalfDivisor > 0; // closer to the UP value
422        }
423        break;
424      default:
425        throw new AssertionError();
426    }
427    return increment ? div + signum : div;
428  }
429
430  /**
431   * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
432   * m}, which might be negative.
433   *
434   * <p>For example:
435   *
436   * <pre>{@code
437   * mod(7, 4) == 3
438   * mod(-7, 4) == 1
439   * mod(-1, 4) == 3
440   * mod(-8, 4) == 0
441   * mod(8, 4) == 0
442   * }</pre>
443   *
444   * @throws ArithmeticException if {@code m <= 0}
445   * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
446   *     Remainder Operator</a>
447   */
448  @GwtIncompatible // TODO
449  public static int mod(long x, int m) {
450    // Cast is safe because the result is guaranteed in the range [0, m)
451    return (int) mod(x, (long) m);
452  }
453
454  /**
455   * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
456   * m}, which might be negative.
457   *
458   * <p>For example:
459   *
460   * <pre>{@code
461   * mod(7, 4) == 3
462   * mod(-7, 4) == 1
463   * mod(-1, 4) == 3
464   * mod(-8, 4) == 0
465   * mod(8, 4) == 0
466   * }</pre>
467   *
468   * @throws ArithmeticException if {@code m <= 0}
469   * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
470   *     Remainder Operator</a>
471   */
472  @GwtIncompatible // TODO
473  public static long mod(long x, long m) {
474    if (m <= 0) {
475      throw new ArithmeticException("Modulus must be positive");
476    }
477    long result = x % m;
478    return (result >= 0) ? result : result + m;
479  }
480
481  /**
482   * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if {@code a == 0 && b ==
483   * 0}.
484   *
485   * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
486   */
487  public static long gcd(long a, long b) {
488    /*
489     * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
490     * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't an
491     * int.
492     */
493    checkNonNegative("a", a);
494    checkNonNegative("b", b);
495    if (a == 0) {
496      // 0 % b == 0, so b divides a, but the converse doesn't hold.
497      // BigInteger.gcd is consistent with this decision.
498      return b;
499    } else if (b == 0) {
500      return a; // similar logic
501    }
502    /*
503     * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. This is
504     * >60% faster than the Euclidean algorithm in benchmarks.
505     */
506    int aTwos = Long.numberOfTrailingZeros(a);
507    a >>= aTwos; // divide out all 2s
508    int bTwos = Long.numberOfTrailingZeros(b);
509    b >>= bTwos; // divide out all 2s
510    while (a != b) { // both a, b are odd
511      // The key to the binary GCD algorithm is as follows:
512      // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
513      // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
514
515      // We bend over backwards to avoid branching, adapting a technique from
516      // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
517
518      long delta = a - b; // can't overflow, since a and b are nonnegative
519
520      long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
521      // equivalent to Math.min(delta, 0)
522
523      a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
524      // a is now nonnegative and even
525
526      b += minDeltaOrZero; // sets b to min(old a, b)
527      a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
528    }
529    return a << min(aTwos, bTwos);
530  }
531
532  /**
533   * Returns the sum of {@code a} and {@code b}, provided it does not overflow.
534   *
535   * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
536   */
537  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
538  @SuppressWarnings("ShortCircuitBoolean")
539  public static long checkedAdd(long a, long b) {
540    long result = a + b;
541    checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0, "checkedAdd", a, b);
542    return result;
543  }
544
545  /**
546   * Returns the difference of {@code a} and {@code b}, provided it does not overflow.
547   *
548   * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
549   */
550  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
551  @SuppressWarnings("ShortCircuitBoolean")
552  public static long checkedSubtract(long a, long b) {
553    long result = a - b;
554    checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0, "checkedSubtract", a, b);
555    return result;
556  }
557
558  /**
559   * Returns the product of {@code a} and {@code b}, provided it does not overflow.
560   *
561   * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
562   */
563  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
564  @SuppressWarnings("ShortCircuitBoolean")
565  public static long checkedMultiply(long a, long b) {
566    // Hacker's Delight, Section 2-12
567    int leadingZeros =
568        Long.numberOfLeadingZeros(a)
569            + Long.numberOfLeadingZeros(~a)
570            + Long.numberOfLeadingZeros(b)
571            + Long.numberOfLeadingZeros(~b);
572    /*
573     * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
574     * bad. We do the leadingZeros check to avoid the division below if at all possible.
575     *
576     * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
577     * care of all a < 0 with their own check, because in particular, the case a == -1 will
578     * incorrectly pass the division check below.
579     *
580     * In all other cases, we check that either a is 0 or the result is consistent with division.
581     */
582    if (leadingZeros > Long.SIZE + 1) {
583      return a * b;
584    }
585    checkNoOverflow(leadingZeros >= Long.SIZE, "checkedMultiply", a, b);
586    checkNoOverflow(a >= 0 | b != Long.MIN_VALUE, "checkedMultiply", a, b);
587    long result = a * b;
588    checkNoOverflow(a == 0 || result / a == b, "checkedMultiply", a, b);
589    return result;
590  }
591
592  /**
593   * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
594   *
595   * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed {@code
596   *     long} arithmetic
597   */
598  @GwtIncompatible // TODO
599  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
600  @SuppressWarnings("ShortCircuitBoolean")
601  public static long checkedPow(long b, int k) {
602    checkNonNegative("exponent", k);
603    if (b >= -2 & b <= 2) {
604      switch ((int) b) {
605        case 0:
606          return (k == 0) ? 1 : 0;
607        case 1:
608          return 1;
609        case (-1):
610          return ((k & 1) == 0) ? 1 : -1;
611        case 2:
612          checkNoOverflow(k < Long.SIZE - 1, "checkedPow", b, k);
613          return 1L << k;
614        case (-2):
615          checkNoOverflow(k < Long.SIZE, "checkedPow", b, k);
616          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
617        default:
618          throw new AssertionError();
619      }
620    }
621    long accum = 1;
622    while (true) {
623      switch (k) {
624        case 0:
625          return accum;
626        case 1:
627          return checkedMultiply(accum, b);
628        default:
629          if ((k & 1) != 0) {
630            accum = checkedMultiply(accum, b);
631          }
632          k >>= 1;
633          if (k > 0) {
634            checkNoOverflow(
635                -FLOOR_SQRT_MAX_LONG <= b && b <= FLOOR_SQRT_MAX_LONG, "checkedPow", b, k);
636            b *= b;
637          }
638      }
639    }
640  }
641
642  /**
643   * Returns the sum of {@code a} and {@code b} unless it would overflow or underflow in which case
644   * {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
645   *
646   * @since 20.0
647   */
648  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
649  @SuppressWarnings("ShortCircuitBoolean")
650  public static long saturatedAdd(long a, long b) {
651    long naiveSum = a + b;
652    if ((a ^ b) < 0 | (a ^ naiveSum) >= 0) {
653      // If a and b have different signs or a has the same sign as the result then there was no
654      // overflow, return.
655      return naiveSum;
656    }
657    // we did over/under flow, if the sign is negative we should return MAX otherwise MIN
658    return Long.MAX_VALUE + ((naiveSum >>> (Long.SIZE - 1)) ^ 1);
659  }
660
661  /**
662   * Returns the difference of {@code a} and {@code b} unless it would overflow or underflow in
663   * which case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
664   *
665   * @since 20.0
666   */
667  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
668  @SuppressWarnings("ShortCircuitBoolean")
669  public static long saturatedSubtract(long a, long b) {
670    long naiveDifference = a - b;
671    if ((a ^ b) >= 0 | (a ^ naiveDifference) >= 0) {
672      // If a and b have the same signs or a has the same sign as the result then there was no
673      // overflow, return.
674      return naiveDifference;
675    }
676    // we did over/under flow
677    return Long.MAX_VALUE + ((naiveDifference >>> (Long.SIZE - 1)) ^ 1);
678  }
679
680  /**
681   * Returns the product of {@code a} and {@code b} unless it would overflow or underflow in which
682   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
683   *
684   * @since 20.0
685   */
686  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
687  @SuppressWarnings("ShortCircuitBoolean")
688  public static long saturatedMultiply(long a, long b) {
689    // see checkedMultiply for explanation
690    int leadingZeros =
691        Long.numberOfLeadingZeros(a)
692            + Long.numberOfLeadingZeros(~a)
693            + Long.numberOfLeadingZeros(b)
694            + Long.numberOfLeadingZeros(~b);
695    if (leadingZeros > Long.SIZE + 1) {
696      return a * b;
697    }
698    // the return value if we will overflow (which we calculate by overflowing a long :) )
699    long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
700    if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
701      // overflow
702      return limit;
703    }
704    long result = a * b;
705    if (a == 0 || result / a == b) {
706      return result;
707    }
708    return limit;
709  }
710
711  /**
712   * Returns the {@code b} to the {@code k}th power, unless it would overflow or underflow in which
713   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
714   *
715   * @since 20.0
716   */
717  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
718  @SuppressWarnings("ShortCircuitBoolean")
719  public static long saturatedPow(long b, int k) {
720    checkNonNegative("exponent", k);
721    if (b >= -2 & b <= 2) {
722      switch ((int) b) {
723        case 0:
724          return (k == 0) ? 1 : 0;
725        case 1:
726          return 1;
727        case (-1):
728          return ((k & 1) == 0) ? 1 : -1;
729        case 2:
730          if (k >= Long.SIZE - 1) {
731            return Long.MAX_VALUE;
732          }
733          return 1L << k;
734        case (-2):
735          if (k >= Long.SIZE) {
736            return Long.MAX_VALUE + (k & 1);
737          }
738          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
739        default:
740          throw new AssertionError();
741      }
742    }
743    long accum = 1;
744    // if b is negative and k is odd then the limit is MIN otherwise the limit is MAX
745    long limit = Long.MAX_VALUE + ((b >>> (Long.SIZE - 1)) & (k & 1));
746    while (true) {
747      switch (k) {
748        case 0:
749          return accum;
750        case 1:
751          return saturatedMultiply(accum, b);
752        default:
753          if ((k & 1) != 0) {
754            accum = saturatedMultiply(accum, b);
755          }
756          k >>= 1;
757          if (k > 0) {
758            if (-FLOOR_SQRT_MAX_LONG > b | b > FLOOR_SQRT_MAX_LONG) {
759              return limit;
760            }
761            b *= b;
762          }
763      }
764    }
765  }
766
767  @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
768
769  /**
770   * Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
771   * {@code n == 0}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
772   *
773   * @throws IllegalArgumentException if {@code n < 0}
774   */
775  @GwtIncompatible // TODO
776  public static long factorial(int n) {
777    checkNonNegative("n", n);
778    return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
779  }
780
781  static final long[] factorials = {
782    1L,
783    1L,
784    1L * 2,
785    1L * 2 * 3,
786    1L * 2 * 3 * 4,
787    1L * 2 * 3 * 4 * 5,
788    1L * 2 * 3 * 4 * 5 * 6,
789    1L * 2 * 3 * 4 * 5 * 6 * 7,
790    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
791    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
792    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
793    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
794    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
795    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
796    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
797    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
798    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
799    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
800    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
801    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
802    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
803  };
804
805  /**
806   * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
807   * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
808   *
809   * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
810   */
811  public static long binomial(int n, int k) {
812    checkNonNegative("n", n);
813    checkNonNegative("k", k);
814    checkArgument(k <= n, "k (%s) > n (%s)", k, n);
815    if (k > (n >> 1)) {
816      k = n - k;
817    }
818    switch (k) {
819      case 0:
820        return 1;
821      case 1:
822        return n;
823      default:
824        if (n < factorials.length) {
825          return factorials[n] / (factorials[k] * factorials[n - k]);
826        } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
827          return Long.MAX_VALUE;
828        } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
829          // guaranteed not to overflow
830          long result = n--;
831          for (int i = 2; i <= k; n--, i++) {
832            result *= n;
833            result /= i;
834          }
835          return result;
836        } else {
837          int nBits = LongMath.log2(n, RoundingMode.CEILING);
838
839          long result = 1;
840          long numerator = n--;
841          long denominator = 1;
842
843          int numeratorBits = nBits;
844          // This is an upper bound on log2(numerator, ceiling).
845
846          /*
847           * We want to do this in long math for speed, but want to avoid overflow. We adapt the
848           * technique previously used by BigIntegerMath: maintain separate numerator and
849           * denominator accumulators, multiplying the fraction into result when near overflow.
850           */
851          for (int i = 2; i <= k; i++, n--) {
852            if (numeratorBits + nBits < Long.SIZE - 1) {
853              // It's definitely safe to multiply into numerator and denominator.
854              numerator *= n;
855              denominator *= i;
856              numeratorBits += nBits;
857            } else {
858              // It might not be safe to multiply into numerator and denominator,
859              // so multiply (numerator / denominator) into result.
860              result = multiplyFraction(result, numerator, denominator);
861              numerator = n;
862              denominator = i;
863              numeratorBits = nBits;
864            }
865          }
866          return multiplyFraction(result, numerator, denominator);
867        }
868    }
869  }
870
871  /** Returns (x * numerator / denominator), which is assumed to come out to an integral value. */
872  static long multiplyFraction(long x, long numerator, long denominator) {
873    if (x == 1) {
874      return numerator / denominator;
875    }
876    long commonDivisor = gcd(x, denominator);
877    x /= commonDivisor;
878    denominator /= commonDivisor;
879    // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
880    // so denominator must be a divisor of numerator.
881    return x * (numerator / denominator);
882  }
883
884  /*
885   * binomial(biggestBinomials[k], k) fits in a long, but not binomial(biggestBinomials[k] + 1, k).
886   */
887  static final int[] biggestBinomials = {
888    Integer.MAX_VALUE,
889    Integer.MAX_VALUE,
890    Integer.MAX_VALUE,
891    3810779,
892    121977,
893    16175,
894    4337,
895    1733,
896    887,
897    534,
898    361,
899    265,
900    206,
901    169,
902    143,
903    125,
904    111,
905    101,
906    94,
907    88,
908    83,
909    79,
910    76,
911    74,
912    72,
913    70,
914    69,
915    68,
916    67,
917    67,
918    66,
919    66,
920    66,
921    66
922  };
923
924  /*
925   * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, but
926   * binomial(biggestSimpleBinomials[k] + 1, k) does.
927   */
928  @VisibleForTesting
929  static final int[] biggestSimpleBinomials = {
930    Integer.MAX_VALUE,
931    Integer.MAX_VALUE,
932    Integer.MAX_VALUE,
933    2642246,
934    86251,
935    11724,
936    3218,
937    1313,
938    684,
939    419,
940    287,
941    214,
942    169,
943    139,
944    119,
945    105,
946    95,
947    87,
948    81,
949    76,
950    73,
951    70,
952    68,
953    66,
954    64,
955    63,
956    62,
957    62,
958    61,
959    61,
960    61
961  };
962  // These values were generated by using checkedMultiply to see when the simple multiply/divide
963  // algorithm would lead to an overflow.
964
965  static boolean fitsInInt(long x) {
966    return (int) x == x;
967  }
968
969  /**
970   * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward negative infinity. This
971   * method is resilient to overflow.
972   *
973   * @since 14.0
974   */
975  public static long mean(long x, long y) {
976    // Efficient method for computing the arithmetic mean.
977    // The alternative (x + y) / 2 fails for large values.
978    // The alternative (x + y) >>> 1 fails for negative values.
979    return (x & y) + ((x ^ y) >> 1);
980  }
981
982  /*
983   * This bitmask is used as an optimization for cheaply testing for divisibility by 2, 3, or 5.
984   * Each bit is set to 1 for all remainders that indicate divisibility by 2, 3, or 5, so
985   * 1, 7, 11, 13, 17, 19, 23, 29 are set to 0. 30 and up don't matter because they won't be hit.
986   */
987  private static final int SIEVE_30 =
988      ~((1 << 1) | (1 << 7) | (1 << 11) | (1 << 13) | (1 << 17) | (1 << 19) | (1 << 23)
989          | (1 << 29));
990
991  /**
992   * Returns {@code true} if {@code n} is a <a
993   * href="http://mathworld.wolfram.com/PrimeNumber.html">prime number</a>: an integer <i>greater
994   * than one</i> that cannot be factored into a product of <i>smaller</i> positive integers.
995   * Returns {@code false} if {@code n} is zero, one, or a composite number (one which <i>can</i> be
996   * factored into smaller positive integers).
997   *
998   * <p>To test larger numbers, use {@link BigInteger#isProbablePrime}.
999   *
1000   * @throws IllegalArgumentException if {@code n} is negative
1001   * @since 20.0
1002   */
1003  @GwtIncompatible // TODO
1004  public static boolean isPrime(long n) {
1005    if (n < 2) {
1006      checkNonNegative("n", n);
1007      return false;
1008    }
1009    if (n < 66) {
1010      // Encode all primes less than 66 into mask without 0 and 1.
1011      long mask =
1012          (1L << (2 - 2))
1013              | (1L << (3 - 2))
1014              | (1L << (5 - 2))
1015              | (1L << (7 - 2))
1016              | (1L << (11 - 2))
1017              | (1L << (13 - 2))
1018              | (1L << (17 - 2))
1019              | (1L << (19 - 2))
1020              | (1L << (23 - 2))
1021              | (1L << (29 - 2))
1022              | (1L << (31 - 2))
1023              | (1L << (37 - 2))
1024              | (1L << (41 - 2))
1025              | (1L << (43 - 2))
1026              | (1L << (47 - 2))
1027              | (1L << (53 - 2))
1028              | (1L << (59 - 2))
1029              | (1L << (61 - 2));
1030      // Look up n within the mask.
1031      return ((mask >> ((int) n - 2)) & 1) != 0;
1032    }
1033
1034    if ((SIEVE_30 & (1 << (n % 30))) != 0) {
1035      return false;
1036    }
1037    if (n % 7 == 0 || n % 11 == 0 || n % 13 == 0) {
1038      return false;
1039    }
1040    if (n < 17 * 17) {
1041      return true;
1042    }
1043
1044    for (long[] baseSet : millerRabinBaseSets) {
1045      if (n <= baseSet[0]) {
1046        for (int i = 1; i < baseSet.length; i++) {
1047          if (!MillerRabinTester.test(baseSet[i], n)) {
1048            return false;
1049          }
1050        }
1051        return true;
1052      }
1053    }
1054    throw new AssertionError();
1055  }
1056
1057  /*
1058   * If n <= millerRabinBases[i][0], then testing n against bases millerRabinBases[i][1..] suffices
1059   * to prove its primality. Values from miller-rabin.appspot.com.
1060   *
1061   * NOTE: We could get slightly better bases that would be treated as unsigned, but benchmarks
1062   * showed negligible performance improvements.
1063   */
1064  private static final long[][] millerRabinBaseSets = {
1065    {291830, 126401071349994536L},
1066    {885594168, 725270293939359937L, 3569819667048198375L},
1067    {273919523040L, 15, 7363882082L, 992620450144556L},
1068    {47636622961200L, 2, 2570940, 211991001, 3749873356L},
1069    {
1070      7999252175582850L,
1071      2,
1072      4130806001517L,
1073      149795463772692060L,
1074      186635894390467037L,
1075      3967304179347715805L
1076    },
1077    {
1078      585226005592931976L,
1079      2,
1080      123635709730000L,
1081      9233062284813009L,
1082      43835965440333360L,
1083      761179012939631437L,
1084      1263739024124850375L
1085    },
1086    {Long.MAX_VALUE, 2, 325, 9375, 28178, 450775, 9780504, 1795265022}
1087  };
1088
1089  private enum MillerRabinTester {
1090    /** Works for inputs ≤ FLOOR_SQRT_MAX_LONG. */
1091    SMALL {
1092      @Override
1093      long mulMod(long a, long b, long m) {
1094        /*
1095         * lowasser, 2015-Feb-12: Benchmarks suggest that changing this to UnsignedLongs.remainder
1096         * and increasing the threshold to 2^32 doesn't pay for itself, and adding another enum
1097         * constant hurts performance further -- I suspect because bimorphic implementation is a
1098         * sweet spot for the JVM.
1099         */
1100        return (a * b) % m;
1101      }
1102
1103      @Override
1104      long squareMod(long a, long m) {
1105        return (a * a) % m;
1106      }
1107    },
1108    /** Works for all nonnegative signed longs. */
1109    LARGE {
1110      /** Returns (a + b) mod m. Precondition: {@code 0 <= a}, {@code b < m < 2^63}. */
1111      private long plusMod(long a, long b, long m) {
1112        return (a >= m - b) ? (a + b - m) : (a + b);
1113      }
1114
1115      /** Returns (a * 2^32) mod m. a may be any unsigned long. */
1116      private long times2ToThe32Mod(long a, long m) {
1117        int remainingPowersOf2 = 32;
1118        do {
1119          int shift = min(remainingPowersOf2, Long.numberOfLeadingZeros(a));
1120          // shift is either the number of powers of 2 left to multiply a by, or the biggest shift
1121          // possible while keeping a in an unsigned long.
1122          a = UnsignedLongs.remainder(a << shift, m);
1123          remainingPowersOf2 -= shift;
1124        } while (remainingPowersOf2 > 0);
1125        return a;
1126      }
1127
1128      @Override
1129      long mulMod(long a, long b, long m) {
1130        long aHi = a >>> 32; // < 2^31
1131        long bHi = b >>> 32; // < 2^31
1132        long aLo = a & 0xFFFFFFFFL; // < 2^32
1133        long bLo = b & 0xFFFFFFFFL; // < 2^32
1134
1135        /*
1136         * a * b == aHi * bHi * 2^64 + (aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo.
1137         *       == (aHi * bHi * 2^32 + aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo
1138         *
1139         * We carry out this computation in modular arithmetic. Since times2ToThe32Mod accepts any
1140         * unsigned long, we don't have to do a mod on every operation, only when intermediate
1141         * results can exceed 2^63.
1142         */
1143        long result = times2ToThe32Mod(aHi * bHi /* < 2^62 */, m); // < m < 2^63
1144        result += aHi * bLo; // aHi * bLo < 2^63, result < 2^64
1145        if (result < 0) {
1146          result = UnsignedLongs.remainder(result, m);
1147        }
1148        // result < 2^63 again
1149        result += aLo * bHi; // aLo * bHi < 2^63, result < 2^64
1150        result = times2ToThe32Mod(result, m); // result < m < 2^63
1151        return plusMod(result, UnsignedLongs.remainder(aLo * bLo /* < 2^64 */, m), m);
1152      }
1153
1154      @Override
1155      long squareMod(long a, long m) {
1156        long aHi = a >>> 32; // < 2^31
1157        long aLo = a & 0xFFFFFFFFL; // < 2^32
1158
1159        /*
1160         * a^2 == aHi^2 * 2^64 + aHi * aLo * 2^33 + aLo^2
1161         *     == (aHi^2 * 2^32 + aHi * aLo * 2) * 2^32 + aLo^2
1162         * We carry out this computation in modular arithmetic.  Since times2ToThe32Mod accepts any
1163         * unsigned long, we don't have to do a mod on every operation, only when intermediate
1164         * results can exceed 2^63.
1165         */
1166        long result = times2ToThe32Mod(aHi * aHi /* < 2^62 */, m); // < m < 2^63
1167        long hiLo = aHi * aLo * 2;
1168        if (hiLo < 0) {
1169          hiLo = UnsignedLongs.remainder(hiLo, m);
1170        }
1171        // hiLo < 2^63
1172        result += hiLo; // result < 2^64
1173        result = times2ToThe32Mod(result, m); // result < m < 2^63
1174        return plusMod(result, UnsignedLongs.remainder(aLo * aLo /* < 2^64 */, m), m);
1175      }
1176    };
1177
1178    static boolean test(long base, long n) {
1179      // Since base will be considered % n, it's okay if base > FLOOR_SQRT_MAX_LONG,
1180      // so long as n <= FLOOR_SQRT_MAX_LONG.
1181      return ((n <= FLOOR_SQRT_MAX_LONG) ? SMALL : LARGE).testWitness(base, n);
1182    }
1183
1184    /** Returns a * b mod m. */
1185    abstract long mulMod(long a, long b, long m);
1186
1187    /** Returns a^2 mod m. */
1188    abstract long squareMod(long a, long m);
1189
1190    /** Returns a^p mod m. */
1191    private long powMod(long a, long p, long m) {
1192      long res = 1;
1193      for (; p != 0; p >>= 1) {
1194        if ((p & 1) != 0) {
1195          res = mulMod(res, a, m);
1196        }
1197        a = squareMod(a, m);
1198      }
1199      return res;
1200    }
1201
1202    /** Returns true if n is a strong probable prime relative to the specified base. */
1203    private boolean testWitness(long base, long n) {
1204      int r = Long.numberOfTrailingZeros(n - 1);
1205      long d = (n - 1) >> r;
1206      base %= n;
1207      if (base == 0) {
1208        return true;
1209      }
1210      // Calculate a := base^d mod n.
1211      long a = powMod(base, d, n);
1212      // n passes this test if
1213      //    base^d = 1 (mod n)
1214      // or base^(2^j * d) = -1 (mod n) for some 0 <= j < r.
1215      if (a == 1) {
1216        return true;
1217      }
1218      int j = 0;
1219      while (a != n - 1) {
1220        if (++j == r) {
1221          return false;
1222        }
1223        a = squareMod(a, n);
1224      }
1225      return true;
1226    }
1227  }
1228
1229  /**
1230   * Returns {@code x}, rounded to a {@code double} with the specified rounding mode. If {@code x}
1231   * is precisely representable as a {@code double}, its {@code double} value will be returned;
1232   * otherwise, the rounding will choose between the two nearest representable values with {@code
1233   * mode}.
1234   *
1235   * <p>For the case of {@link RoundingMode#HALF_EVEN}, this implementation uses the IEEE 754
1236   * default rounding mode: if the two nearest representable values are equally near, the one with
1237   * the least significant bit zero is chosen. (In such cases, both of the nearest representable
1238   * values are even integers; this method returns the one that is a multiple of a greater power of
1239   * two.)
1240   *
1241   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
1242   *     is not precisely representable as a {@code double}
1243   * @since 30.0
1244   */
1245  @SuppressWarnings("deprecation")
1246  @GwtIncompatible
1247  public static double roundToDouble(long x, RoundingMode mode) {
1248    // Logic adapted from ToDoubleRounder.
1249    double roundArbitrarily = (double) x;
1250    long roundArbitrarilyAsLong = (long) roundArbitrarily;
1251    int cmpXToRoundArbitrarily;
1252
1253    if (roundArbitrarilyAsLong == Long.MAX_VALUE) {
1254      /*
1255       * For most values, the conversion from roundArbitrarily to roundArbitrarilyAsLong is
1256       * lossless. In that case we can compare x to roundArbitrarily using Long.compare(x,
1257       * roundArbitrarilyAsLong). The exception is for values where the conversion to double rounds
1258       * up to give roundArbitrarily equal to 2^63, so the conversion back to long overflows and
1259       * roundArbitrarilyAsLong is Long.MAX_VALUE. (This is the only way this condition can occur as
1260       * otherwise the conversion back to long pads with zero bits.) In this case we know that
1261       * roundArbitrarily > x. (This is important when x == Long.MAX_VALUE ==
1262       * roundArbitrarilyAsLong.)
1263       */
1264      cmpXToRoundArbitrarily = -1;
1265    } else {
1266      cmpXToRoundArbitrarily = Long.compare(x, roundArbitrarilyAsLong);
1267    }
1268
1269    switch (mode) {
1270      case UNNECESSARY:
1271        checkRoundingUnnecessary(cmpXToRoundArbitrarily == 0);
1272        return roundArbitrarily;
1273      case FLOOR:
1274        return (cmpXToRoundArbitrarily >= 0)
1275            ? roundArbitrarily
1276            : DoubleUtils.nextDown(roundArbitrarily);
1277      case CEILING:
1278        return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1279      case DOWN:
1280        if (x >= 0) {
1281          return (cmpXToRoundArbitrarily >= 0)
1282              ? roundArbitrarily
1283              : DoubleUtils.nextDown(roundArbitrarily);
1284        } else {
1285          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1286        }
1287      case UP:
1288        if (x >= 0) {
1289          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1290        } else {
1291          return (cmpXToRoundArbitrarily >= 0)
1292              ? roundArbitrarily
1293              : DoubleUtils.nextDown(roundArbitrarily);
1294        }
1295      case HALF_DOWN:
1296      case HALF_UP:
1297      case HALF_EVEN:
1298        {
1299          long roundFloor;
1300          double roundFloorAsDouble;
1301          long roundCeiling;
1302          double roundCeilingAsDouble;
1303
1304          if (cmpXToRoundArbitrarily >= 0) {
1305            roundFloorAsDouble = roundArbitrarily;
1306            roundFloor = roundArbitrarilyAsLong;
1307            roundCeilingAsDouble = Math.nextUp(roundArbitrarily);
1308            roundCeiling = (long) Math.ceil(roundCeilingAsDouble);
1309          } else {
1310            roundCeilingAsDouble = roundArbitrarily;
1311            roundCeiling = roundArbitrarilyAsLong;
1312            roundFloorAsDouble = DoubleUtils.nextDown(roundArbitrarily);
1313            roundFloor = (long) Math.floor(roundFloorAsDouble);
1314          }
1315
1316          long deltaToFloor = x - roundFloor;
1317          long deltaToCeiling = roundCeiling - x;
1318
1319          if (roundCeiling == Long.MAX_VALUE) {
1320            // correct for Long.MAX_VALUE as discussed above: roundCeilingAsDouble must be 2^63, but
1321            // roundCeiling is 2^63-1.
1322            deltaToCeiling++;
1323          }
1324
1325          int diff = Long.compare(deltaToFloor, deltaToCeiling);
1326          if (diff < 0) { // closer to floor
1327            return roundFloorAsDouble;
1328          } else if (diff > 0) { // closer to ceiling
1329            return roundCeilingAsDouble;
1330          }
1331          // halfway between the representable values; do the half-whatever logic
1332          switch (mode) {
1333            case HALF_EVEN:
1334              return ((DoubleUtils.getSignificand(roundFloorAsDouble) & 1L) == 0)
1335                  ? roundFloorAsDouble
1336                  : roundCeilingAsDouble;
1337            case HALF_DOWN:
1338              return (x >= 0) ? roundFloorAsDouble : roundCeilingAsDouble;
1339            case HALF_UP:
1340              return (x >= 0) ? roundCeilingAsDouble : roundFloorAsDouble;
1341            default:
1342              throw new AssertionError("impossible");
1343          }
1344        }
1345    }
1346    throw new AssertionError("impossible");
1347  }
1348
1349  private LongMath() {}
1350}