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  @GwtIncompatible // TODO
551  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
552  @SuppressWarnings("ShortCircuitBoolean")
553  public static long checkedSubtract(long a, long b) {
554    long result = a - b;
555    checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0, "checkedSubtract", a, b);
556    return result;
557  }
558
559  /**
560   * Returns the product of {@code a} and {@code b}, provided it does not overflow.
561   *
562   * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
563   */
564  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
565  @SuppressWarnings("ShortCircuitBoolean")
566  public static long checkedMultiply(long a, long b) {
567    // Hacker's Delight, Section 2-12
568    int leadingZeros =
569        Long.numberOfLeadingZeros(a)
570            + Long.numberOfLeadingZeros(~a)
571            + Long.numberOfLeadingZeros(b)
572            + Long.numberOfLeadingZeros(~b);
573    /*
574     * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
575     * bad. We do the leadingZeros check to avoid the division below if at all possible.
576     *
577     * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
578     * care of all a < 0 with their own check, because in particular, the case a == -1 will
579     * incorrectly pass the division check below.
580     *
581     * In all other cases, we check that either a is 0 or the result is consistent with division.
582     */
583    if (leadingZeros > Long.SIZE + 1) {
584      return a * b;
585    }
586    checkNoOverflow(leadingZeros >= Long.SIZE, "checkedMultiply", a, b);
587    checkNoOverflow(a >= 0 | b != Long.MIN_VALUE, "checkedMultiply", a, b);
588    long result = a * b;
589    checkNoOverflow(a == 0 || result / a == b, "checkedMultiply", a, b);
590    return result;
591  }
592
593  /**
594   * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
595   *
596   * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed {@code
597   *     long} arithmetic
598   */
599  @GwtIncompatible // TODO
600  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
601  @SuppressWarnings("ShortCircuitBoolean")
602  public static long checkedPow(long b, int k) {
603    checkNonNegative("exponent", k);
604    if (b >= -2 & b <= 2) {
605      switch ((int) b) {
606        case 0:
607          return (k == 0) ? 1 : 0;
608        case 1:
609          return 1;
610        case (-1):
611          return ((k & 1) == 0) ? 1 : -1;
612        case 2:
613          checkNoOverflow(k < Long.SIZE - 1, "checkedPow", b, k);
614          return 1L << k;
615        case (-2):
616          checkNoOverflow(k < Long.SIZE, "checkedPow", b, k);
617          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
618        default:
619          throw new AssertionError();
620      }
621    }
622    long accum = 1;
623    while (true) {
624      switch (k) {
625        case 0:
626          return accum;
627        case 1:
628          return checkedMultiply(accum, b);
629        default:
630          if ((k & 1) != 0) {
631            accum = checkedMultiply(accum, b);
632          }
633          k >>= 1;
634          if (k > 0) {
635            checkNoOverflow(
636                -FLOOR_SQRT_MAX_LONG <= b && b <= FLOOR_SQRT_MAX_LONG, "checkedPow", b, k);
637            b *= b;
638          }
639      }
640    }
641  }
642
643  /**
644   * Returns the sum of {@code a} and {@code b} unless it would overflow or underflow in which case
645   * {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
646   *
647   * @since 20.0
648   */
649  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
650  @SuppressWarnings("ShortCircuitBoolean")
651  public static long saturatedAdd(long a, long b) {
652    long naiveSum = a + b;
653    if ((a ^ b) < 0 | (a ^ naiveSum) >= 0) {
654      // If a and b have different signs or a has the same sign as the result then there was no
655      // overflow, return.
656      return naiveSum;
657    }
658    // we did over/under flow, if the sign is negative we should return MAX otherwise MIN
659    return Long.MAX_VALUE + ((naiveSum >>> (Long.SIZE - 1)) ^ 1);
660  }
661
662  /**
663   * Returns the difference of {@code a} and {@code b} unless it would overflow or underflow in
664   * which case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
665   *
666   * @since 20.0
667   */
668  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
669  @SuppressWarnings("ShortCircuitBoolean")
670  public static long saturatedSubtract(long a, long b) {
671    long naiveDifference = a - b;
672    if ((a ^ b) >= 0 | (a ^ naiveDifference) >= 0) {
673      // If a and b have the same signs or a has the same sign as the result then there was no
674      // overflow, return.
675      return naiveDifference;
676    }
677    // we did over/under flow
678    return Long.MAX_VALUE + ((naiveDifference >>> (Long.SIZE - 1)) ^ 1);
679  }
680
681  /**
682   * Returns the product of {@code a} and {@code b} unless it would overflow or underflow in which
683   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
684   *
685   * @since 20.0
686   */
687  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
688  @SuppressWarnings("ShortCircuitBoolean")
689  public static long saturatedMultiply(long a, long b) {
690    // see checkedMultiply for explanation
691    int leadingZeros =
692        Long.numberOfLeadingZeros(a)
693            + Long.numberOfLeadingZeros(~a)
694            + Long.numberOfLeadingZeros(b)
695            + Long.numberOfLeadingZeros(~b);
696    if (leadingZeros > Long.SIZE + 1) {
697      return a * b;
698    }
699    // the return value if we will overflow (which we calculate by overflowing a long :) )
700    long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
701    if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
702      // overflow
703      return limit;
704    }
705    long result = a * b;
706    if (a == 0 || result / a == b) {
707      return result;
708    }
709    return limit;
710  }
711
712  /**
713   * Returns the {@code b} to the {@code k}th power, unless it would overflow or underflow in which
714   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
715   *
716   * @since 20.0
717   */
718  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
719  @SuppressWarnings("ShortCircuitBoolean")
720  public static long saturatedPow(long b, int k) {
721    checkNonNegative("exponent", k);
722    if (b >= -2 & b <= 2) {
723      switch ((int) b) {
724        case 0:
725          return (k == 0) ? 1 : 0;
726        case 1:
727          return 1;
728        case (-1):
729          return ((k & 1) == 0) ? 1 : -1;
730        case 2:
731          if (k >= Long.SIZE - 1) {
732            return Long.MAX_VALUE;
733          }
734          return 1L << k;
735        case (-2):
736          if (k >= Long.SIZE) {
737            return Long.MAX_VALUE + (k & 1);
738          }
739          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
740        default:
741          throw new AssertionError();
742      }
743    }
744    long accum = 1;
745    // if b is negative and k is odd then the limit is MIN otherwise the limit is MAX
746    long limit = Long.MAX_VALUE + ((b >>> (Long.SIZE - 1)) & (k & 1));
747    while (true) {
748      switch (k) {
749        case 0:
750          return accum;
751        case 1:
752          return saturatedMultiply(accum, b);
753        default:
754          if ((k & 1) != 0) {
755            accum = saturatedMultiply(accum, b);
756          }
757          k >>= 1;
758          if (k > 0) {
759            if (-FLOOR_SQRT_MAX_LONG > b | b > FLOOR_SQRT_MAX_LONG) {
760              return limit;
761            }
762            b *= b;
763          }
764      }
765    }
766  }
767
768  @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
769
770  /**
771   * Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
772   * {@code n == 0}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
773   *
774   * @throws IllegalArgumentException if {@code n < 0}
775   */
776  @GwtIncompatible // TODO
777  public static long factorial(int n) {
778    checkNonNegative("n", n);
779    return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
780  }
781
782  static final long[] factorials = {
783    1L,
784    1L,
785    1L * 2,
786    1L * 2 * 3,
787    1L * 2 * 3 * 4,
788    1L * 2 * 3 * 4 * 5,
789    1L * 2 * 3 * 4 * 5 * 6,
790    1L * 2 * 3 * 4 * 5 * 6 * 7,
791    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
792    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
793    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
794    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
795    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
796    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
797    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
798    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
799    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
800    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
801    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
802    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
803    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
804  };
805
806  /**
807   * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
808   * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
809   *
810   * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
811   */
812  public static long binomial(int n, int k) {
813    checkNonNegative("n", n);
814    checkNonNegative("k", k);
815    checkArgument(k <= n, "k (%s) > n (%s)", k, n);
816    if (k > (n >> 1)) {
817      k = n - k;
818    }
819    switch (k) {
820      case 0:
821        return 1;
822      case 1:
823        return n;
824      default:
825        if (n < factorials.length) {
826          return factorials[n] / (factorials[k] * factorials[n - k]);
827        } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
828          return Long.MAX_VALUE;
829        } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
830          // guaranteed not to overflow
831          long result = n--;
832          for (int i = 2; i <= k; n--, i++) {
833            result *= n;
834            result /= i;
835          }
836          return result;
837        } else {
838          int nBits = LongMath.log2(n, RoundingMode.CEILING);
839
840          long result = 1;
841          long numerator = n--;
842          long denominator = 1;
843
844          int numeratorBits = nBits;
845          // This is an upper bound on log2(numerator, ceiling).
846
847          /*
848           * We want to do this in long math for speed, but want to avoid overflow. We adapt the
849           * technique previously used by BigIntegerMath: maintain separate numerator and
850           * denominator accumulators, multiplying the fraction into result when near overflow.
851           */
852          for (int i = 2; i <= k; i++, n--) {
853            if (numeratorBits + nBits < Long.SIZE - 1) {
854              // It's definitely safe to multiply into numerator and denominator.
855              numerator *= n;
856              denominator *= i;
857              numeratorBits += nBits;
858            } else {
859              // It might not be safe to multiply into numerator and denominator,
860              // so multiply (numerator / denominator) into result.
861              result = multiplyFraction(result, numerator, denominator);
862              numerator = n;
863              denominator = i;
864              numeratorBits = nBits;
865            }
866          }
867          return multiplyFraction(result, numerator, denominator);
868        }
869    }
870  }
871
872  /** Returns (x * numerator / denominator), which is assumed to come out to an integral value. */
873  static long multiplyFraction(long x, long numerator, long denominator) {
874    if (x == 1) {
875      return numerator / denominator;
876    }
877    long commonDivisor = gcd(x, denominator);
878    x /= commonDivisor;
879    denominator /= commonDivisor;
880    // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
881    // so denominator must be a divisor of numerator.
882    return x * (numerator / denominator);
883  }
884
885  /*
886   * binomial(biggestBinomials[k], k) fits in a long, but not binomial(biggestBinomials[k] + 1, k).
887   */
888  static final int[] biggestBinomials = {
889    Integer.MAX_VALUE,
890    Integer.MAX_VALUE,
891    Integer.MAX_VALUE,
892    3810779,
893    121977,
894    16175,
895    4337,
896    1733,
897    887,
898    534,
899    361,
900    265,
901    206,
902    169,
903    143,
904    125,
905    111,
906    101,
907    94,
908    88,
909    83,
910    79,
911    76,
912    74,
913    72,
914    70,
915    69,
916    68,
917    67,
918    67,
919    66,
920    66,
921    66,
922    66
923  };
924
925  /*
926   * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, but
927   * binomial(biggestSimpleBinomials[k] + 1, k) does.
928   */
929  @VisibleForTesting
930  static final int[] biggestSimpleBinomials = {
931    Integer.MAX_VALUE,
932    Integer.MAX_VALUE,
933    Integer.MAX_VALUE,
934    2642246,
935    86251,
936    11724,
937    3218,
938    1313,
939    684,
940    419,
941    287,
942    214,
943    169,
944    139,
945    119,
946    105,
947    95,
948    87,
949    81,
950    76,
951    73,
952    70,
953    68,
954    66,
955    64,
956    63,
957    62,
958    62,
959    61,
960    61,
961    61
962  };
963  // These values were generated by using checkedMultiply to see when the simple multiply/divide
964  // algorithm would lead to an overflow.
965
966  static boolean fitsInInt(long x) {
967    return (int) x == x;
968  }
969
970  /**
971   * Returns the arithmetic mean of {@code x} and {@code y}, rounded toward negative infinity. This
972   * method is resilient to overflow.
973   *
974   * @since 14.0
975   */
976  public static long mean(long x, long y) {
977    // Efficient method for computing the arithmetic mean.
978    // The alternative (x + y) / 2 fails for large values.
979    // The alternative (x + y) >>> 1 fails for negative values.
980    return (x & y) + ((x ^ y) >> 1);
981  }
982
983  /*
984   * This bitmask is used as an optimization for cheaply testing for divisibility by 2, 3, or 5.
985   * Each bit is set to 1 for all remainders that indicate divisibility by 2, 3, or 5, so
986   * 1, 7, 11, 13, 17, 19, 23, 29 are set to 0. 30 and up don't matter because they won't be hit.
987   */
988  private static final int SIEVE_30 =
989      ~((1 << 1) | (1 << 7) | (1 << 11) | (1 << 13) | (1 << 17) | (1 << 19) | (1 << 23)
990          | (1 << 29));
991
992  /**
993   * Returns {@code true} if {@code n} is a <a
994   * href="http://mathworld.wolfram.com/PrimeNumber.html">prime number</a>: an integer <i>greater
995   * than one</i> that cannot be factored into a product of <i>smaller</i> positive integers.
996   * Returns {@code false} if {@code n} is zero, one, or a composite number (one which <i>can</i> be
997   * factored into smaller positive integers).
998   *
999   * <p>To test larger numbers, use {@link BigInteger#isProbablePrime}.
1000   *
1001   * @throws IllegalArgumentException if {@code n} is negative
1002   * @since 20.0
1003   */
1004  @GwtIncompatible // TODO
1005  public static boolean isPrime(long n) {
1006    if (n < 2) {
1007      checkNonNegative("n", n);
1008      return false;
1009    }
1010    if (n < 66) {
1011      // Encode all primes less than 66 into mask without 0 and 1.
1012      long mask =
1013          (1L << (2 - 2))
1014              | (1L << (3 - 2))
1015              | (1L << (5 - 2))
1016              | (1L << (7 - 2))
1017              | (1L << (11 - 2))
1018              | (1L << (13 - 2))
1019              | (1L << (17 - 2))
1020              | (1L << (19 - 2))
1021              | (1L << (23 - 2))
1022              | (1L << (29 - 2))
1023              | (1L << (31 - 2))
1024              | (1L << (37 - 2))
1025              | (1L << (41 - 2))
1026              | (1L << (43 - 2))
1027              | (1L << (47 - 2))
1028              | (1L << (53 - 2))
1029              | (1L << (59 - 2))
1030              | (1L << (61 - 2));
1031      // Look up n within the mask.
1032      return ((mask >> ((int) n - 2)) & 1) != 0;
1033    }
1034
1035    if ((SIEVE_30 & (1 << (n % 30))) != 0) {
1036      return false;
1037    }
1038    if (n % 7 == 0 || n % 11 == 0 || n % 13 == 0) {
1039      return false;
1040    }
1041    if (n < 17 * 17) {
1042      return true;
1043    }
1044
1045    for (long[] baseSet : millerRabinBaseSets) {
1046      if (n <= baseSet[0]) {
1047        for (int i = 1; i < baseSet.length; i++) {
1048          if (!MillerRabinTester.test(baseSet[i], n)) {
1049            return false;
1050          }
1051        }
1052        return true;
1053      }
1054    }
1055    throw new AssertionError();
1056  }
1057
1058  /*
1059   * If n <= millerRabinBases[i][0], then testing n against bases millerRabinBases[i][1..] suffices
1060   * to prove its primality. Values from miller-rabin.appspot.com.
1061   *
1062   * NOTE: We could get slightly better bases that would be treated as unsigned, but benchmarks
1063   * showed negligible performance improvements.
1064   */
1065  private static final long[][] millerRabinBaseSets = {
1066    {291830, 126401071349994536L},
1067    {885594168, 725270293939359937L, 3569819667048198375L},
1068    {273919523040L, 15, 7363882082L, 992620450144556L},
1069    {47636622961200L, 2, 2570940, 211991001, 3749873356L},
1070    {
1071      7999252175582850L,
1072      2,
1073      4130806001517L,
1074      149795463772692060L,
1075      186635894390467037L,
1076      3967304179347715805L
1077    },
1078    {
1079      585226005592931976L,
1080      2,
1081      123635709730000L,
1082      9233062284813009L,
1083      43835965440333360L,
1084      761179012939631437L,
1085      1263739024124850375L
1086    },
1087    {Long.MAX_VALUE, 2, 325, 9375, 28178, 450775, 9780504, 1795265022}
1088  };
1089
1090  private enum MillerRabinTester {
1091    /** Works for inputs ≤ FLOOR_SQRT_MAX_LONG. */
1092    SMALL {
1093      @Override
1094      long mulMod(long a, long b, long m) {
1095        /*
1096         * lowasser, 2015-Feb-12: Benchmarks suggest that changing this to UnsignedLongs.remainder
1097         * and increasing the threshold to 2^32 doesn't pay for itself, and adding another enum
1098         * constant hurts performance further -- I suspect because bimorphic implementation is a
1099         * sweet spot for the JVM.
1100         */
1101        return (a * b) % m;
1102      }
1103
1104      @Override
1105      long squareMod(long a, long m) {
1106        return (a * a) % m;
1107      }
1108    },
1109    /** Works for all nonnegative signed longs. */
1110    LARGE {
1111      /** Returns (a + b) mod m. Precondition: {@code 0 <= a}, {@code b < m < 2^63}. */
1112      private long plusMod(long a, long b, long m) {
1113        return (a >= m - b) ? (a + b - m) : (a + b);
1114      }
1115
1116      /** Returns (a * 2^32) mod m. a may be any unsigned long. */
1117      private long times2ToThe32Mod(long a, long m) {
1118        int remainingPowersOf2 = 32;
1119        do {
1120          int shift = min(remainingPowersOf2, Long.numberOfLeadingZeros(a));
1121          // shift is either the number of powers of 2 left to multiply a by, or the biggest shift
1122          // possible while keeping a in an unsigned long.
1123          a = UnsignedLongs.remainder(a << shift, m);
1124          remainingPowersOf2 -= shift;
1125        } while (remainingPowersOf2 > 0);
1126        return a;
1127      }
1128
1129      @Override
1130      long mulMod(long a, long b, long m) {
1131        long aHi = a >>> 32; // < 2^31
1132        long bHi = b >>> 32; // < 2^31
1133        long aLo = a & 0xFFFFFFFFL; // < 2^32
1134        long bLo = b & 0xFFFFFFFFL; // < 2^32
1135
1136        /*
1137         * a * b == aHi * bHi * 2^64 + (aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo.
1138         *       == (aHi * bHi * 2^32 + aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo
1139         *
1140         * We carry out this computation in modular arithmetic. Since times2ToThe32Mod accepts any
1141         * unsigned long, we don't have to do a mod on every operation, only when intermediate
1142         * results can exceed 2^63.
1143         */
1144        long result = times2ToThe32Mod(aHi * bHi /* < 2^62 */, m); // < m < 2^63
1145        result += aHi * bLo; // aHi * bLo < 2^63, result < 2^64
1146        if (result < 0) {
1147          result = UnsignedLongs.remainder(result, m);
1148        }
1149        // result < 2^63 again
1150        result += aLo * bHi; // aLo * bHi < 2^63, result < 2^64
1151        result = times2ToThe32Mod(result, m); // result < m < 2^63
1152        return plusMod(result, UnsignedLongs.remainder(aLo * bLo /* < 2^64 */, m), m);
1153      }
1154
1155      @Override
1156      long squareMod(long a, long m) {
1157        long aHi = a >>> 32; // < 2^31
1158        long aLo = a & 0xFFFFFFFFL; // < 2^32
1159
1160        /*
1161         * a^2 == aHi^2 * 2^64 + aHi * aLo * 2^33 + aLo^2
1162         *     == (aHi^2 * 2^32 + aHi * aLo * 2) * 2^32 + aLo^2
1163         * We carry out this computation in modular arithmetic.  Since times2ToThe32Mod accepts any
1164         * unsigned long, we don't have to do a mod on every operation, only when intermediate
1165         * results can exceed 2^63.
1166         */
1167        long result = times2ToThe32Mod(aHi * aHi /* < 2^62 */, m); // < m < 2^63
1168        long hiLo = aHi * aLo * 2;
1169        if (hiLo < 0) {
1170          hiLo = UnsignedLongs.remainder(hiLo, m);
1171        }
1172        // hiLo < 2^63
1173        result += hiLo; // result < 2^64
1174        result = times2ToThe32Mod(result, m); // result < m < 2^63
1175        return plusMod(result, UnsignedLongs.remainder(aLo * aLo /* < 2^64 */, m), m);
1176      }
1177    };
1178
1179    static boolean test(long base, long n) {
1180      // Since base will be considered % n, it's okay if base > FLOOR_SQRT_MAX_LONG,
1181      // so long as n <= FLOOR_SQRT_MAX_LONG.
1182      return ((n <= FLOOR_SQRT_MAX_LONG) ? SMALL : LARGE).testWitness(base, n);
1183    }
1184
1185    /** Returns a * b mod m. */
1186    abstract long mulMod(long a, long b, long m);
1187
1188    /** Returns a^2 mod m. */
1189    abstract long squareMod(long a, long m);
1190
1191    /** Returns a^p mod m. */
1192    private long powMod(long a, long p, long m) {
1193      long res = 1;
1194      for (; p != 0; p >>= 1) {
1195        if ((p & 1) != 0) {
1196          res = mulMod(res, a, m);
1197        }
1198        a = squareMod(a, m);
1199      }
1200      return res;
1201    }
1202
1203    /** Returns true if n is a strong probable prime relative to the specified base. */
1204    private boolean testWitness(long base, long n) {
1205      int r = Long.numberOfTrailingZeros(n - 1);
1206      long d = (n - 1) >> r;
1207      base %= n;
1208      if (base == 0) {
1209        return true;
1210      }
1211      // Calculate a := base^d mod n.
1212      long a = powMod(base, d, n);
1213      // n passes this test if
1214      //    base^d = 1 (mod n)
1215      // or base^(2^j * d) = -1 (mod n) for some 0 <= j < r.
1216      if (a == 1) {
1217        return true;
1218      }
1219      int j = 0;
1220      while (a != n - 1) {
1221        if (++j == r) {
1222          return false;
1223        }
1224        a = squareMod(a, n);
1225      }
1226      return true;
1227    }
1228  }
1229
1230  /**
1231   * Returns {@code x}, rounded to a {@code double} with the specified rounding mode. If {@code x}
1232   * is precisely representable as a {@code double}, its {@code double} value will be returned;
1233   * otherwise, the rounding will choose between the two nearest representable values with {@code
1234   * mode}.
1235   *
1236   * <p>For the case of {@link RoundingMode#HALF_EVEN}, this implementation uses the IEEE 754
1237   * default rounding mode: if the two nearest representable values are equally near, the one with
1238   * the least significant bit zero is chosen. (In such cases, both of the nearest representable
1239   * values are even integers; this method returns the one that is a multiple of a greater power of
1240   * two.)
1241   *
1242   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
1243   *     is not precisely representable as a {@code double}
1244   * @since 30.0
1245   */
1246  @SuppressWarnings("deprecation")
1247  @GwtIncompatible
1248  public static double roundToDouble(long x, RoundingMode mode) {
1249    // Logic adapted from ToDoubleRounder.
1250    double roundArbitrarily = (double) x;
1251    long roundArbitrarilyAsLong = (long) roundArbitrarily;
1252    int cmpXToRoundArbitrarily;
1253
1254    if (roundArbitrarilyAsLong == Long.MAX_VALUE) {
1255      /*
1256       * For most values, the conversion from roundArbitrarily to roundArbitrarilyAsLong is
1257       * lossless. In that case we can compare x to roundArbitrarily using Long.compare(x,
1258       * roundArbitrarilyAsLong). The exception is for values where the conversion to double rounds
1259       * up to give roundArbitrarily equal to 2^63, so the conversion back to long overflows and
1260       * roundArbitrarilyAsLong is Long.MAX_VALUE. (This is the only way this condition can occur as
1261       * otherwise the conversion back to long pads with zero bits.) In this case we know that
1262       * roundArbitrarily > x. (This is important when x == Long.MAX_VALUE ==
1263       * roundArbitrarilyAsLong.)
1264       */
1265      cmpXToRoundArbitrarily = -1;
1266    } else {
1267      cmpXToRoundArbitrarily = Long.compare(x, roundArbitrarilyAsLong);
1268    }
1269
1270    switch (mode) {
1271      case UNNECESSARY:
1272        checkRoundingUnnecessary(cmpXToRoundArbitrarily == 0);
1273        return roundArbitrarily;
1274      case FLOOR:
1275        return (cmpXToRoundArbitrarily >= 0)
1276            ? roundArbitrarily
1277            : DoubleUtils.nextDown(roundArbitrarily);
1278      case CEILING:
1279        return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1280      case DOWN:
1281        if (x >= 0) {
1282          return (cmpXToRoundArbitrarily >= 0)
1283              ? roundArbitrarily
1284              : DoubleUtils.nextDown(roundArbitrarily);
1285        } else {
1286          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1287        }
1288      case UP:
1289        if (x >= 0) {
1290          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1291        } else {
1292          return (cmpXToRoundArbitrarily >= 0)
1293              ? roundArbitrarily
1294              : DoubleUtils.nextDown(roundArbitrarily);
1295        }
1296      case HALF_DOWN:
1297      case HALF_UP:
1298      case HALF_EVEN:
1299        {
1300          long roundFloor;
1301          double roundFloorAsDouble;
1302          long roundCeiling;
1303          double roundCeilingAsDouble;
1304
1305          if (cmpXToRoundArbitrarily >= 0) {
1306            roundFloorAsDouble = roundArbitrarily;
1307            roundFloor = roundArbitrarilyAsLong;
1308            roundCeilingAsDouble = Math.nextUp(roundArbitrarily);
1309            roundCeiling = (long) Math.ceil(roundCeilingAsDouble);
1310          } else {
1311            roundCeilingAsDouble = roundArbitrarily;
1312            roundCeiling = roundArbitrarilyAsLong;
1313            roundFloorAsDouble = DoubleUtils.nextDown(roundArbitrarily);
1314            roundFloor = (long) Math.floor(roundFloorAsDouble);
1315          }
1316
1317          long deltaToFloor = x - roundFloor;
1318          long deltaToCeiling = roundCeiling - x;
1319
1320          if (roundCeiling == Long.MAX_VALUE) {
1321            // correct for Long.MAX_VALUE as discussed above: roundCeilingAsDouble must be 2^63, but
1322            // roundCeiling is 2^63-1.
1323            deltaToCeiling++;
1324          }
1325
1326          int diff = Long.compare(deltaToFloor, deltaToCeiling);
1327          if (diff < 0) { // closer to floor
1328            return roundFloorAsDouble;
1329          } else if (diff > 0) { // closer to ceiling
1330            return roundCeilingAsDouble;
1331          }
1332          // halfway between the representable values; do the half-whatever logic
1333          switch (mode) {
1334            case HALF_EVEN:
1335              return ((DoubleUtils.getSignificand(roundFloorAsDouble) & 1L) == 0)
1336                  ? roundFloorAsDouble
1337                  : roundCeilingAsDouble;
1338            case HALF_DOWN:
1339              return (x >= 0) ? roundFloorAsDouble : roundCeilingAsDouble;
1340            case HALF_UP:
1341              return (x >= 0) ? roundCeilingAsDouble : roundFloorAsDouble;
1342            default:
1343              throw new AssertionError("impossible");
1344          }
1345        }
1346    }
1347    throw new AssertionError("impossible");
1348  }
1349
1350  private LongMath() {}
1351}