Source code

001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.core;
018
019import java.math.BigInteger;
020import java.text.MessageFormat;
021
022/**
023 * Some useful, arithmetics related, additions to the built-in functions in
024 * {@link Math}.
025 *
026 */
027public final class ArithmeticUtils {
028
029    /** Overflow gcd exception message for 2^63. */
030    private static final String OVERFLOW_GCD_MESSAGE_2_POWER_63 = "overflow: gcd({0}, {1}) is 2^63";
031
032    /** Negative exponent exception message part 1. */
033    private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({";
034    /** Negative exponent exception message part 2. */
035    private static final String NEGATIVE_EXPONENT_2 = "})";
036
037    /** Private constructor. */
038    private ArithmeticUtils() {
039        // intentionally empty.
040    }
041
042    /**
043     * Computes the greatest common divisor of the absolute value of two
044     * numbers, using a modified version of the "binary gcd" method.
045     * See Knuth 4.5.2 algorithm B.
046     * The algorithm is due to Josef Stein (1961).
047     * <br>
048     * Special cases:
049     * <ul>
050     *  <li>The invocations
051     *   {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
052     *   {@code gcd(Integer.MIN_VALUE, 0)} and
053     *   {@code gcd(0, Integer.MIN_VALUE)} throw an
054     *   {@code ArithmeticException}, because the result would be 2^31, which
055     *   is too large for an int value.</li>
056     *  <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
057     *   {@code gcd(x, 0)} is the absolute value of {@code x}, except
058     *   for the special cases above.</li>
059     *  <li>The invocation {@code gcd(0, 0)} is the only one which returns
060     *   {@code 0}.</li>
061     * </ul>
062     *
063     * @param p Number.
064     * @param q Number.
065     * @return the greatest common divisor (never negative).
066     * @throws ArithmeticException if the result cannot be represented as
067     * a non-negative {@code int} value.
068     */
069    public static int gcd(int p, int q) {
070        // Perform the gcd algorithm on negative numbers, so that -2^31 does not
071        // need to be handled separately
072        int a = p > 0 ? -p : p;
073        int b = q > 0 ? -q : q;
074
075        int negatedGcd;
076        if (a == 0) {
077            negatedGcd = b;
078        } else if (b == 0) {
079            negatedGcd = a;
080        } else {
081            // Make "a" and "b" odd, keeping track of common power of 2.
082            final int aTwos = Integer.numberOfTrailingZeros(a);
083            final int bTwos = Integer.numberOfTrailingZeros(b);
084            a >>= aTwos;
085            b >>= bTwos;
086            final int shift = Math.min(aTwos, bTwos);
087
088            // "a" and "b" are negative and odd.
089            // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
090            // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
091            // Hence, in the successive iterations:
092            //  "a" becomes the negative absolute difference of the current values,
093            //  "b" becomes that value of the two that is closer to zero.
094            while (a != b) {
095                final int delta = a - b;
096                b = Math.max(a, b);
097                a = delta > 0 ? -delta : delta;
098
099                // 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}