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.ceil;
025import static java.lang.Math.floor;
026import static java.lang.Math.min;
027import static java.lang.Math.nextDown;
028import static java.lang.Math.nextUp;
029import static java.math.RoundingMode.HALF_EVEN;
030import static java.math.RoundingMode.HALF_UP;
031
032import com.google.common.annotations.GwtCompatible;
033import com.google.common.annotations.GwtIncompatible;
034import com.google.common.annotations.VisibleForTesting;
035import com.google.common.primitives.UnsignedLongs;
036import com.google.errorprone.annotations.InlineMe;
037import java.math.BigInteger;
038import java.math.RoundingMode;
039
040/**
041 * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
042 * named analogously to their {@code BigInteger} counterparts.
043 *
044 * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
045 * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
046 *
047 * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in {@link
048 * IntMath} and {@link BigIntegerMath} respectively. For other common operations on {@code long}
049 * values, see {@link com.google.common.primitives.Longs}.
050 *
051 * @author Louis Wasserman
052 * @since 11.0
053 */
054@GwtCompatible
055public final class LongMath {
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  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
097  @SuppressWarnings("ShortCircuitBoolean")
098  public static boolean isPowerOfTwo(long x) {
099    return x > 0 & (x & (x - 1)) == 0;
100  }
101
102  /**
103   * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
104   * signed long. The implementation is branch-free, and benchmarks suggest it is measurably faster
105   * than the straightforward ternary expression.
106   */
107  @VisibleForTesting
108  static int lessThanBranchFree(long x, long y) {
109    // The double negation is optimized away by normal Java, but is necessary for GWT
110    // to make sure bit twiddling works as expected.
111    return (int) (~~(x - y) >>> (Long.SIZE - 1));
112  }
113
114  /**
115   * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
116   *
117   * @throws IllegalArgumentException if {@code x <= 0}
118   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
119   *     is not a power of two
120   */
121  @SuppressWarnings("fallthrough")
122  // TODO(kevinb): remove after this warning is disabled globally
123  public static int log2(long x, RoundingMode mode) {
124    checkPositive("x", x);
125    switch (mode) {
126      case UNNECESSARY:
127        checkRoundingUnnecessary(isPowerOfTwo(x));
128      // fall through
129      case DOWN:
130      case FLOOR:
131        return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
132
133      case UP:
134      case CEILING:
135        return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
136
137      case HALF_DOWN:
138      case HALF_UP:
139      case HALF_EVEN:
140        // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
141        int leadingZeros = Long.numberOfLeadingZeros(x);
142        long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
143        // floor(2^(logFloor + 0.5))
144        int logFloor = (Long.SIZE - 1) - leadingZeros;
145        return logFloor + lessThanBranchFree(cmp, x);
146    }
147    throw new AssertionError("impossible");
148  }
149
150  /** The biggest half power of two that fits into an unsigned long */
151  @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
152
153  /**
154   * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
155   *
156   * @throws IllegalArgumentException if {@code x <= 0}
157   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
158   *     is not a power of ten
159   */
160  @GwtIncompatible // TODO
161  @SuppressWarnings("fallthrough")
162  // TODO(kevinb): remove after this warning is disabled globally
163  public static int log10(long x, RoundingMode mode) {
164    checkPositive("x", x);
165    int logFloor = log10Floor(x);
166    long floorPow = powersOf10[logFloor];
167    switch (mode) {
168      case UNNECESSARY:
169        checkRoundingUnnecessary(x == floorPow);
170      // fall through
171      case FLOOR:
172      case DOWN:
173        return logFloor;
174      case CEILING:
175      case UP:
176        return logFloor + lessThanBranchFree(floorPow, x);
177      case HALF_DOWN:
178      case HALF_UP:
179      case HALF_EVEN:
180        // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
181        return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x);
182    }
183    throw new AssertionError();
184  }
185
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  @GwtIncompatible // TODO
212  @VisibleForTesting
213  static final long[] powersOf10 = {
214    1L,
215    10L,
216    100L,
217    1000L,
218    10000L,
219    100000L,
220    1000000L,
221    10000000L,
222    100000000L,
223    1000000000L,
224    10000000000L,
225    100000000000L,
226    1000000000000L,
227    10000000000000L,
228    100000000000000L,
229    1000000000000000L,
230    10000000000000000L,
231    100000000000000000L,
232    1000000000000000000L
233  };
234
235  // halfPowersOf10[i] = largest long less than 10^(i + 0.5)
236  @GwtIncompatible // TODO
237  @VisibleForTesting
238  static final long[] halfPowersOf10 = {
239    3L,
240    31L,
241    316L,
242    3162L,
243    31622L,
244    316227L,
245    3162277L,
246    31622776L,
247    316227766L,
248    3162277660L,
249    31622776601L,
250    316227766016L,
251    3162277660168L,
252    31622776601683L,
253    316227766016837L,
254    3162277660168379L,
255    31622776601683793L,
256    316227766016837933L,
257    3162277660168379331L
258  };
259
260  /**
261   * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
262   * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
263   * time.
264   *
265   * @throws IllegalArgumentException if {@code k < 0}
266   */
267  @GwtIncompatible // TODO
268  public static long pow(long b, int k) {
269    checkNonNegative("exponent", k);
270    if (-2 <= b && b <= 2) {
271      switch ((int) b) {
272        case 0:
273          return (k == 0) ? 1 : 0;
274        case 1:
275          return 1;
276        case -1:
277          return ((k & 1) == 0) ? 1 : -1;
278        case 2:
279          return (k < Long.SIZE) ? 1L << k : 0;
280        case -2:
281          if (k < Long.SIZE) {
282            return ((k & 1) == 0) ? 1L << k : -(1L << k);
283          } else {
284            return 0;
285          }
286        default:
287          throw new AssertionError();
288      }
289    }
290    for (long accum = 1; ; k >>= 1) {
291      switch (k) {
292        case 0:
293          return accum;
294        case 1:
295          return accum * b;
296        default:
297          accum *= ((k & 1) == 0) ? 1 : b;
298          b *= b;
299      }
300    }
301  }
302
303  /**
304   * Returns the square root of {@code x}, rounded with the specified rounding mode.
305   *
306   * @throws IllegalArgumentException if {@code x < 0}
307   * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
308   *     sqrt(x)} is not an integer
309   */
310  @GwtIncompatible // TODO
311  public static long sqrt(long x, RoundingMode mode) {
312    checkNonNegative("x", x);
313    if (fitsInInt(x)) {
314      return IntMath.sqrt((int) x, mode);
315    }
316    /*
317     * Let k be the true value of floor(sqrt(x)), so that
318     *
319     *            k * k <= x          <  (k + 1) * (k + 1)
320     * (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1))
321     *          since casting to double is nondecreasing.
322     *          Note that the right-hand inequality is no longer strict.
323     * Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1))
324     *          since Math.sqrt is monotonic.
325     * (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1))
326     *          since casting to long is monotonic
327     * k <= (long) Math.sqrt(x) <= k + 1
328     *          since (long) Math.sqrt(k * k) == k, as checked exhaustively in
329     *          {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect}
330     */
331    long guess = (long) Math.sqrt((double) x);
332    // Note: guess is always <= FLOOR_SQRT_MAX_LONG.
333    long guessSquared = guess * guess;
334    // Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is
335    // faster here than using lessThanBranchFree.
336    switch (mode) {
337      case UNNECESSARY:
338        checkRoundingUnnecessary(guessSquared == x);
339        return guess;
340      case FLOOR:
341      case DOWN:
342        if (x < guessSquared) {
343          return guess - 1;
344        }
345        return guess;
346      case CEILING:
347      case UP:
348        if (x > guessSquared) {
349          return guess + 1;
350        }
351        return guess;
352      case HALF_DOWN:
353      case HALF_UP:
354      case HALF_EVEN:
355        long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0);
356        long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
357        /*
358         * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
359         * and halfSquare are integers, this is equivalent to testing whether or not x <=
360         * halfSquare. (We have to deal with overflow, though.)
361         *
362         * If we treat halfSquare as an unsigned long, we know that
363         *            sqrtFloor^2 <= x < (sqrtFloor + 1)^2
364         * halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1
365         * so |x - halfSquare| <= sqrtFloor.  Therefore, it's safe to treat x - halfSquare as a
366         * signed long, so lessThanBranchFree is safe for use.
367         */
368        return sqrtFloor + lessThanBranchFree(halfSquare, x);
369    }
370    throw new AssertionError();
371  }
372
373  /**
374   * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
375   * RoundingMode}. If the {@code RoundingMode} is {@link RoundingMode#DOWN}, then this method is
376   * equivalent to regular Java division, {@code p / q}; and if it is {@link RoundingMode#FLOOR},
377   * then this method is equivalent to {@link Math#floorDiv(long,long) Math.floorDiv}{@code (p, q)}.
378   *
379   * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
380   *     is not an integer multiple of {@code b}
381   */
382  @GwtIncompatible // TODO
383  @SuppressWarnings("fallthrough")
384  public static long divide(long p, long q, RoundingMode mode) {
385    checkNotNull(mode);
386    long div = p / q; // throws if q == 0
387    long rem = p - q * div; // equals p % q
388
389    if (rem == 0) {
390      return div;
391    }
392
393    /*
394     * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
395     * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
396     * p / q.
397     *
398     * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
399     */
400    int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
401    boolean increment;
402    switch (mode) {
403      case UNNECESSARY:
404        checkRoundingUnnecessary(rem == 0);
405      // fall through
406      case DOWN:
407        increment = false;
408        break;
409      case UP:
410        increment = true;
411        break;
412      case CEILING:
413        increment = signum > 0;
414        break;
415      case FLOOR:
416        increment = signum < 0;
417        break;
418      case HALF_EVEN:
419      case HALF_DOWN:
420      case HALF_UP:
421        long absRem = abs(rem);
422        long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
423        // subtracting two nonnegative longs can't overflow
424        // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
425        if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
426          increment = (mode == HALF_UP || (mode == HALF_EVEN && (div & 1) != 0));
427        } else {
428          increment = cmpRemToHalfDivisor > 0; // closer to the UP value
429        }
430        break;
431      default:
432        throw new AssertionError();
433    }
434    return increment ? div + signum : div;
435  }
436
437  /**
438   * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
439   * m}, which might be negative.
440   *
441   * <p>For example:
442   *
443   * {@snippet :
444   * mod(7, 4) == 3
445   * mod(-7, 4) == 1
446   * mod(-1, 4) == 3
447   * mod(-8, 4) == 0
448   * mod(8, 4) == 0
449   * }
450   *
451   * @throws ArithmeticException if {@code m <= 0}
452   * @see <a href="https://docs.oracle.com/javase/specs/jls/se27/html/jls-15.html#jls-15.17.3">
453   *     Remainder Operator</a>
454   */
455  @GwtIncompatible // TODO
456  public static int mod(long x, int m) {
457    // Cast is safe because the result is guaranteed in the range [0, m)
458    return (int) mod(x, (long) m);
459  }
460
461  /**
462   * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
463   * m}, which might be negative.
464   *
465   * <p>For example:
466   *
467   * {@snippet :
468   * mod(7, 4) == 3
469   * mod(-7, 4) == 1
470   * mod(-1, 4) == 3
471   * mod(-8, 4) == 0
472   * mod(8, 4) == 0
473   * }
474   *
475   * @throws ArithmeticException if {@code m <= 0}
476   * @see <a href="https://docs.oracle.com/javase/specs/jls/se27/html/jls-15.html#jls-15.17.3">
477   *     Remainder Operator</a>
478   */
479  @GwtIncompatible // TODO
480  public static long mod(long x, long m) {
481    if (m <= 0) {
482      throw new ArithmeticException("Modulus must be positive");
483    }
484    return Math.floorMod(x, m);
485  }
486
487  /**
488   * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if {@code a == 0 && b ==
489   * 0}.
490   *
491   * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
492   */
493  public static long gcd(long a, long b) {
494    /*
495     * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
496     * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't an
497     * int.
498     */
499    checkNonNegative("a", a);
500    checkNonNegative("b", b);
501    if (a == 0) {
502      // 0 % b == 0, so b divides a, but the converse doesn't hold.
503      // BigInteger.gcd is consistent with this decision.
504      return b;
505    } else if (b == 0) {
506      return a; // similar logic
507    }
508    /*
509     * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. This is
510     * >60% faster than the Euclidean algorithm in benchmarks.
511     */
512    int aTwos = Long.numberOfTrailingZeros(a);
513    a >>= aTwos; // divide out all 2s
514    int bTwos = Long.numberOfTrailingZeros(b);
515    b >>= bTwos; // divide out all 2s
516    while (a != b) { // both a, b are odd
517      // The key to the binary GCD algorithm is as follows:
518      // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
519      // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
520
521      // We bend over backwards to avoid branching, adapting a technique from
522      // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
523
524      long delta = a - b; // can't overflow, since a and b are nonnegative
525
526      long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
527      // equivalent to Math.min(delta, 0)
528
529      a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
530      // a is now nonnegative and even
531
532      b += minDeltaOrZero; // sets b to min(old a, b)
533      a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
534    }
535    return a << min(aTwos, bTwos);
536  }
537
538  /**
539   * Returns the sum of {@code a} and {@code b}, provided it does not overflow.
540   *
541   * <p><b>Note:</b> this method is now unnecessary and should be treated as deprecated; use {@link
542   * Math#addExact(long, long)} instead. Note that if both arguments are {@code int} values, writing
543   * {@code Math.addExact(a, b)} will call the {@link Math#addExact(int, int)} overload, not {@link
544   * Math#addExact(long, long)}. Also note that adding two {@code int} values can <b>never</b>
545   * overflow a {@code long}, so you can just write {@code (long) a + b}.
546   *
547   * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
548   */
549  @InlineMe(replacement = "Math.addExact(a, b)")
550  public static long checkedAdd(long a, long b) {
551    return Math.addExact(a, b);
552  }
553
554  /**
555   * Returns the difference of {@code a} and {@code b}, provided it does not overflow.
556   *
557   * <p><b>Note:</b> this method is now unnecessary and should be treated as deprecated; use {@link
558   * Math#subtractExact(long, long)} instead. Note that if both arguments are {@code int} values,
559   * writing {@code Math.subtractExact(a, b)} will call the {@link Math#subtractExact(int, int)}
560   * overload, not {@link Math#subtractExact(long, long)}. Also note that subtracting two {@code
561   * int} values can <b>never</b> overflow a {@code long}, so you can just write {@code (long) a -
562   * b}.
563   *
564   * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
565   */
566  @InlineMe(replacement = "Math.subtractExact(a, b)")
567  public static long checkedSubtract(long a, long b) {
568    return Math.subtractExact(a, b);
569  }
570
571  /**
572   * Returns the product of {@code a} and {@code b}, provided it does not overflow.
573   *
574   * <p><b>Note:</b> this method is now unnecessary and should be treated as deprecated; use {@link
575   * Math#multiplyExact(long, long)} instead. Note that if both arguments are {@code int} values,
576   * writing {@code Math.multiplyExact(a, b)} will call the {@link Math#multiplyExact(int, int)}
577   * overload, not {@link Math#multiplyExact(long, long)}. Also note that multiplying two {@code
578   * int} values can <b>never</b> overflow a {@code long}, so you can just write {@code (long) a *
579   * b}.
580   *
581   * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
582   */
583  @InlineMe(replacement = "Math.multiplyExact(a, b)")
584  public static long checkedMultiply(long a, long b) {
585    return Math.multiplyExact(a, b);
586  }
587
588  /**
589   * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
590   *
591   * <p><b>Java 25+ users:</b> Prefer {@link Math#powExact(long, int)}. Note that if {@code b} is an
592   * {@code int} value, writing {@code Math.powExact(b, k)} will call the {@link Math#powExact(int,
593   * int)} overload, not {@link Math#powExact(long, int)}.
594   *
595   * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed {@code
596   *     long} arithmetic
597   */
598  @GwtIncompatible // TODO
599  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
600  @SuppressWarnings("ShortCircuitBoolean")
601  public static long checkedPow(long b, int k) {
602    checkNonNegative("exponent", k);
603    if (b >= -2 & b <= 2) {
604      switch ((int) b) {
605        case 0:
606          return (k == 0) ? 1 : 0;
607        case 1:
608          return 1;
609        case -1:
610          return ((k & 1) == 0) ? 1 : -1;
611        case 2:
612          checkNoOverflow(k < Long.SIZE - 1, "checkedPow", b, k);
613          return 1L << k;
614        case -2:
615          checkNoOverflow(k < Long.SIZE, "checkedPow", b, k);
616          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
617        default:
618          throw new AssertionError();
619      }
620    }
621    long accum = 1;
622    while (true) {
623      switch (k) {
624        case 0:
625          return accum;
626        case 1:
627          return Math.multiplyExact(accum, b);
628        default:
629          if ((k & 1) != 0) {
630            accum = Math.multiplyExact(accum, b);
631          }
632          k >>= 1;
633          if (k > 0) {
634            checkNoOverflow(
635                -FLOOR_SQRT_MAX_LONG <= b && b <= FLOOR_SQRT_MAX_LONG, "checkedPow", b, k);
636            b *= b;
637          }
638      }
639    }
640  }
641
642  /**
643   * Returns the sum of {@code a} and {@code b} unless it would overflow or underflow in which case
644   * {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
645   *
646   * @since 20.0
647   */
648  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
649  @SuppressWarnings("ShortCircuitBoolean")
650  public static long saturatedAdd(long a, long b) {
651    long naiveSum = a + b;
652    if ((a ^ b) < 0 | (a ^ naiveSum) >= 0) {
653      // If a and b have different signs or a has the same sign as the result then there was no
654      // overflow, return.
655      return naiveSum;
656    }
657    // we did over/under flow, if the sign is negative we should return MAX otherwise MIN
658    return Long.MAX_VALUE + ((naiveSum >>> (Long.SIZE - 1)) ^ 1);
659  }
660
661  /**
662   * Returns the difference of {@code a} and {@code b} unless it would overflow or underflow in
663   * which case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
664   *
665   * @since 20.0
666   */
667  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
668  @SuppressWarnings("ShortCircuitBoolean")
669  public static long saturatedSubtract(long a, long b) {
670    long naiveDifference = a - b;
671    if ((a ^ b) >= 0 | (a ^ naiveDifference) >= 0) {
672      // If a and b have the same signs or a has the same sign as the result then there was no
673      // overflow, return.
674      return naiveDifference;
675    }
676    // we did over/under flow
677    return Long.MAX_VALUE + ((naiveDifference >>> (Long.SIZE - 1)) ^ 1);
678  }
679
680  /**
681   * Returns the product of {@code a} and {@code b} unless it would overflow or underflow in which
682   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
683   *
684   * @since 20.0
685   */
686  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
687  @SuppressWarnings("ShortCircuitBoolean")
688  public static long saturatedMultiply(long a, long b) {
689    // see checkedMultiply for explanation
690    int leadingZeros =
691        Long.numberOfLeadingZeros(a)
692            + Long.numberOfLeadingZeros(~a)
693            + Long.numberOfLeadingZeros(b)
694            + Long.numberOfLeadingZeros(~b);
695    if (leadingZeros > Long.SIZE + 1) {
696      return a * b;
697    }
698    // the return value if we will overflow (which we calculate by overflowing a long :) )
699    long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
700    if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
701      // overflow
702      return limit;
703    }
704    long result = a * b;
705    if (a == 0 || result / a == b) {
706      return result;
707    }
708    return limit;
709  }
710
711  /**
712   * Returns the {@code b} to the {@code k}th power, unless it would overflow or underflow in which
713   * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
714   *
715   * @since 20.0
716   */
717  // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
718  @SuppressWarnings("ShortCircuitBoolean")
719  public static long saturatedPow(long b, int k) {
720    checkNonNegative("exponent", k);
721    if (b >= -2 & b <= 2) {
722      switch ((int) b) {
723        case 0:
724          return (k == 0) ? 1 : 0;
725        case 1:
726          return 1;
727        case -1:
728          return ((k & 1) == 0) ? 1 : -1;
729        case 2:
730          if (k >= Long.SIZE - 1) {
731            return Long.MAX_VALUE;
732          }
733          return 1L << k;
734        case -2:
735          if (k >= Long.SIZE) {
736            return Long.MAX_VALUE + (k & 1);
737          }
738          return ((k & 1) == 0) ? (1L << k) : (-1L << k);
739        default:
740          throw new AssertionError();
741      }
742    }
743    long accum = 1;
744    // if b is negative and k is odd then the limit is MIN otherwise the limit is MAX
745    long limit = Long.MAX_VALUE + ((b >>> (Long.SIZE - 1)) & (k & 1));
746    while (true) {
747      switch (k) {
748        case 0:
749          return accum;
750        case 1:
751          return saturatedMultiply(accum, b);
752        default:
753          if ((k & 1) != 0) {
754            accum = saturatedMultiply(accum, b);
755          }
756          k >>= 1;
757          if (k > 0) {
758            if (-FLOOR_SQRT_MAX_LONG > b | b > FLOOR_SQRT_MAX_LONG) {
759              return limit;
760            }
761            b *= b;
762          }
763      }
764    }
765  }
766
767  @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
768
769  /**
770   * Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
771   * {@code n == 0}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
772   *
773   * @throws IllegalArgumentException if {@code n < 0}
774   */
775  @GwtIncompatible // TODO
776  public static long factorial(int n) {
777    checkNonNegative("n", n);
778    return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
779  }
780
781  static final long[] factorials = {
782    1L,
783    1L,
784    1L * 2,
785    1L * 2 * 3,
786    1L * 2 * 3 * 4,
787    1L * 2 * 3 * 4 * 5,
788    1L * 2 * 3 * 4 * 5 * 6,
789    1L * 2 * 3 * 4 * 5 * 6 * 7,
790    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
791    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
792    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
793    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
794    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
795    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
796    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
797    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
798    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
799    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
800    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
801    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
802    1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
803  };
804
805  /**
806   * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
807   * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
808   *
809   * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
810   */
811  public static long binomial(int n, int k) {
812    checkNonNegative("n", n);
813    checkNonNegative("k", k);
814    checkArgument(k <= n, "k (%s) > n (%s)", k, n);
815    if (k > (n >> 1)) {
816      k = n - k;
817    }
818    switch (k) {
819      case 0:
820        return 1;
821      case 1:
822        return n;
823      default:
824        if (n < factorials.length) {
825          return factorials[n] / (factorials[k] * factorials[n - k]);
826        } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
827          return Long.MAX_VALUE;
828        } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
829          // guaranteed not to overflow
830          long result = n--;
831          for (int i = 2; i <= k; n--, i++) {
832            result *= n;
833            result /= i;
834          }
835          return result;
836        } else {
837          int nBits = log2(n, RoundingMode.CEILING);
838
839          long result = 1;
840          long numerator = n--;
841          long denominator = 1;
842
843          int numeratorBits = nBits;
844          // This is an upper bound on log2(numerator, ceiling).
845
846          /*
847           * We want to do this in long math for speed, but want to avoid overflow. We adapt the
848           * technique previously used by BigIntegerMath: maintain separate numerator and
849           * denominator accumulators, multiplying the fraction into result when near overflow.
850           */
851          for (int i = 2; i <= k; i++, n--) {
852            if (numeratorBits + nBits < Long.SIZE - 1) {
853              // It's definitely safe to multiply into numerator and denominator.
854              numerator *= n;
855              denominator *= i;
856              numeratorBits += nBits;
857            } else {
858              // It might not be safe to multiply into numerator and denominator,
859              // so multiply (numerator / denominator) into result.
860              result = multiplyFraction(result, numerator, denominator);
861              numerator = n;
862              denominator = i;
863              numeratorBits = nBits;
864            }
865          }
866          return multiplyFraction(result, numerator, denominator);
867        }
868    }
869  }
870
871  /** Returns (x * numerator / denominator), which is assumed to come out to an integral value. */
872  static long multiplyFraction(long x, long numerator, long denominator) {
873    if (x == 1) {
874      return numerator / denominator;
875    }
876    long commonDivisor = gcd(x, denominator);
877    x /= commonDivisor;
878    denominator /= commonDivisor;
879    // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
880    // so denominator must be a divisor of numerator.
881    return x * (numerator / denominator);
882  }
883
884  /*
885   * binomial(biggestBinomials[k], k) fits in a long, but not binomial(biggestBinomials[k] + 1, k).
886   */
887  static final int[] biggestBinomials = {
888    Integer.MAX_VALUE,
889    Integer.MAX_VALUE,
890    Integer.MAX_VALUE,
891    3810779,
892    121977,
893    16175,
894    4337,
895    1733,
896    887,
897    534,
898    361,
899    265,
900    206,
901    169,
902    143,
903    125,
904    111,
905    101,
906    94,
907    88,
908    83,
909    79,
910    76,
911    74,
912    72,
913    70,
914    69,
915    68,
916    67,
917    67,
918    66,
919    66,
920    66,
921    66
922  };
923
924  /*
925   * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, but
926   * binomial(biggestSimpleBinomials[k] + 1, k) does.
927   */
928  @VisibleForTesting
929  static final int[] biggestSimpleBinomials = {
930    Integer.MAX_VALUE,
931    Integer.MAX_VALUE,
932    Integer.MAX_VALUE,
933    2642246,
934    86251,
935    11724,
936    3218,
937    1313,
938    684,
939    419,
940    287,
941    214,
942    169,
943    139,
944    119,
945    105,
946    95,
947    87,
948    81,
949    76,
950    73,
951    70,
952    68,
953    66,
954    64,
955    63,
956    62,
957    62,
958    61,
959    61,
960    61
961  };
962
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      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      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  @GwtIncompatible
1247  public static double roundToDouble(long x, RoundingMode mode) {
1248    // Logic adapted from ToDoubleRounder.
1249    double roundArbitrarily = (double) x;
1250    long roundArbitrarilyAsLong = (long) roundArbitrarily;
1251    int cmpXToRoundArbitrarily;
1252
1253    if (roundArbitrarilyAsLong == Long.MAX_VALUE) {
1254      /*
1255       * For most values, the conversion from roundArbitrarily to roundArbitrarilyAsLong is
1256       * lossless. In that case we can compare x to roundArbitrarily using Long.compare(x,
1257       * roundArbitrarilyAsLong). The exception is for values where the conversion to double rounds
1258       * up to give roundArbitrarily equal to 2^63, so the conversion back to long overflows and
1259       * roundArbitrarilyAsLong is Long.MAX_VALUE. (This is the only way this condition can occur as
1260       * otherwise the conversion back to long pads with zero bits.) In this case we know that
1261       * roundArbitrarily > x. (This is important when x == Long.MAX_VALUE ==
1262       * roundArbitrarilyAsLong.)
1263       */
1264      cmpXToRoundArbitrarily = -1;
1265    } else {
1266      cmpXToRoundArbitrarily = Long.compare(x, roundArbitrarilyAsLong);
1267    }
1268
1269    switch (mode) {
1270      case UNNECESSARY:
1271        checkRoundingUnnecessary(cmpXToRoundArbitrarily == 0);
1272        return roundArbitrarily;
1273      case FLOOR:
1274        return (cmpXToRoundArbitrarily >= 0) ? roundArbitrarily : nextDown(roundArbitrarily);
1275      case CEILING:
1276        return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : nextUp(roundArbitrarily);
1277      case DOWN:
1278        if (x >= 0) {
1279          return (cmpXToRoundArbitrarily >= 0) ? roundArbitrarily : nextDown(roundArbitrarily);
1280        } else {
1281          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : nextUp(roundArbitrarily);
1282        }
1283      case UP:
1284        if (x >= 0) {
1285          return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : nextUp(roundArbitrarily);
1286        } else {
1287          return (cmpXToRoundArbitrarily >= 0) ? roundArbitrarily : nextDown(roundArbitrarily);
1288        }
1289      case HALF_DOWN:
1290      case HALF_UP:
1291      case HALF_EVEN:
1292        {
1293          long roundFloor;
1294          double roundFloorAsDouble;
1295          long roundCeiling;
1296          double roundCeilingAsDouble;
1297
1298          if (cmpXToRoundArbitrarily >= 0) {
1299            roundFloorAsDouble = roundArbitrarily;
1300            roundFloor = roundArbitrarilyAsLong;
1301            roundCeilingAsDouble = nextUp(roundArbitrarily);
1302            roundCeiling = (long) ceil(roundCeilingAsDouble);
1303          } else {
1304            roundCeilingAsDouble = roundArbitrarily;
1305            roundCeiling = roundArbitrarilyAsLong;
1306            roundFloorAsDouble = nextDown(roundArbitrarily);
1307            roundFloor = (long) floor(roundFloorAsDouble);
1308          }
1309
1310          long deltaToFloor = x - roundFloor;
1311          long deltaToCeiling = roundCeiling - x;
1312
1313          if (roundCeiling == Long.MAX_VALUE) {
1314            // correct for Long.MAX_VALUE as discussed above: roundCeilingAsDouble must be 2^63, but
1315            // roundCeiling is 2^63-1.
1316            deltaToCeiling++;
1317          }
1318
1319          int diff = Long.compare(deltaToFloor, deltaToCeiling);
1320          if (diff < 0) { // closer to floor
1321            return roundFloorAsDouble;
1322          } else if (diff > 0) { // closer to ceiling
1323            return roundCeilingAsDouble;
1324          }
1325          // halfway between the representable values; do the half-whatever logic
1326          switch (mode) {
1327            case HALF_EVEN:
1328              return ((DoubleUtils.getSignificand(roundFloorAsDouble) & 1L) == 0)
1329                  ? roundFloorAsDouble
1330                  : roundCeilingAsDouble;
1331            case HALF_DOWN:
1332              return (x >= 0) ? roundFloorAsDouble : roundCeilingAsDouble;
1333            case HALF_UP:
1334              return (x >= 0) ? roundCeilingAsDouble : roundFloorAsDouble;
1335            default:
1336              throw new AssertionError("impossible");
1337          }
1338        }
1339    }
1340    throw new AssertionError("impossible");
1341  }
1342
1343  /**
1344   * Returns the closest representable {@code long} to the absolute value of {@code x}.
1345   *
1346   * <p>This is the same thing as the true absolute value of {@code x} except in the case when
1347   * {@code x} is {@link Long#MIN_VALUE}, in which case this returns {@link Long#MAX_VALUE}. (Note
1348   * that {@code Long.MAX_VALUE} is mathematically equal to {@code -Long.MIN_VALUE - 1}.)
1349   *
1350   * <p>There are three common APIs for determining the absolute value of a long, all of which
1351   * behave identically except when passed {@code Long.MIN_VALUE}. Those methods are:
1352   *
1353   * <ul>
1354   *   <li>{@link Math#abs(long)}, which returns {@code Long.MIN_VALUE} when passed {@code
1355   *       Long.MIN_VALUE}
1356   *   <li>{@link Math#absExact(long)}, which throws {@link ArithmeticException} when passed {@code
1357   *       Long.MIN_VALUE}
1358   *   <li>this method, {@code LongMath.saturatedAbs(long)}, which returns {@code Long.MAX_VALUE}
1359   *       when passed {@code Long.MIN_VALUE}
1360   * </ul>
1361   *
1362   * <p>Note that if your only goal is to turn a well-distributed {@code long} (such as a random
1363   * number) into a well-distributed nonnegative number, the most even distribution is achieved not
1364   * by this method or other absolute value methods, but by {@code x & Long.MAX_VALUE}.
1365   *
1366   * @since 33.5.0
1367   */
1368  public static long saturatedAbs(long x) {
1369    return (x == Long.MIN_VALUE) ? Long.MAX_VALUE : Math.abs(x);
1370  }
1371
1372  private LongMath() {}
1373}