View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.numbers.core;
18  
19  import java.math.BigInteger;
20  import java.text.MessageFormat;
21  
22  /**
23   * Some useful, arithmetics related, additions to the built-in functions in
24   * {@link Math}.
25   *
26   */
27  public final class ArithmeticUtils {
28  
29      /** Overflow gcd exception message for 2^63. */
30      private static final String OVERFLOW_GCD_MESSAGE_2_POWER_63 = "overflow: gcd({0}, {1}) is 2^63";
31  
32      /** Negative exponent exception message part 1. */
33      private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({";
34      /** Negative exponent exception message part 2. */
35      private static final String NEGATIVE_EXPONENT_2 = "})";
36  
37      /** Private constructor. */
38      private ArithmeticUtils() {
39          // intentionally empty.
40      }
41  
42      /**
43       * Computes the greatest common divisor of the absolute value of two
44       * numbers, using a modified version of the "binary gcd" method.
45       * See Knuth 4.5.2 algorithm B.
46       * The algorithm is due to Josef Stein (1961).
47       * <br>
48       * Special cases:
49       * <ul>
50       *  <li>The invocations
51       *   {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
52       *   {@code gcd(Integer.MIN_VALUE, 0)} and
53       *   {@code gcd(0, Integer.MIN_VALUE)} throw an
54       *   {@code ArithmeticException}, because the result would be 2^31, which
55       *   is too large for an int value.</li>
56       *  <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
57       *   {@code gcd(x, 0)} is the absolute value of {@code x}, except
58       *   for the special cases above.</li>
59       *  <li>The invocation {@code gcd(0, 0)} is the only one which returns
60       *   {@code 0}.</li>
61       * </ul>
62       *
63       * @param p Number.
64       * @param q Number.
65       * @return the greatest common divisor (never negative).
66       * @throws ArithmeticException if the result cannot be represented as
67       * a non-negative {@code int} value.
68       */
69      public static int gcd(int p, int q) {
70          // Perform the gcd algorithm on negative numbers, so that -2^31 does not
71          // need to be handled separately
72          int a = p > 0 ? -p : p;
73          int b = q > 0 ? -q : q;
74  
75          int negatedGcd;
76          if (a == 0) {
77              negatedGcd = b;
78          } else if (b == 0) {
79              negatedGcd = a;
80          } else {
81              // Make "a" and "b" odd, keeping track of common power of 2.
82              final int aTwos = Integer.numberOfTrailingZeros(a);
83              final int bTwos = Integer.numberOfTrailingZeros(b);
84              a >>= aTwos;
85              b >>= bTwos;
86              final int shift = Math.min(aTwos, bTwos);
87  
88              // "a" and "b" are negative and odd.
89              // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
90              // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
91              // Hence, in the successive iterations:
92              //  "a" becomes the negative absolute difference of the current values,
93              //  "b" becomes that value of the two that is closer to zero.
94              while (a != b) {
95                  final int delta = a - b;
96                  b = Math.max(a, b);
97                  a = delta > 0 ? -delta : delta;
98  
99                  // Remove any power of 2 in "a" ("b" is guaranteed to be odd).
100                 a >>= Integer.numberOfTrailingZeros(a);
101             }
102 
103             // Recover the common power of 2.
104             negatedGcd = a << shift;
105         }
106         if (negatedGcd == Integer.MIN_VALUE) {
107             throw new NumbersArithmeticException("overflow: gcd({0}, {1}) is 2^31",
108                                               p, q);
109         }
110         return -negatedGcd;
111     }
112 
113     /**
114      * <p>
115      * Gets the greatest common divisor of the absolute value of two numbers,
116      * using the "binary gcd" method which avoids division and modulo
117      * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
118      * Stein (1961).
119      * </p>
120      * Special cases:
121      * <ul>
122      * <li>The invocations
123      * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)},
124      * {@code gcd(Long.MIN_VALUE, 0L)} and
125      * {@code gcd(0L, Long.MIN_VALUE)} throw an
126      * {@code ArithmeticException}, because the result would be 2^63, which
127      * is too large for a long value.</li>
128      * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and
129      * {@code gcd(x, 0L)} is the absolute value of {@code x}, except
130      * for the special cases above.
131      * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns
132      * {@code 0L}.</li>
133      * </ul>
134      *
135      * @param p Number.
136      * @param q Number.
137      * @return the greatest common divisor, never negative.
138      * @throws ArithmeticException if the result cannot be represented as
139      * a non-negative {@code long} value.
140      */
141     public static long gcd(final long p, final long q) {
142         long u = p;
143         long v = q;
144         if (u == 0 || v == 0) {
145             if (u == Long.MIN_VALUE || v == Long.MIN_VALUE) {
146                 throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63,
147                                                   p, q);
148             }
149             return Math.abs(u) + Math.abs(v);
150         }
151         // keep u and v negative, as negative integers range down to
152         // -2^63, while positive numbers can only be as large as 2^63-1
153         // (i.e. we can't necessarily negate a negative number without
154         // overflow)
155         /* assert u!=0 && v!=0; */
156         if (u > 0) {
157             u = -u;
158         } // make u negative
159         if (v > 0) {
160             v = -v;
161         } // make v negative
162         // B1. [Find power of 2]
163         int k = 0;
164         while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
165                                                             // both even...
166             u /= 2;
167             v /= 2;
168             k++; // cast out twos.
169         }
170         if (k == 63) {
171             throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63,
172                                               p, q);
173         }
174         // B2. Initialize: u and v have been divided by 2^k and at least
175         // one is odd.
176         long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
177         // t negative: u was odd, v may be even (t replaces v)
178         // t positive: u was even, v is odd (t replaces u)
179         do {
180             /* assert u<0 && v<0; */
181             // B4/B3: cast out twos from t.
182             while ((t & 1) == 0) { // while t is even..
183                 t /= 2; // cast out twos
184             }
185             // B5 [reset max(u,v)]
186             if (t > 0) {
187                 u = -t;
188             } else {
189                 v = t;
190             }
191             // B6/B3. at this point both u and v should be odd.
192             t = (v - u) / 2;
193             // |u| larger: t positive (replace u)
194             // |v| larger: t negative (replace v)
195         } while (t != 0);
196         return -u * (1L << k); // gcd is u*2^k
197     }
198 
199     /**
200      * <p>
201      * Returns the least common multiple of the absolute value of two numbers,
202      * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
203      * </p>
204      * Special cases:
205      * <ul>
206      * <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and
207      * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
208      * power of 2, throw an {@code ArithmeticException}, because the result
209      * would be 2^31, which is too large for an int value.</li>
210      * <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
211      * {@code 0} for any {@code x}.
212      * </ul>
213      *
214      * @param a Number.
215      * @param b Number.
216      * @return the least common multiple, never negative.
217      * @throws ArithmeticException if the result cannot be represented as
218      * a non-negative {@code int} value.
219      */
220     public static int lcm(int a, int b) {
221         if (a == 0 || b == 0) {
222             return 0;
223         }
224         final int lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
225         if (lcm == Integer.MIN_VALUE) {
226             throw new NumbersArithmeticException("overflow: lcm({0}, {1}) is 2^31",
227                                               a, b);
228         }
229         return lcm;
230     }
231 
232     /**
233      * <p>
234      * Returns the least common multiple of the absolute value of two numbers,
235      * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
236      * </p>
237      * Special cases:
238      * <ul>
239      * <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and
240      * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
241      * power of 2, throw an {@code ArithmeticException}, because the result
242      * would be 2^63, which is too large for an int value.</li>
243      * <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
244      * {@code 0L} for any {@code x}.
245      * </ul>
246      *
247      * @param a Number.
248      * @param b Number.
249      * @return the least common multiple, never negative.
250      * @throws ArithmeticException if the result cannot be represented
251      * as a non-negative {@code long} value.
252      */
253     public static long lcm(long a, long b) {
254         if (a == 0 || b == 0) {
255             return 0;
256         }
257         final long lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
258         if (lcm == Long.MIN_VALUE) {
259             throw new NumbersArithmeticException("overflow: lcm({0}, {1}) is 2^63",
260                                               a, b);
261         }
262         return lcm;
263     }
264 
265     /**
266      * Raise an int to an int power.
267      *
268      * <p>Special cases:</p>
269      * <ul>
270      *   <li>{@code k^0} returns {@code 1} (including {@code k=0})
271      *   <li>{@code k^1} returns {@code k} (including {@code k=0})
272      *   <li>{@code 0^0} returns {@code 1}
273      *   <li>{@code 0^e} returns {@code 0}
274      *   <li>{@code 1^e} returns {@code 1}
275      *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
276      * </ul>
277      *
278      * @param k Number to raise.
279      * @param e Exponent (must be positive or zero).
280      * @return \( k^e \)
281      * @throws IllegalArgumentException if {@code e < 0}.
282      * @throws ArithmeticException if the result would overflow.
283      */
284     public static int pow(final int k,
285                           final int e) {
286         if (e < 0) {
287             throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
288         }
289 
290         if (k == 0) {
291             return e == 0 ? 1 : 0;
292         }
293 
294         if (k == 1) {
295             return 1;
296         }
297 
298         if (k == -1) {
299             return (e & 1) == 0 ? 1 : -1;
300         }
301 
302         if (e >= 31) {
303             throw new ArithmeticException("integer overflow");
304         }
305 
306         int exp = e;
307         int result = 1;
308         int k2p    = k;
309         while (true) {
310             if ((exp & 0x1) != 0) {
311                 result = Math.multiplyExact(result, k2p);
312             }
313 
314             exp >>= 1;
315             if (exp == 0) {
316                 break;
317             }
318 
319             k2p = Math.multiplyExact(k2p, k2p);
320         }
321 
322         return result;
323     }
324 
325     /**
326      * Raise a long to an int power.
327      *
328      * <p>Special cases:</p>
329      * <ul>
330      *   <li>{@code k^0} returns {@code 1} (including {@code k=0})
331      *   <li>{@code k^1} returns {@code k} (including {@code k=0})
332      *   <li>{@code 0^0} returns {@code 1}
333      *   <li>{@code 0^e} returns {@code 0}
334      *   <li>{@code 1^e} returns {@code 1}
335      *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
336      * </ul>
337      *
338      * @param k Number to raise.
339      * @param e Exponent (must be positive or zero).
340      * @return \( k^e \)
341      * @throws IllegalArgumentException if {@code e < 0}.
342      * @throws ArithmeticException if the result would overflow.
343      */
344     public static long pow(final long k,
345                            final int e) {
346         if (e < 0) {
347             throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
348         }
349 
350         if (k == 0L) {
351             return e == 0 ? 1L : 0L;
352         }
353 
354         if (k == 1L) {
355             return 1L;
356         }
357 
358         if (k == -1L) {
359             return (e & 1) == 0 ? 1L : -1L;
360         }
361 
362         if (e >= 63) {
363             throw new ArithmeticException("long overflow");
364         }
365 
366         int exp = e;
367         long result = 1;
368         long k2p    = k;
369         while (true) {
370             if ((exp & 0x1) != 0) {
371                 result = Math.multiplyExact(result, k2p);
372             }
373 
374             exp >>= 1;
375             if (exp == 0) {
376                 break;
377             }
378 
379             k2p = Math.multiplyExact(k2p, k2p);
380         }
381 
382         return result;
383     }
384 
385     /**
386      * Raise a BigInteger to an int power.
387      *
388      * @param k Number to raise.
389      * @param e Exponent (must be positive or zero).
390      * @return k<sup>e</sup>
391      * @throws IllegalArgumentException if {@code e < 0}.
392      */
393     public static BigInteger pow(final BigInteger k, int e) {
394         if (e < 0) {
395             throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
396         }
397 
398         return k.pow(e);
399     }
400 
401     /**
402      * Raise a BigInteger to a long power.
403      *
404      * @param k Number to raise.
405      * @param e Exponent (must be positive or zero).
406      * @return k<sup>e</sup>
407      * @throws IllegalArgumentException if {@code e < 0}.
408      */
409     public static BigInteger pow(final BigInteger k, final long e) {
410         if (e < 0) {
411             throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
412         }
413 
414         long exp = e;
415         BigInteger result = BigInteger.ONE;
416         BigInteger k2p    = k;
417         while (exp != 0) {
418             if ((exp & 0x1) != 0) {
419                 result = result.multiply(k2p);
420             }
421             k2p = k2p.multiply(k2p);
422             exp >>= 1;
423         }
424 
425         return result;
426 
427     }
428 
429     /**
430      * Raise a BigInteger to a BigInteger power.
431      *
432      * @param k Number to raise.
433      * @param e Exponent (must be positive or zero).
434      * @return k<sup>e</sup>
435      * @throws IllegalArgumentException if {@code e < 0}.
436      */
437     public static BigInteger pow(final BigInteger k, final BigInteger e) {
438         if (e.compareTo(BigInteger.ZERO) < 0) {
439             throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
440         }
441 
442         BigInteger exp = e;
443         BigInteger result = BigInteger.ONE;
444         BigInteger k2p    = k;
445         while (!BigInteger.ZERO.equals(exp)) {
446             if (exp.testBit(0)) {
447                 result = result.multiply(k2p);
448             }
449             k2p = k2p.multiply(k2p);
450             exp = exp.shiftRight(1);
451         }
452 
453         return result;
454     }
455 
456     /**
457      * Returns true if the argument is a power of two.
458      *
459      * @param n the number to test
460      * @return true if the argument is a power of two
461      */
462     public static boolean isPowerOfTwo(long n) {
463         return n > 0 && (n & (n - 1)) == 0;
464     }
465 
466     /**
467      * Returns the unsigned remainder from dividing the first argument
468      * by the second where each argument and the result is interpreted
469      * as an unsigned value.
470      * <p>This method does not use the {@code long} datatype.</p>
471      *
472      * @param dividend the value to be divided
473      * @param divisor the value doing the dividing
474      * @return the unsigned remainder of the first argument divided by
475      * the second argument.
476      */
477     public static int remainderUnsigned(int dividend, int divisor) {
478         if (divisor >= 0) {
479             if (dividend >= 0) {
480                 return dividend % divisor;
481             }
482             // The implementation is a Java port of algorithm described in the book
483             // "Hacker's Delight" (section "Unsigned short division from signed division").
484             final int q = ((dividend >>> 1) / divisor) << 1;
485             dividend -= q * divisor;
486             if (dividend < 0 || dividend >= divisor) {
487                 dividend -= divisor;
488             }
489             return dividend;
490         }
491         return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor;
492     }
493 
494     /**
495      * Returns the unsigned remainder from dividing the first argument
496      * by the second where each argument and the result is interpreted
497      * as an unsigned value.
498      * <p>This method does not use the {@code BigInteger} datatype.</p>
499      *
500      * @param dividend the value to be divided
501      * @param divisor the value doing the dividing
502      * @return the unsigned remainder of the first argument divided by
503      * the second argument.
504      */
505     public static long remainderUnsigned(long dividend, long divisor) {
506         if (divisor >= 0L) {
507             if (dividend >= 0L) {
508                 return dividend % divisor;
509             }
510             // The implementation is a Java port of algorithm described in the book
511             // "Hacker's Delight" (section "Unsigned short division from signed division").
512             final long q = ((dividend >>> 1) / divisor) << 1;
513             dividend -= q * divisor;
514             if (dividend < 0L || dividend >= divisor) {
515                 dividend -= divisor;
516             }
517             return dividend;
518         }
519         return dividend >= 0L || dividend < divisor ? dividend : dividend - divisor;
520     }
521 
522     /**
523      * Returns the unsigned quotient of dividing the first argument by
524      * the second where each argument and the result is interpreted as
525      * an unsigned value.
526      * <p>Note that in two's complement arithmetic, the three other
527      * basic arithmetic operations of add, subtract, and multiply are
528      * bit-wise identical if the two operands are regarded as both
529      * being signed or both being unsigned. Therefore separate {@code
530      * addUnsigned}, etc. methods are not provided.</p>
531      * <p>This method does not use the {@code long} datatype.</p>
532      *
533      * @param dividend the value to be divided
534      * @param divisor the value doing the dividing
535      * @return the unsigned quotient of the first argument divided by
536      * the second argument
537      */
538     public static int divideUnsigned(int dividend, int divisor) {
539         if (divisor >= 0) {
540             if (dividend >= 0) {
541                 return dividend / divisor;
542             }
543             // The implementation is a Java port of algorithm described in the book
544             // "Hacker's Delight" (section "Unsigned short division from signed division").
545             int q = ((dividend >>> 1) / divisor) << 1;
546             dividend -= q * divisor;
547             if (dividend < 0L || dividend >= divisor) {
548                 q++;
549             }
550             return q;
551         }
552         return dividend >= 0 || dividend < divisor ? 0 : 1;
553     }
554 
555     /**
556      * Returns the unsigned quotient of dividing the first argument by
557      * the second where each argument and the result is interpreted as
558      * an unsigned value.
559      * <p>Note that in two's complement arithmetic, the three other
560      * basic arithmetic operations of add, subtract, and multiply are
561      * bit-wise identical if the two operands are regarded as both
562      * being signed or both being unsigned. Therefore separate {@code
563      * addUnsigned}, etc. methods are not provided.</p>
564      * <p>This method does not use the {@code BigInteger} datatype.</p>
565      *
566      * @param dividend the value to be divided
567      * @param divisor the value doing the dividing
568      * @return the unsigned quotient of the first argument divided by
569      * the second argument.
570      */
571     public static long divideUnsigned(long dividend, long divisor) {
572         if (divisor >= 0L) {
573             if (dividend >= 0L) {
574                 return dividend / divisor;
575             }
576             // The implementation is a Java port of algorithm described in the book
577             // "Hacker's Delight" (section "Unsigned short division from signed division").
578             long q = ((dividend >>> 1) / divisor) << 1;
579             dividend -= q * divisor;
580             if (dividend < 0L || dividend >= divisor) {
581                 q++;
582             }
583             return q;
584         }
585         return dividend >= 0L || dividend < divisor ? 0L : 1L;
586     }
587 
588     /**
589      * Exception.
590      */
591     private static class NumbersArithmeticException extends ArithmeticException {
592         /** Serializable version Id. */
593         private static final long serialVersionUID = 20180130L;
594 
595         /**
596          * Constructor with a specific message.
597          *
598          * @param message Message pattern providing the specific context of
599          * the error.
600          * @param args Arguments.
601          */
602         NumbersArithmeticException(String message, Object... args) {
603             super(MessageFormat.format(message, args));
604         }
605     }
606 }