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;
020
021/**
022 * Some useful, arithmetics related, additions to the built-in functions in
023 * {@link Math}.
024 *
025 */
026public final class ArithmeticUtils {
027
028    /** Negative exponent exception message part 1. */
029    private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({";
030    /** Negative exponent exception message part 2. */
031    private static final String NEGATIVE_EXPONENT_2 = "})";
032
033    /** Private constructor. */
034    private ArithmeticUtils() {
035        // intentionally empty.
036    }
037
038    /**
039     * Computes the greatest common divisor of the absolute value of two
040     * numbers, using a modified version of the "binary gcd" method.
041     * See Knuth 4.5.2 algorithm B.
042     * The algorithm is due to Josef Stein (1961).
043     * <br>
044     * Special cases:
045     * <ul>
046     *  <li>The invocations
047     *   {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
048     *   {@code gcd(Integer.MIN_VALUE, 0)} and
049     *   {@code gcd(0, Integer.MIN_VALUE)} throw an
050     *   {@code ArithmeticException}, because the result would be 2^31, which
051     *   is too large for an int value.</li>
052     *  <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
053     *   {@code gcd(x, 0)} is the absolute value of {@code x}, except
054     *   for the special cases above.</li>
055     *  <li>The invocation {@code gcd(0, 0)} is the only one which returns
056     *   {@code 0}.</li>
057     * </ul>
058     *
059     * <p>Two numbers are relatively prime, or coprime, if their gcd is 1.</p>
060     *
061     * @param p Number.
062     * @param q Number.
063     * @return the greatest common divisor (never negative).
064     * @throws ArithmeticException if the result cannot be represented as
065     * a non-negative {@code int} value.
066     */
067    public static int gcd(int p, int q) {
068        // Perform the gcd algorithm on negative numbers, so that -2^31 does not
069        // need to be handled separately
070        int a = p > 0 ? -p : p;
071        int b = q > 0 ? -q : q;
072
073        final int negatedGcd;
074        if (a == 0) {
075            negatedGcd = b;
076        } else if (b == 0) {
077            negatedGcd = a;
078        } else {
079            // Make "a" and "b" odd, keeping track of common power of 2.
080            final int aTwos = Integer.numberOfTrailingZeros(a);
081            final int bTwos = Integer.numberOfTrailingZeros(b);
082            a >>= aTwos;
083            b >>= bTwos;
084            final int shift = Math.min(aTwos, bTwos);
085
086            // "a" and "b" are negative and odd.
087            // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
088            // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
089            // Hence, in the successive iterations:
090            //  "a" becomes the negative absolute difference of the current values,
091            //  "b" becomes that value of the two that is closer to zero.
092            while (a != b) {
093                final int delta = a - b;
094                b = Math.max(a, b);
095                a = delta > 0 ? -delta : delta;
096
097                // Remove any power of 2 in "a" ("b" is guaranteed to be odd).
098                a >>= Integer.numberOfTrailingZeros(a);
099            }
100
101            // Recover the common power of 2.
102            negatedGcd = a << shift;
103        }
104        if (negatedGcd == Integer.MIN_VALUE) {
105            throw new NumbersArithmeticException("overflow: gcd(%d, %d) is 2^31",
106                                                 p, q);
107        }
108        return -negatedGcd;
109    }
110
111    /**
112     * <p>
113     * Gets the greatest common divisor of the absolute value of two numbers,
114     * using the "binary gcd" method which avoids division and modulo
115     * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
116     * Stein (1961).
117     * </p>
118     * Special cases:
119     * <ul>
120     * <li>The invocations
121     * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)},
122     * {@code gcd(Long.MIN_VALUE, 0L)} and
123     * {@code gcd(0L, Long.MIN_VALUE)} throw an
124     * {@code ArithmeticException}, because the result would be 2^63, which
125     * is too large for a long value.</li>
126     * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and
127     * {@code gcd(x, 0L)} is the absolute value of {@code x}, except
128     * for the special cases above.</li>
129     * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns
130     * {@code 0L}.</li>
131     * </ul>
132     *
133     * <p>Two numbers are relatively prime, or coprime, if their gcd is 1.</p>
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(long p, long q) {
142        // Perform the gcd algorithm on negative numbers, so that -2^63 does not
143        // need to be handled separately
144        long a = p > 0 ? -p : p;
145        long b = q > 0 ? -q : q;
146
147        final long negatedGcd;
148        if (a == 0) {
149            negatedGcd = b;
150        } else if (b == 0) {
151            negatedGcd = a;
152        } else {
153            // Make "a" and "b" odd, keeping track of common power of 2.
154            final int aTwos = Long.numberOfTrailingZeros(a);
155            final int bTwos = Long.numberOfTrailingZeros(b);
156            a >>= aTwos;
157            b >>= bTwos;
158            final int shift = Math.min(aTwos, bTwos);
159
160            // "a" and "b" are negative and odd.
161            // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
162            // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
163            // Hence, in the successive iterations:
164            //  "a" becomes the negative absolute difference of the current values,
165            //  "b" becomes that value of the two that is closer to zero.
166            while (true) {
167                final long delta = a - b;
168
169                if (delta == 0) {
170                    // This way of terminating the loop is intentionally different from the int gcd implementation.
171                    // Benchmarking shows that testing for long inequality (a != b) is slow compared to
172                    // testing the delta against zero. The same change on the int gcd reduces performance there,
173                    // hence we have two variants of this loop.
174                    break;
175                }
176
177                b = Math.max(a, b);
178                a = delta > 0 ? -delta : delta;
179
180                // Remove any power of 2 in "a" ("b" is guaranteed to be odd).
181                a >>= Long.numberOfTrailingZeros(a);
182            }
183
184            // Recover the common power of 2.
185            negatedGcd = a << shift;
186        }
187        if (negatedGcd == Long.MIN_VALUE) {
188            throw new NumbersArithmeticException("overflow: gcd(%d, %d) is 2^63",
189                    p, q);
190        }
191        return -negatedGcd;
192    }
193
194    /**
195     * <p>
196     * Returns the least common multiple of the absolute value of two numbers,
197     * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
198     * </p>
199     * Special cases:
200     * <ul>
201     * <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and
202     * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
203     * power of 2, throw an {@code ArithmeticException}, because the result
204     * would be 2^31, which is too large for an int value.</li>
205     * <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
206     * {@code 0} for any {@code x}.</li>
207     * </ul>
208     *
209     * @param a Number.
210     * @param b Number.
211     * @return the least common multiple, never negative.
212     * @throws ArithmeticException if the result cannot be represented as
213     * a non-negative {@code int} value.
214     */
215    public static int lcm(int a, int b) {
216        if (a == 0 || b == 0) {
217            return 0;
218        }
219        final int lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
220        if (lcm == Integer.MIN_VALUE) {
221            throw new NumbersArithmeticException("overflow: lcm(%d, %d) is 2^31",
222                                                 a, b);
223        }
224        return lcm;
225    }
226
227    /**
228     * <p>
229     * Returns the least common multiple of the absolute value of two numbers,
230     * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
231     * </p>
232     * Special cases:
233     * <ul>
234     * <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and
235     * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
236     * power of 2, throw an {@code ArithmeticException}, because the result
237     * would be 2^63, which is too large for an int value.</li>
238     * <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
239     * {@code 0L} for any {@code x}.</li>
240     * </ul>
241     *
242     * @param a Number.
243     * @param b Number.
244     * @return the least common multiple, never negative.
245     * @throws ArithmeticException if the result cannot be represented
246     * as a non-negative {@code long} value.
247     */
248    public static long lcm(long a, long b) {
249        if (a == 0 || b == 0) {
250            return 0;
251        }
252        final long lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
253        if (lcm == Long.MIN_VALUE) {
254            throw new NumbersArithmeticException("overflow: lcm(%d, %d) is 2^63",
255                                                 a, b);
256        }
257        return lcm;
258    }
259
260    /**
261     * Raise an int to an int power.
262     *
263     * <p>Special cases:</p>
264     * <ul>
265     *   <li>{@code k^0} returns {@code 1} (including {@code k=0})</li>
266     *   <li>{@code k^1} returns {@code k} (including {@code k=0})</li>
267     *   <li>{@code 0^0} returns {@code 1}</li>
268     *   <li>{@code 0^e} returns {@code 0}</li>
269     *   <li>{@code 1^e} returns {@code 1}</li>
270     *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even</li>
271     * </ul>
272     *
273     * @param k Number to raise.
274     * @param e Exponent (must be positive or zero).
275     * @return \( k^e \)
276     * @throws IllegalArgumentException if {@code e < 0}.
277     * @throws ArithmeticException if the result would overflow.
278     */
279    public static int pow(final int k,
280                          final int e) {
281        if (e < 0) {
282            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
283        }
284
285        if (k == 0) {
286            return e == 0 ? 1 : 0;
287        }
288
289        if (k == 1) {
290            return 1;
291        }
292
293        if (k == -1) {
294            return (e & 1) == 0 ? 1 : -1;
295        }
296
297        if (e >= 31) {
298            throw new ArithmeticException("integer overflow");
299        }
300
301        int exp = e;
302        int result = 1;
303        int k2p    = k;
304        while (true) {
305            if ((exp & 0x1) != 0) {
306                result = Math.multiplyExact(result, k2p);
307            }
308
309            exp >>= 1;
310            if (exp == 0) {
311                break;
312            }
313
314            k2p = Math.multiplyExact(k2p, k2p);
315        }
316
317        return result;
318    }
319
320    /**
321     * Raise a long to an int power.
322     *
323     * <p>Special cases:</p>
324     * <ul>
325     *   <li>{@code k^0} returns {@code 1} (including {@code k=0})</li>
326     *   <li>{@code k^1} returns {@code k} (including {@code k=0})</li>
327     *   <li>{@code 0^0} returns {@code 1}</li>
328     *   <li>{@code 0^e} returns {@code 0}</li>
329     *   <li>{@code 1^e} returns {@code 1}</li>
330     *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even</li>
331     * </ul>
332     *
333     * @param k Number to raise.
334     * @param e Exponent (must be positive or zero).
335     * @return \( k^e \)
336     * @throws IllegalArgumentException if {@code e < 0}.
337     * @throws ArithmeticException if the result would overflow.
338     */
339    public static long pow(final long k,
340                           final int e) {
341        if (e < 0) {
342            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
343        }
344
345        if (k == 0L) {
346            return e == 0 ? 1L : 0L;
347        }
348
349        if (k == 1L) {
350            return 1L;
351        }
352
353        if (k == -1L) {
354            return (e & 1) == 0 ? 1L : -1L;
355        }
356
357        if (e >= 63) {
358            throw new ArithmeticException("long overflow");
359        }
360
361        int exp = e;
362        long result = 1;
363        long k2p    = k;
364        while (true) {
365            if ((exp & 0x1) != 0) {
366                result = Math.multiplyExact(result, k2p);
367            }
368
369            exp >>= 1;
370            if (exp == 0) {
371                break;
372            }
373
374            k2p = Math.multiplyExact(k2p, k2p);
375        }
376
377        return result;
378    }
379
380    /**
381     * Raise a BigInteger to an int power.
382     *
383     * @param k Number to raise.
384     * @param e Exponent (must be positive or zero).
385     * @return k<sup>e</sup>
386     * @throws IllegalArgumentException if {@code e < 0}.
387     */
388    public static BigInteger pow(final BigInteger k, int e) {
389        if (e < 0) {
390            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
391        }
392
393        return k.pow(e);
394    }
395
396    /**
397     * Raise a BigInteger to a long power.
398     *
399     * @param k Number to raise.
400     * @param e Exponent (must be positive or zero).
401     * @return k<sup>e</sup>
402     * @throws IllegalArgumentException if {@code e < 0}.
403     */
404    public static BigInteger pow(final BigInteger k, final long e) {
405        if (e < 0) {
406            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
407        }
408
409        long exp = e;
410        BigInteger result = BigInteger.ONE;
411        BigInteger k2p    = k;
412        while (exp != 0) {
413            if ((exp & 0x1) != 0) {
414                result = result.multiply(k2p);
415            }
416            k2p = k2p.multiply(k2p);
417            exp >>= 1;
418        }
419
420        return result;
421    }
422
423    /**
424     * Raise a BigInteger to a BigInteger power.
425     *
426     * @param k Number to raise.
427     * @param e Exponent (must be positive or zero).
428     * @return k<sup>e</sup>
429     * @throws IllegalArgumentException if {@code e < 0}.
430     */
431    public static BigInteger pow(final BigInteger k, final BigInteger e) {
432        if (e.compareTo(BigInteger.ZERO) < 0) {
433            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
434        }
435
436        BigInteger exp = e;
437        BigInteger result = BigInteger.ONE;
438        BigInteger k2p    = k;
439        while (!BigInteger.ZERO.equals(exp)) {
440            if (exp.testBit(0)) {
441                result = result.multiply(k2p);
442            }
443            k2p = k2p.multiply(k2p);
444            exp = exp.shiftRight(1);
445        }
446
447        return result;
448    }
449
450    /**
451     * Returns true if the argument is a power of two.
452     *
453     * @param n the number to test
454     * @return true if the argument is a power of two
455     */
456    public static boolean isPowerOfTwo(long n) {
457        return n > 0 && (n & (n - 1)) == 0;
458    }
459
460    /**
461     * Returns the unsigned remainder from dividing the first argument
462     * by the second where each argument and the result is interpreted
463     * as an unsigned value.
464     *
465     * <p>Implementation note
466     *
467     * <p>In v1.0 this method did not use the {@code long} datatype.
468     * Modern 64-bit processors make use of the {@code long} datatype
469     * faster than an algorithm using the {@code int} datatype. This method
470     * now delegates to {@link Integer#remainderUnsigned(int, int)}
471     * which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
472     *
473     * @param dividend the value to be divided
474     * @param divisor the value doing the dividing
475     * @return the unsigned remainder of the first argument divided by
476     * the second argument.
477     * @see Integer#remainderUnsigned(int, int)
478     */
479    public static int remainderUnsigned(int dividend, int divisor) {
480        return Integer.remainderUnsigned(dividend, divisor);
481    }
482
483    /**
484     * Returns the unsigned remainder from dividing the first argument
485     * by the second where each argument and the result is interpreted
486     * as an unsigned value.
487     *
488     * <p>Implementation note
489     *
490     * <p>This method does not use the {@code BigInteger} datatype.
491     * The JDK implementation of {@link Long#remainderUnsigned(long, long)}
492     * uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
493     * From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
494     * even faster due to use of an intrinsic method.
495     *
496     * @param dividend the value to be divided
497     * @param divisor the value doing the dividing
498     * @return the unsigned remainder of the first argument divided by
499     * the second argument.
500     * @see Long#remainderUnsigned(long, long)
501     */
502    public static long remainderUnsigned(long dividend, long divisor) {
503        // Adapts the divideUnsigned method to compute the remainder.
504        if (divisor < 0) {
505            // Using unsigned compare:
506            // if dividend < divisor: return dividend
507            // else: return dividend - divisor
508
509            // Subtracting divisor using masking is more complex in this case
510            // and we use a condition
511            return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor;
512        }
513        // From Hacker's Delight 2.0, section 9.3
514        final long q = ((dividend >>> 1) / divisor) << 1;
515        final long r = dividend - q * divisor;
516        // unsigned r: 0 <= r < 2 * divisor
517        // if (r < divisor): r
518        // else: r - divisor
519
520        // The compare of unsigned r can be done using:
521        // return (r + Long.MIN_VALUE) < (divisor | Long.MIN_VALUE) ? r : r - divisor
522
523        // Here we subtract divisor if (r - divisor) is positive, else the result is r.
524        // This can be done by flipping the sign bit and
525        // creating a mask as -1 or 0 by signed shift.
526        return r - (divisor & (~(r - divisor) >> 63));
527    }
528
529    /**
530     * Returns the unsigned quotient of dividing the first argument by
531     * the second where each argument and the result is interpreted as
532     * an unsigned value.
533     * <p>Note that in two's complement arithmetic, the three other
534     * basic arithmetic operations of add, subtract, and multiply are
535     * bit-wise identical if the two operands are regarded as both
536     * being signed or both being unsigned. Therefore separate {@code
537     * addUnsigned}, etc. methods are not provided.</p>
538     *
539     * <p>Implementation note
540     *
541     * <p>In v1.0 this method did not use the {@code long} datatype.
542     * Modern 64-bit processors make use of the {@code long} datatype
543     * faster than an algorithm using the {@code int} datatype. This method
544     * now delegates to {@link Integer#divideUnsigned(int, int)}
545     * which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
546     *
547     * @param dividend the value to be divided
548     * @param divisor the value doing the dividing
549     * @return the unsigned quotient of the first argument divided by
550     * the second argument
551     * @see Integer#divideUnsigned(int, int)
552     */
553    public static int divideUnsigned(int dividend, int divisor) {
554        return Integer.divideUnsigned(dividend, divisor);
555    }
556
557    /**
558     * Returns the unsigned quotient of dividing the first argument by
559     * the second where each argument and the result is interpreted as
560     * an unsigned value.
561     * <p>Note that in two's complement arithmetic, the three other
562     * basic arithmetic operations of add, subtract, and multiply are
563     * bit-wise identical if the two operands are regarded as both
564     * being signed or both being unsigned. Therefore separate {@code
565     * addUnsigned}, etc. methods are not provided.</p>
566     *
567     * <p>Implementation note
568     *
569     * <p>This method does not use the {@code BigInteger} datatype.
570     * The JDK implementation of {@link Long#divideUnsigned(long, long)}
571     * uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
572     * From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
573     * even faster due to use of an intrinsic method.
574     *
575     * @param dividend the value to be divided
576     * @param divisor the value doing the dividing
577     * @return the unsigned quotient of the first argument divided by
578     * the second argument.
579     * @see Long#divideUnsigned(long, long)
580     */
581    public static long divideUnsigned(long dividend, long divisor) {
582        // The implementation is a Java port of algorithm described in the book
583        // "Hacker's Delight 2.0" (section 9.3 "Unsigned short division from signed division").
584        // Adapts 6-line predicate expressions program with (u >=) an unsigned compare
585        // using the provided branchless variants.
586        if (divisor < 0) {
587            // line 1 branchless:
588            // q <- (dividend (u >=) divisor)
589            return (dividend & ~(dividend - divisor)) >>> 63;
590        }
591        final long q = ((dividend >>> 1) / divisor) << 1;
592        final long r = dividend - q * divisor;
593        // line 5 branchless:
594        // q <- q + (r (u >=) divisor)
595        return q + ((r | ~(r - divisor)) >>> 63);
596    }
597
598    /**
599     * Exception.
600     */
601    private static class NumbersArithmeticException extends ArithmeticException {
602        /** Serializable version Id. */
603        private static final long serialVersionUID = 20180130L;
604
605        /**
606         * Create an exception where the message is constructed by applying
607         * {@link String#format(String, Object...)}.
608         *
609         * @param message Exception message format string
610         * @param args Arguments for formatting the message
611         */
612        NumbersArithmeticException(String message, Object... args) {
613            super(String.format(message, args));
614        }
615    }
616}