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