• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (C) 2011 The Guava Authors
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
5  * in compliance with the License. You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software distributed under the License
10  * is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
11  * or implied. See the License for the specific language governing permissions and limitations under
12  * the License.
13  */
14 
15 package com.google.common.math;
16 
17 import static com.google.common.base.Preconditions.checkArgument;
18 import static com.google.common.base.Preconditions.checkNotNull;
19 import static com.google.common.math.MathPreconditions.checkNoOverflow;
20 import static com.google.common.math.MathPreconditions.checkNonNegative;
21 import static com.google.common.math.MathPreconditions.checkPositive;
22 import static com.google.common.math.MathPreconditions.checkRoundingUnnecessary;
23 import static java.lang.Math.abs;
24 import static java.lang.Math.min;
25 import static java.math.RoundingMode.HALF_EVEN;
26 import static java.math.RoundingMode.HALF_UP;
27 
28 import com.google.common.annotations.GwtCompatible;
29 import com.google.common.annotations.GwtIncompatible;
30 import com.google.common.annotations.VisibleForTesting;
31 import com.google.common.primitives.UnsignedLongs;
32 import java.math.BigInteger;
33 import java.math.RoundingMode;
34 
35 /**
36  * A class for arithmetic on values of type {@code long}. Where possible, methods are defined and
37  * named analogously to their {@code BigInteger} counterparts.
38  *
39  * <p>The implementations of many methods in this class are based on material from Henry S. Warren,
40  * Jr.'s <i>Hacker's Delight</i>, (Addison Wesley, 2002).
41  *
42  * <p>Similar functionality for {@code int} and for {@link BigInteger} can be found in {@link
43  * IntMath} and {@link BigIntegerMath} respectively. For other common operations on {@code long}
44  * values, see {@link com.google.common.primitives.Longs}.
45  *
46  * @author Louis Wasserman
47  * @since 11.0
48  */
49 @GwtCompatible(emulated = true)
50 @ElementTypesAreNonnullByDefault
51 public final class LongMath {
52   @VisibleForTesting static final long MAX_SIGNED_POWER_OF_TWO = 1L << (Long.SIZE - 2);
53 
54   /**
55    * Returns the smallest power of two greater than or equal to {@code x}. This is equivalent to
56    * {@code checkedPow(2, log2(x, CEILING))}.
57    *
58    * @throws IllegalArgumentException if {@code x <= 0}
59    * @throws ArithmeticException of the next-higher power of two is not representable as a {@code
60    *     long}, i.e. when {@code x > 2^62}
61    * @since 20.0
62    */
ceilingPowerOfTwo(long x)63   public static long ceilingPowerOfTwo(long x) {
64     checkPositive("x", x);
65     if (x > MAX_SIGNED_POWER_OF_TWO) {
66       throw new ArithmeticException("ceilingPowerOfTwo(" + x + ") is not representable as a long");
67     }
68     return 1L << -Long.numberOfLeadingZeros(x - 1);
69   }
70 
71   /**
72    * Returns the largest power of two less than or equal to {@code x}. This is equivalent to {@code
73    * checkedPow(2, log2(x, FLOOR))}.
74    *
75    * @throws IllegalArgumentException if {@code x <= 0}
76    * @since 20.0
77    */
floorPowerOfTwo(long x)78   public static long floorPowerOfTwo(long x) {
79     checkPositive("x", x);
80 
81     // Long.highestOneBit was buggy on GWT.  We've fixed it, but I'm not certain when the fix will
82     // be released.
83     return 1L << ((Long.SIZE - 1) - Long.numberOfLeadingZeros(x));
84   }
85 
86   /**
87    * Returns {@code true} if {@code x} represents a power of two.
88    *
89    * <p>This differs from {@code Long.bitCount(x) == 1}, because {@code
90    * Long.bitCount(Long.MIN_VALUE) == 1}, but {@link Long#MIN_VALUE} is not a power of two.
91    */
92   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
93   @SuppressWarnings("ShortCircuitBoolean")
isPowerOfTwo(long x)94   public static boolean isPowerOfTwo(long x) {
95     return x > 0 & (x & (x - 1)) == 0;
96   }
97 
98   /**
99    * Returns 1 if {@code x < y} as unsigned longs, and 0 otherwise. Assumes that x - y fits into a
100    * signed long. The implementation is branch-free, and benchmarks suggest it is measurably faster
101    * than the straightforward ternary expression.
102    */
103   @VisibleForTesting
lessThanBranchFree(long x, long y)104   static int lessThanBranchFree(long x, long y) {
105     // Returns the sign bit of x - y.
106     return (int) (~~(x - y) >>> (Long.SIZE - 1));
107   }
108 
109   /**
110    * Returns the base-2 logarithm of {@code x}, rounded according to the specified rounding mode.
111    *
112    * @throws IllegalArgumentException if {@code x <= 0}
113    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
114    *     is not a power of two
115    */
116   @SuppressWarnings("fallthrough")
117   // TODO(kevinb): remove after this warning is disabled globally
log2(long x, RoundingMode mode)118   public static int log2(long x, RoundingMode mode) {
119     checkPositive("x", x);
120     switch (mode) {
121       case UNNECESSARY:
122         checkRoundingUnnecessary(isPowerOfTwo(x));
123         // fall through
124       case DOWN:
125       case FLOOR:
126         return (Long.SIZE - 1) - Long.numberOfLeadingZeros(x);
127 
128       case UP:
129       case CEILING:
130         return Long.SIZE - Long.numberOfLeadingZeros(x - 1);
131 
132       case HALF_DOWN:
133       case HALF_UP:
134       case HALF_EVEN:
135         // Since sqrt(2) is irrational, log2(x) - logFloor cannot be exactly 0.5
136         int leadingZeros = Long.numberOfLeadingZeros(x);
137         long cmp = MAX_POWER_OF_SQRT2_UNSIGNED >>> leadingZeros;
138         // floor(2^(logFloor + 0.5))
139         int logFloor = (Long.SIZE - 1) - leadingZeros;
140         return logFloor + lessThanBranchFree(cmp, x);
141     }
142     throw new AssertionError("impossible");
143   }
144 
145   /** The biggest half power of two that fits into an unsigned long */
146   @VisibleForTesting static final long MAX_POWER_OF_SQRT2_UNSIGNED = 0xB504F333F9DE6484L;
147 
148   /**
149    * Returns the base-10 logarithm of {@code x}, rounded according to the specified rounding mode.
150    *
151    * @throws IllegalArgumentException if {@code x <= 0}
152    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
153    *     is not a power of ten
154    */
155   @GwtIncompatible // TODO
156   @SuppressWarnings("fallthrough")
157   // TODO(kevinb): remove after this warning is disabled globally
log10(long x, RoundingMode mode)158   public static int log10(long x, RoundingMode mode) {
159     checkPositive("x", x);
160     int logFloor = log10Floor(x);
161     long floorPow = powersOf10[logFloor];
162     switch (mode) {
163       case UNNECESSARY:
164         checkRoundingUnnecessary(x == floorPow);
165         // fall through
166       case FLOOR:
167       case DOWN:
168         return logFloor;
169       case CEILING:
170       case UP:
171         return logFloor + lessThanBranchFree(floorPow, x);
172       case HALF_DOWN:
173       case HALF_UP:
174       case HALF_EVEN:
175         // sqrt(10) is irrational, so log10(x)-logFloor is never exactly 0.5
176         return logFloor + lessThanBranchFree(halfPowersOf10[logFloor], x);
177     }
178     throw new AssertionError();
179   }
180 
181   @GwtIncompatible // TODO
log10Floor(long x)182   static int log10Floor(long x) {
183     /*
184      * Based on Hacker's Delight Fig. 11-5, the two-table-lookup, branch-free implementation.
185      *
186      * The key idea is that based on the number of leading zeros (equivalently, floor(log2(x))), we
187      * can narrow the possible floor(log10(x)) values to two. For example, if floor(log2(x)) is 6,
188      * then 64 <= x < 128, so floor(log10(x)) is either 1 or 2.
189      */
190     int y = maxLog10ForLeadingZeros[Long.numberOfLeadingZeros(x)];
191     /*
192      * y is the higher of the two possible values of floor(log10(x)). If x < 10^y, then we want the
193      * lower of the two possible values, or y - 1, otherwise, we want y.
194      */
195     return y - lessThanBranchFree(x, powersOf10[y]);
196   }
197 
198   // maxLog10ForLeadingZeros[i] == floor(log10(2^(Long.SIZE - i)))
199   @VisibleForTesting
200   static final byte[] maxLog10ForLeadingZeros = {
201     19, 18, 18, 18, 18, 17, 17, 17, 16, 16, 16, 15, 15, 15, 15, 14, 14, 14, 13, 13, 13, 12, 12, 12,
202     12, 11, 11, 11, 10, 10, 10, 9, 9, 9, 9, 8, 8, 8, 7, 7, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3,
203     3, 2, 2, 2, 1, 1, 1, 0, 0, 0
204   };
205 
206   @GwtIncompatible // TODO
207   @VisibleForTesting
208   static final long[] powersOf10 = {
209     1L,
210     10L,
211     100L,
212     1000L,
213     10000L,
214     100000L,
215     1000000L,
216     10000000L,
217     100000000L,
218     1000000000L,
219     10000000000L,
220     100000000000L,
221     1000000000000L,
222     10000000000000L,
223     100000000000000L,
224     1000000000000000L,
225     10000000000000000L,
226     100000000000000000L,
227     1000000000000000000L
228   };
229 
230   // halfPowersOf10[i] = largest long less than 10^(i + 0.5)
231   @GwtIncompatible // TODO
232   @VisibleForTesting
233   static final long[] halfPowersOf10 = {
234     3L,
235     31L,
236     316L,
237     3162L,
238     31622L,
239     316227L,
240     3162277L,
241     31622776L,
242     316227766L,
243     3162277660L,
244     31622776601L,
245     316227766016L,
246     3162277660168L,
247     31622776601683L,
248     316227766016837L,
249     3162277660168379L,
250     31622776601683793L,
251     316227766016837933L,
252     3162277660168379331L
253   };
254 
255   /**
256    * Returns {@code b} to the {@code k}th power. Even if the result overflows, it will be equal to
257    * {@code BigInteger.valueOf(b).pow(k).longValue()}. This implementation runs in {@code O(log k)}
258    * time.
259    *
260    * @throws IllegalArgumentException if {@code k < 0}
261    */
262   @GwtIncompatible // TODO
pow(long b, int k)263   public static long pow(long b, int k) {
264     checkNonNegative("exponent", k);
265     if (-2 <= b && b <= 2) {
266       switch ((int) b) {
267         case 0:
268           return (k == 0) ? 1 : 0;
269         case 1:
270           return 1;
271         case (-1):
272           return ((k & 1) == 0) ? 1 : -1;
273         case 2:
274           return (k < Long.SIZE) ? 1L << k : 0;
275         case (-2):
276           if (k < Long.SIZE) {
277             return ((k & 1) == 0) ? 1L << k : -(1L << k);
278           } else {
279             return 0;
280           }
281         default:
282           throw new AssertionError();
283       }
284     }
285     for (long accum = 1; ; k >>= 1) {
286       switch (k) {
287         case 0:
288           return accum;
289         case 1:
290           return accum * b;
291         default:
292           accum *= ((k & 1) == 0) ? 1 : b;
293           b *= b;
294       }
295     }
296   }
297 
298   /**
299    * Returns the square root of {@code x}, rounded with the specified rounding mode.
300    *
301    * @throws IllegalArgumentException if {@code x < 0}
302    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code
303    *     sqrt(x)} is not an integer
304    */
305   @GwtIncompatible // TODO
sqrt(long x, RoundingMode mode)306   public static long sqrt(long x, RoundingMode mode) {
307     checkNonNegative("x", x);
308     if (fitsInInt(x)) {
309       return IntMath.sqrt((int) x, mode);
310     }
311     /*
312      * Let k be the true value of floor(sqrt(x)), so that
313      *
314      *            k * k <= x          <  (k + 1) * (k + 1)
315      * (double) (k * k) <= (double) x <= (double) ((k + 1) * (k + 1))
316      *          since casting to double is nondecreasing.
317      *          Note that the right-hand inequality is no longer strict.
318      * Math.sqrt(k * k) <= Math.sqrt(x) <= Math.sqrt((k + 1) * (k + 1))
319      *          since Math.sqrt is monotonic.
320      * (long) Math.sqrt(k * k) <= (long) Math.sqrt(x) <= (long) Math.sqrt((k + 1) * (k + 1))
321      *          since casting to long is monotonic
322      * k <= (long) Math.sqrt(x) <= k + 1
323      *          since (long) Math.sqrt(k * k) == k, as checked exhaustively in
324      *          {@link LongMathTest#testSqrtOfPerfectSquareAsDoubleIsPerfect}
325      */
326     long guess = (long) Math.sqrt((double) x);
327     // Note: guess is always <= FLOOR_SQRT_MAX_LONG.
328     long guessSquared = guess * guess;
329     // Note (2013-2-26): benchmarks indicate that, inscrutably enough, using if statements is
330     // faster here than using lessThanBranchFree.
331     switch (mode) {
332       case UNNECESSARY:
333         checkRoundingUnnecessary(guessSquared == x);
334         return guess;
335       case FLOOR:
336       case DOWN:
337         if (x < guessSquared) {
338           return guess - 1;
339         }
340         return guess;
341       case CEILING:
342       case UP:
343         if (x > guessSquared) {
344           return guess + 1;
345         }
346         return guess;
347       case HALF_DOWN:
348       case HALF_UP:
349       case HALF_EVEN:
350         long sqrtFloor = guess - ((x < guessSquared) ? 1 : 0);
351         long halfSquare = sqrtFloor * sqrtFloor + sqrtFloor;
352         /*
353          * We wish to test whether or not x <= (sqrtFloor + 0.5)^2 = halfSquare + 0.25. Since both x
354          * and halfSquare are integers, this is equivalent to testing whether or not x <=
355          * halfSquare. (We have to deal with overflow, though.)
356          *
357          * If we treat halfSquare as an unsigned long, we know that
358          *            sqrtFloor^2 <= x < (sqrtFloor + 1)^2
359          * halfSquare - sqrtFloor <= x < halfSquare + sqrtFloor + 1
360          * so |x - halfSquare| <= sqrtFloor.  Therefore, it's safe to treat x - halfSquare as a
361          * signed long, so lessThanBranchFree is safe for use.
362          */
363         return sqrtFloor + lessThanBranchFree(halfSquare, x);
364     }
365     throw new AssertionError();
366   }
367 
368   /**
369    * Returns the result of dividing {@code p} by {@code q}, rounding using the specified {@code
370    * RoundingMode}.
371    *
372    * @throws ArithmeticException if {@code q == 0}, or if {@code mode == UNNECESSARY} and {@code a}
373    *     is not an integer multiple of {@code b}
374    */
375   @GwtIncompatible // TODO
376   @SuppressWarnings("fallthrough")
divide(long p, long q, RoundingMode mode)377   public static long divide(long p, long q, RoundingMode mode) {
378     checkNotNull(mode);
379     long div = p / q; // throws if q == 0
380     long rem = p - q * div; // equals p % q
381 
382     if (rem == 0) {
383       return div;
384     }
385 
386     /*
387      * Normal Java division rounds towards 0, consistently with RoundingMode.DOWN. We just have to
388      * deal with the cases where rounding towards 0 is wrong, which typically depends on the sign of
389      * p / q.
390      *
391      * signum is 1 if p and q are both nonnegative or both negative, and -1 otherwise.
392      */
393     int signum = 1 | (int) ((p ^ q) >> (Long.SIZE - 1));
394     boolean increment;
395     switch (mode) {
396       case UNNECESSARY:
397         checkRoundingUnnecessary(rem == 0);
398         // fall through
399       case DOWN:
400         increment = false;
401         break;
402       case UP:
403         increment = true;
404         break;
405       case CEILING:
406         increment = signum > 0;
407         break;
408       case FLOOR:
409         increment = signum < 0;
410         break;
411       case HALF_EVEN:
412       case HALF_DOWN:
413       case HALF_UP:
414         long absRem = abs(rem);
415         long cmpRemToHalfDivisor = absRem - (abs(q) - absRem);
416         // subtracting two nonnegative longs can't overflow
417         // cmpRemToHalfDivisor has the same sign as compare(abs(rem), abs(q) / 2).
418         if (cmpRemToHalfDivisor == 0) { // exactly on the half mark
419           increment = (mode == HALF_UP || (mode == HALF_EVEN && (div & 1) != 0));
420         } else {
421           increment = cmpRemToHalfDivisor > 0; // closer to the UP value
422         }
423         break;
424       default:
425         throw new AssertionError();
426     }
427     return increment ? div + signum : div;
428   }
429 
430   /**
431    * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
432    * m}, which might be negative.
433    *
434    * <p>For example:
435    *
436    * <pre>{@code
437    * mod(7, 4) == 3
438    * mod(-7, 4) == 1
439    * mod(-1, 4) == 3
440    * mod(-8, 4) == 0
441    * mod(8, 4) == 0
442    * }</pre>
443    *
444    * @throws ArithmeticException if {@code m <= 0}
445    * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
446    *     Remainder Operator</a>
447    */
448   @GwtIncompatible // TODO
mod(long x, int m)449   public static int mod(long x, int m) {
450     // Cast is safe because the result is guaranteed in the range [0, m)
451     return (int) mod(x, (long) m);
452   }
453 
454   /**
455    * Returns {@code x mod m}, a non-negative value less than {@code m}. This differs from {@code x %
456    * m}, which might be negative.
457    *
458    * <p>For example:
459    *
460    * <pre>{@code
461    * mod(7, 4) == 3
462    * mod(-7, 4) == 1
463    * mod(-1, 4) == 3
464    * mod(-8, 4) == 0
465    * mod(8, 4) == 0
466    * }</pre>
467    *
468    * @throws ArithmeticException if {@code m <= 0}
469    * @see <a href="http://docs.oracle.com/javase/specs/jls/se7/html/jls-15.html#jls-15.17.3">
470    *     Remainder Operator</a>
471    */
472   @GwtIncompatible // TODO
mod(long x, long m)473   public static long mod(long x, long m) {
474     if (m <= 0) {
475       throw new ArithmeticException("Modulus must be positive");
476     }
477     long result = x % m;
478     return (result >= 0) ? result : result + m;
479   }
480 
481   /**
482    * Returns the greatest common divisor of {@code a, b}. Returns {@code 0} if {@code a == 0 && b ==
483    * 0}.
484    *
485    * @throws IllegalArgumentException if {@code a < 0} or {@code b < 0}
486    */
gcd(long a, long b)487   public static long gcd(long a, long b) {
488     /*
489      * The reason we require both arguments to be >= 0 is because otherwise, what do you return on
490      * gcd(0, Long.MIN_VALUE)? BigInteger.gcd would return positive 2^63, but positive 2^63 isn't an
491      * int.
492      */
493     checkNonNegative("a", a);
494     checkNonNegative("b", b);
495     if (a == 0) {
496       // 0 % b == 0, so b divides a, but the converse doesn't hold.
497       // BigInteger.gcd is consistent with this decision.
498       return b;
499     } else if (b == 0) {
500       return a; // similar logic
501     }
502     /*
503      * Uses the binary GCD algorithm; see http://en.wikipedia.org/wiki/Binary_GCD_algorithm. This is
504      * >60% faster than the Euclidean algorithm in benchmarks.
505      */
506     int aTwos = Long.numberOfTrailingZeros(a);
507     a >>= aTwos; // divide out all 2s
508     int bTwos = Long.numberOfTrailingZeros(b);
509     b >>= bTwos; // divide out all 2s
510     while (a != b) { // both a, b are odd
511       // The key to the binary GCD algorithm is as follows:
512       // Both a and b are odd. Assume a > b; then gcd(a - b, b) = gcd(a, b).
513       // But in gcd(a - b, b), a - b is even and b is odd, so we can divide out powers of two.
514 
515       // We bend over backwards to avoid branching, adapting a technique from
516       // http://graphics.stanford.edu/~seander/bithacks.html#IntegerMinOrMax
517 
518       long delta = a - b; // can't overflow, since a and b are nonnegative
519 
520       long minDeltaOrZero = delta & (delta >> (Long.SIZE - 1));
521       // equivalent to Math.min(delta, 0)
522 
523       a = delta - minDeltaOrZero - minDeltaOrZero; // sets a to Math.abs(a - b)
524       // a is now nonnegative and even
525 
526       b += minDeltaOrZero; // sets b to min(old a, b)
527       a >>= Long.numberOfTrailingZeros(a); // divide out all 2s, since 2 doesn't divide b
528     }
529     return a << min(aTwos, bTwos);
530   }
531 
532   /**
533    * Returns the sum of {@code a} and {@code b}, provided it does not overflow.
534    *
535    * @throws ArithmeticException if {@code a + b} overflows in signed {@code long} arithmetic
536    */
537   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
538   @SuppressWarnings("ShortCircuitBoolean")
checkedAdd(long a, long b)539   public static long checkedAdd(long a, long b) {
540     long result = a + b;
541     checkNoOverflow((a ^ b) < 0 | (a ^ result) >= 0, "checkedAdd", a, b);
542     return result;
543   }
544 
545   /**
546    * Returns the difference of {@code a} and {@code b}, provided it does not overflow.
547    *
548    * @throws ArithmeticException if {@code a - b} overflows in signed {@code long} arithmetic
549    */
550   @GwtIncompatible // TODO
551   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
552   @SuppressWarnings("ShortCircuitBoolean")
checkedSubtract(long a, long b)553   public static long checkedSubtract(long a, long b) {
554     long result = a - b;
555     checkNoOverflow((a ^ b) >= 0 | (a ^ result) >= 0, "checkedSubtract", a, b);
556     return result;
557   }
558 
559   /**
560    * Returns the product of {@code a} and {@code b}, provided it does not overflow.
561    *
562    * @throws ArithmeticException if {@code a * b} overflows in signed {@code long} arithmetic
563    */
564   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
565   @SuppressWarnings("ShortCircuitBoolean")
checkedMultiply(long a, long b)566   public static long checkedMultiply(long a, long b) {
567     // Hacker's Delight, Section 2-12
568     int leadingZeros =
569         Long.numberOfLeadingZeros(a)
570             + Long.numberOfLeadingZeros(~a)
571             + Long.numberOfLeadingZeros(b)
572             + Long.numberOfLeadingZeros(~b);
573     /*
574      * If leadingZeros > Long.SIZE + 1 it's definitely fine, if it's < Long.SIZE it's definitely
575      * bad. We do the leadingZeros check to avoid the division below if at all possible.
576      *
577      * Otherwise, if b == Long.MIN_VALUE, then the only allowed values of a are 0 and 1. We take
578      * care of all a < 0 with their own check, because in particular, the case a == -1 will
579      * incorrectly pass the division check below.
580      *
581      * In all other cases, we check that either a is 0 or the result is consistent with division.
582      */
583     if (leadingZeros > Long.SIZE + 1) {
584       return a * b;
585     }
586     checkNoOverflow(leadingZeros >= Long.SIZE, "checkedMultiply", a, b);
587     checkNoOverflow(a >= 0 | b != Long.MIN_VALUE, "checkedMultiply", a, b);
588     long result = a * b;
589     checkNoOverflow(a == 0 || result / a == b, "checkedMultiply", a, b);
590     return result;
591   }
592 
593   /**
594    * Returns the {@code b} to the {@code k}th power, provided it does not overflow.
595    *
596    * @throws ArithmeticException if {@code b} to the {@code k}th power overflows in signed {@code
597    *     long} arithmetic
598    */
599   @GwtIncompatible // TODO
600   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
601   @SuppressWarnings("ShortCircuitBoolean")
checkedPow(long b, int k)602   public static long checkedPow(long b, int k) {
603     checkNonNegative("exponent", k);
604     if (b >= -2 & b <= 2) {
605       switch ((int) b) {
606         case 0:
607           return (k == 0) ? 1 : 0;
608         case 1:
609           return 1;
610         case (-1):
611           return ((k & 1) == 0) ? 1 : -1;
612         case 2:
613           checkNoOverflow(k < Long.SIZE - 1, "checkedPow", b, k);
614           return 1L << k;
615         case (-2):
616           checkNoOverflow(k < Long.SIZE, "checkedPow", b, k);
617           return ((k & 1) == 0) ? (1L << k) : (-1L << k);
618         default:
619           throw new AssertionError();
620       }
621     }
622     long accum = 1;
623     while (true) {
624       switch (k) {
625         case 0:
626           return accum;
627         case 1:
628           return checkedMultiply(accum, b);
629         default:
630           if ((k & 1) != 0) {
631             accum = checkedMultiply(accum, b);
632           }
633           k >>= 1;
634           if (k > 0) {
635             checkNoOverflow(
636                 -FLOOR_SQRT_MAX_LONG <= b && b <= FLOOR_SQRT_MAX_LONG, "checkedPow", b, k);
637             b *= b;
638           }
639       }
640     }
641   }
642 
643   /**
644    * Returns the sum of {@code a} and {@code b} unless it would overflow or underflow in which case
645    * {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
646    *
647    * @since 20.0
648    */
649   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
650   @SuppressWarnings("ShortCircuitBoolean")
651   public static long saturatedAdd(long a, long b) {
652     long naiveSum = a + b;
653     if ((a ^ b) < 0 | (a ^ naiveSum) >= 0) {
654       // If a and b have different signs or a has the same sign as the result then there was no
655       // overflow, return.
656       return naiveSum;
657     }
658     // we did over/under flow, if the sign is negative we should return MAX otherwise MIN
659     return Long.MAX_VALUE + ((naiveSum >>> (Long.SIZE - 1)) ^ 1);
660   }
661 
662   /**
663    * Returns the difference of {@code a} and {@code b} unless it would overflow or underflow in
664    * which case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
665    *
666    * @since 20.0
667    */
668   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
669   @SuppressWarnings("ShortCircuitBoolean")
670   public static long saturatedSubtract(long a, long b) {
671     long naiveDifference = a - b;
672     if ((a ^ b) >= 0 | (a ^ naiveDifference) >= 0) {
673       // If a and b have the same signs or a has the same sign as the result then there was no
674       // overflow, return.
675       return naiveDifference;
676     }
677     // we did over/under flow
678     return Long.MAX_VALUE + ((naiveDifference >>> (Long.SIZE - 1)) ^ 1);
679   }
680 
681   /**
682    * Returns the product of {@code a} and {@code b} unless it would overflow or underflow in which
683    * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
684    *
685    * @since 20.0
686    */
687   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
688   @SuppressWarnings("ShortCircuitBoolean")
saturatedMultiply(long a, long b)689   public static long saturatedMultiply(long a, long b) {
690     // see checkedMultiply for explanation
691     int leadingZeros =
692         Long.numberOfLeadingZeros(a)
693             + Long.numberOfLeadingZeros(~a)
694             + Long.numberOfLeadingZeros(b)
695             + Long.numberOfLeadingZeros(~b);
696     if (leadingZeros > Long.SIZE + 1) {
697       return a * b;
698     }
699     // the return value if we will overflow (which we calculate by overflowing a long :) )
700     long limit = Long.MAX_VALUE + ((a ^ b) >>> (Long.SIZE - 1));
701     if (leadingZeros < Long.SIZE | (a < 0 & b == Long.MIN_VALUE)) {
702       // overflow
703       return limit;
704     }
705     long result = a * b;
706     if (a == 0 || result / a == b) {
707       return result;
708     }
709     return limit;
710   }
711 
712   /**
713    * Returns the {@code b} to the {@code k}th power, unless it would overflow or underflow in which
714    * case {@code Long.MAX_VALUE} or {@code Long.MIN_VALUE} is returned, respectively.
715    *
716    * @since 20.0
717    */
718   // Whenever both tests are cheap and functional, it's faster to use &, | instead of &&, ||
719   @SuppressWarnings("ShortCircuitBoolean")
saturatedPow(long b, int k)720   public static long saturatedPow(long b, int k) {
721     checkNonNegative("exponent", k);
722     if (b >= -2 & b <= 2) {
723       switch ((int) b) {
724         case 0:
725           return (k == 0) ? 1 : 0;
726         case 1:
727           return 1;
728         case (-1):
729           return ((k & 1) == 0) ? 1 : -1;
730         case 2:
731           if (k >= Long.SIZE - 1) {
732             return Long.MAX_VALUE;
733           }
734           return 1L << k;
735         case (-2):
736           if (k >= Long.SIZE) {
737             return Long.MAX_VALUE + (k & 1);
738           }
739           return ((k & 1) == 0) ? (1L << k) : (-1L << k);
740         default:
741           throw new AssertionError();
742       }
743     }
744     long accum = 1;
745     // if b is negative and k is odd then the limit is MIN otherwise the limit is MAX
746     long limit = Long.MAX_VALUE + ((b >>> (Long.SIZE - 1)) & (k & 1));
747     while (true) {
748       switch (k) {
749         case 0:
750           return accum;
751         case 1:
752           return saturatedMultiply(accum, b);
753         default:
754           if ((k & 1) != 0) {
755             accum = saturatedMultiply(accum, b);
756           }
757           k >>= 1;
758           if (k > 0) {
759             if (-FLOOR_SQRT_MAX_LONG > b | b > FLOOR_SQRT_MAX_LONG) {
760               return limit;
761             }
762             b *= b;
763           }
764       }
765     }
766   }
767 
768   @VisibleForTesting static final long FLOOR_SQRT_MAX_LONG = 3037000499L;
769 
770   /**
771    * Returns {@code n!}, that is, the product of the first {@code n} positive integers, {@code 1} if
772    * {@code n == 0}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
773    *
774    * @throws IllegalArgumentException if {@code n < 0}
775    */
776   @GwtIncompatible // TODO
factorial(int n)777   public static long factorial(int n) {
778     checkNonNegative("n", n);
779     return (n < factorials.length) ? factorials[n] : Long.MAX_VALUE;
780   }
781 
782   static final long[] factorials = {
783     1L,
784     1L,
785     1L * 2,
786     1L * 2 * 3,
787     1L * 2 * 3 * 4,
788     1L * 2 * 3 * 4 * 5,
789     1L * 2 * 3 * 4 * 5 * 6,
790     1L * 2 * 3 * 4 * 5 * 6 * 7,
791     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8,
792     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9,
793     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10,
794     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11,
795     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12,
796     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13,
797     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14,
798     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15,
799     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16,
800     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17,
801     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18,
802     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19,
803     1L * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * 10 * 11 * 12 * 13 * 14 * 15 * 16 * 17 * 18 * 19 * 20
804   };
805 
806   /**
807    * Returns {@code n} choose {@code k}, also known as the binomial coefficient of {@code n} and
808    * {@code k}, or {@link Long#MAX_VALUE} if the result does not fit in a {@code long}.
809    *
810    * @throws IllegalArgumentException if {@code n < 0}, {@code k < 0}, or {@code k > n}
811    */
binomial(int n, int k)812   public static long binomial(int n, int k) {
813     checkNonNegative("n", n);
814     checkNonNegative("k", k);
815     checkArgument(k <= n, "k (%s) > n (%s)", k, n);
816     if (k > (n >> 1)) {
817       k = n - k;
818     }
819     switch (k) {
820       case 0:
821         return 1;
822       case 1:
823         return n;
824       default:
825         if (n < factorials.length) {
826           return factorials[n] / (factorials[k] * factorials[n - k]);
827         } else if (k >= biggestBinomials.length || n > biggestBinomials[k]) {
828           return Long.MAX_VALUE;
829         } else if (k < biggestSimpleBinomials.length && n <= biggestSimpleBinomials[k]) {
830           // guaranteed not to overflow
831           long result = n--;
832           for (int i = 2; i <= k; n--, i++) {
833             result *= n;
834             result /= i;
835           }
836           return result;
837         } else {
838           int nBits = LongMath.log2(n, RoundingMode.CEILING);
839 
840           long result = 1;
841           long numerator = n--;
842           long denominator = 1;
843 
844           int numeratorBits = nBits;
845           // This is an upper bound on log2(numerator, ceiling).
846 
847           /*
848            * We want to do this in long math for speed, but want to avoid overflow. We adapt the
849            * technique previously used by BigIntegerMath: maintain separate numerator and
850            * denominator accumulators, multiplying the fraction into result when near overflow.
851            */
852           for (int i = 2; i <= k; i++, n--) {
853             if (numeratorBits + nBits < Long.SIZE - 1) {
854               // It's definitely safe to multiply into numerator and denominator.
855               numerator *= n;
856               denominator *= i;
857               numeratorBits += nBits;
858             } else {
859               // It might not be safe to multiply into numerator and denominator,
860               // so multiply (numerator / denominator) into result.
861               result = multiplyFraction(result, numerator, denominator);
862               numerator = n;
863               denominator = i;
864               numeratorBits = nBits;
865             }
866           }
867           return multiplyFraction(result, numerator, denominator);
868         }
869     }
870   }
871 
872   /** Returns (x * numerator / denominator), which is assumed to come out to an integral value. */
multiplyFraction(long x, long numerator, long denominator)873   static long multiplyFraction(long x, long numerator, long denominator) {
874     if (x == 1) {
875       return numerator / denominator;
876     }
877     long commonDivisor = gcd(x, denominator);
878     x /= commonDivisor;
879     denominator /= commonDivisor;
880     // We know gcd(x, denominator) = 1, and x * numerator / denominator is exact,
881     // so denominator must be a divisor of numerator.
882     return x * (numerator / denominator);
883   }
884 
885   /*
886    * binomial(biggestBinomials[k], k) fits in a long, but not binomial(biggestBinomials[k] + 1, k).
887    */
888   static final int[] biggestBinomials = {
889     Integer.MAX_VALUE,
890     Integer.MAX_VALUE,
891     Integer.MAX_VALUE,
892     3810779,
893     121977,
894     16175,
895     4337,
896     1733,
897     887,
898     534,
899     361,
900     265,
901     206,
902     169,
903     143,
904     125,
905     111,
906     101,
907     94,
908     88,
909     83,
910     79,
911     76,
912     74,
913     72,
914     70,
915     69,
916     68,
917     67,
918     67,
919     66,
920     66,
921     66,
922     66
923   };
924 
925   /*
926    * binomial(biggestSimpleBinomials[k], k) doesn't need to use the slower GCD-based impl, but
927    * binomial(biggestSimpleBinomials[k] + 1, k) does.
928    */
929   @VisibleForTesting
930   static final int[] biggestSimpleBinomials = {
931     Integer.MAX_VALUE,
932     Integer.MAX_VALUE,
933     Integer.MAX_VALUE,
934     2642246,
935     86251,
936     11724,
937     3218,
938     1313,
939     684,
940     419,
941     287,
942     214,
943     169,
944     139,
945     119,
946     105,
947     95,
948     87,
949     81,
950     76,
951     73,
952     70,
953     68,
954     66,
955     64,
956     63,
957     62,
958     62,
959     61,
960     61,
961     61
962   };
963   // These values were generated by using checkedMultiply to see when the simple multiply/divide
964   // algorithm would lead to an overflow.
965 
fitsInInt(long x)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    */
mean(long x, long y)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
isPrime(long n)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
mulMod(long a, long b, long m)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
squareMod(long a, long m)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}. */
plusMod(long a, long b, long m)1112       private long plusMod(long a, long b, long m) {
1113         return (a >= m - b) ? (a + b - m) : (a + b);
1114       }
1115 
1116       /** Returns (a * 2^32) mod m. a may be any unsigned long. */
times2ToThe32Mod(long a, long m)1117       private long times2ToThe32Mod(long a, long m) {
1118         int remainingPowersOf2 = 32;
1119         do {
1120           int shift = min(remainingPowersOf2, Long.numberOfLeadingZeros(a));
1121           // shift is either the number of powers of 2 left to multiply a by, or the biggest shift
1122           // possible while keeping a in an unsigned long.
1123           a = UnsignedLongs.remainder(a << shift, m);
1124           remainingPowersOf2 -= shift;
1125         } while (remainingPowersOf2 > 0);
1126         return a;
1127       }
1128 
1129       @Override
mulMod(long a, long b, long m)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
squareMod(long a, long m)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 
test(long base, long n)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. */
powMod(long a, long p, long 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. */
testWitness(long base, long n)1204     private boolean testWitness(long base, long n) {
1205       int r = Long.numberOfTrailingZeros(n - 1);
1206       long d = (n - 1) >> r;
1207       base %= n;
1208       if (base == 0) {
1209         return true;
1210       }
1211       // Calculate a := base^d mod n.
1212       long a = powMod(base, d, n);
1213       // n passes this test if
1214       //    base^d = 1 (mod n)
1215       // or base^(2^j * d) = -1 (mod n) for some 0 <= j < r.
1216       if (a == 1) {
1217         return true;
1218       }
1219       int j = 0;
1220       while (a != n - 1) {
1221         if (++j == r) {
1222           return false;
1223         }
1224         a = squareMod(a, n);
1225       }
1226       return true;
1227     }
1228   }
1229 
1230   /**
1231    * Returns {@code x}, rounded to a {@code double} with the specified rounding mode. If {@code x}
1232    * is precisely representable as a {@code double}, its {@code double} value will be returned;
1233    * otherwise, the rounding will choose between the two nearest representable values with {@code
1234    * mode}.
1235    *
1236    * <p>For the case of {@link RoundingMode#HALF_EVEN}, this implementation uses the IEEE 754
1237    * default rounding mode: if the two nearest representable values are equally near, the one with
1238    * the least significant bit zero is chosen. (In such cases, both of the nearest representable
1239    * values are even integers; this method returns the one that is a multiple of a greater power of
1240    * two.)
1241    *
1242    * @throws ArithmeticException if {@code mode} is {@link RoundingMode#UNNECESSARY} and {@code x}
1243    *     is not precisely representable as a {@code double}
1244    * @since 30.0
1245    */
1246   @SuppressWarnings("deprecation")
1247   @GwtIncompatible
roundToDouble(long x, RoundingMode mode)1248   public static double roundToDouble(long x, RoundingMode mode) {
1249     // Logic adapted from ToDoubleRounder.
1250     double roundArbitrarily = (double) x;
1251     long roundArbitrarilyAsLong = (long) roundArbitrarily;
1252     int cmpXToRoundArbitrarily;
1253 
1254     if (roundArbitrarilyAsLong == Long.MAX_VALUE) {
1255       /*
1256        * For most values, the conversion from roundArbitrarily to roundArbitrarilyAsLong is
1257        * lossless. In that case we can compare x to roundArbitrarily using Long.compare(x,
1258        * roundArbitrarilyAsLong). The exception is for values where the conversion to double rounds
1259        * up to give roundArbitrarily equal to 2^63, so the conversion back to long overflows and
1260        * roundArbitrarilyAsLong is Long.MAX_VALUE. (This is the only way this condition can occur as
1261        * otherwise the conversion back to long pads with zero bits.) In this case we know that
1262        * roundArbitrarily > x. (This is important when x == Long.MAX_VALUE ==
1263        * roundArbitrarilyAsLong.)
1264        */
1265       cmpXToRoundArbitrarily = -1;
1266     } else {
1267       cmpXToRoundArbitrarily = Long.compare(x, roundArbitrarilyAsLong);
1268     }
1269 
1270     switch (mode) {
1271       case UNNECESSARY:
1272         checkRoundingUnnecessary(cmpXToRoundArbitrarily == 0);
1273         return roundArbitrarily;
1274       case FLOOR:
1275         return (cmpXToRoundArbitrarily >= 0)
1276             ? roundArbitrarily
1277             : DoubleUtils.nextDown(roundArbitrarily);
1278       case CEILING:
1279         return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1280       case DOWN:
1281         if (x >= 0) {
1282           return (cmpXToRoundArbitrarily >= 0)
1283               ? roundArbitrarily
1284               : DoubleUtils.nextDown(roundArbitrarily);
1285         } else {
1286           return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1287         }
1288       case UP:
1289         if (x >= 0) {
1290           return (cmpXToRoundArbitrarily <= 0) ? roundArbitrarily : Math.nextUp(roundArbitrarily);
1291         } else {
1292           return (cmpXToRoundArbitrarily >= 0)
1293               ? roundArbitrarily
1294               : DoubleUtils.nextDown(roundArbitrarily);
1295         }
1296       case HALF_DOWN:
1297       case HALF_UP:
1298       case HALF_EVEN:
1299         {
1300           long roundFloor;
1301           double roundFloorAsDouble;
1302           long roundCeiling;
1303           double roundCeilingAsDouble;
1304 
1305           if (cmpXToRoundArbitrarily >= 0) {
1306             roundFloorAsDouble = roundArbitrarily;
1307             roundFloor = roundArbitrarilyAsLong;
1308             roundCeilingAsDouble = Math.nextUp(roundArbitrarily);
1309             roundCeiling = (long) Math.ceil(roundCeilingAsDouble);
1310           } else {
1311             roundCeilingAsDouble = roundArbitrarily;
1312             roundCeiling = roundArbitrarilyAsLong;
1313             roundFloorAsDouble = DoubleUtils.nextDown(roundArbitrarily);
1314             roundFloor = (long) Math.floor(roundFloorAsDouble);
1315           }
1316 
1317           long deltaToFloor = x - roundFloor;
1318           long deltaToCeiling = roundCeiling - x;
1319 
1320           if (roundCeiling == Long.MAX_VALUE) {
1321             // correct for Long.MAX_VALUE as discussed above: roundCeilingAsDouble must be 2^63, but
1322             // roundCeiling is 2^63-1.
1323             deltaToCeiling++;
1324           }
1325 
1326           int diff = Long.compare(deltaToFloor, deltaToCeiling);
1327           if (diff < 0) { // closer to floor
1328             return roundFloorAsDouble;
1329           } else if (diff > 0) { // closer to ceiling
1330             return roundCeilingAsDouble;
1331           }
1332           // halfway between the representable values; do the half-whatever logic
1333           switch (mode) {
1334             case HALF_EVEN:
1335               return ((DoubleUtils.getSignificand(roundFloorAsDouble) & 1L) == 0)
1336                   ? roundFloorAsDouble
1337                   : roundCeilingAsDouble;
1338             case HALF_DOWN:
1339               return (x >= 0) ? roundFloorAsDouble : roundCeilingAsDouble;
1340             case HALF_UP:
1341               return (x >= 0) ? roundCeilingAsDouble : roundFloorAsDouble;
1342             default:
1343               throw new AssertionError("impossible");
1344           }
1345         }
1346     }
1347     throw new AssertionError("impossible");
1348   }
1349 
LongMath()1350   private LongMath() {}
1351 }
1352