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        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.
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        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}.
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}.
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})
266     *   <li>{@code k^1} returns {@code k} (including {@code k=0})
267     *   <li>{@code 0^0} returns {@code 1}
268     *   <li>{@code 0^e} returns {@code 0}
269     *   <li>{@code 1^e} returns {@code 1}
270     *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
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})
326     *   <li>{@code k^1} returns {@code k} (including {@code k=0})
327     *   <li>{@code 0^0} returns {@code 1}
328     *   <li>{@code 0^e} returns {@code 0}
329     *   <li>{@code 1^e} returns {@code 1}
330     *   <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
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    /**
425     * Raise a BigInteger to a BigInteger power.
426     *
427     * @param k Number to raise.
428     * @param e Exponent (must be positive or zero).
429     * @return k<sup>e</sup>
430     * @throws IllegalArgumentException if {@code e < 0}.
431     */
432    public static BigInteger pow(final BigInteger k, final BigInteger e) {
433        if (e.compareTo(BigInteger.ZERO) < 0) {
434            throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
435        }
436
437        BigInteger exp = e;
438        BigInteger result = BigInteger.ONE;
439        BigInteger k2p    = k;
440        while (!BigInteger.ZERO.equals(exp)) {
441            if (exp.testBit(0)) {
442                result = result.multiply(k2p);
443            }
444            k2p = k2p.multiply(k2p);
445            exp = exp.shiftRight(1);
446        }
447
448        return result;
449    }
450
451    /**
452     * Returns true if the argument is a power of two.
453     *
454     * @param n the number to test
455     * @return true if the argument is a power of two
456     */
457    public static boolean isPowerOfTwo(long n) {
458        return n > 0 && (n & (n - 1)) == 0;
459    }
460
461    /**
462     * Returns the unsigned remainder from dividing the first argument
463     * by the second where each argument and the result is interpreted
464     * as an unsigned value.
465     *
466     * <p>Implementation note
467     *
468     * <p>In v1.0 this method did not use the {@code long} datatype.
469     * Modern 64-bit processors make use of the {@code long} datatype
470     * faster than an algorithm using the {@code int} datatype. This method
471     * now delegates to {@link Integer#remainderUnsigned(int, int)}
472     * which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
473     *
474     * @param dividend the value to be divided
475     * @param divisor the value doing the dividing
476     * @return the unsigned remainder of the first argument divided by
477     * the second argument.
478     * @see Integer#remainderUnsigned(int, int)
479     */
480    public static int remainderUnsigned(int dividend, int divisor) {
481        return Integer.remainderUnsigned(dividend, divisor);
482    }
483
484    /**
485     * Returns the unsigned remainder from dividing the first argument
486     * by the second where each argument and the result is interpreted
487     * as an unsigned value.
488     *
489     * <p>Implementation note
490     *
491     * <p>This method does not use the {@code BigInteger} datatype.
492     * The JDK implementation of {@link Long#remainderUnsigned(long, long)}
493     * uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
494     * From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
495     * even faster due to use of an intrinsic method.
496     *
497     * @param dividend the value to be divided
498     * @param divisor the value doing the dividing
499     * @return the unsigned remainder of the first argument divided by
500     * the second argument.
501     * @see Long#remainderUnsigned(long, long)
502     */
503    public static long remainderUnsigned(long dividend, long divisor) {
504        // Adapts the divideUnsigned method to compute the remainder.
505        if (divisor < 0) {
506            // Using unsigned compare:
507            // if dividend < divisor: return dividend
508            // else: return dividend - divisor
509
510            // Subtracting divisor using masking is more complex in this case
511            // and we use a condition
512            return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor;
513
514        }
515        // From Hacker's Delight 2.0, section 9.3
516        final long q = ((dividend >>> 1) / divisor) << 1;
517        final long r = dividend - q * divisor;
518        // unsigned r: 0 <= r < 2 * divisor
519        // if (r < divisor): r
520        // else: r - divisor
521
522        // The compare of unsigned r can be done using:
523        // return (r + Long.MIN_VALUE) < (divisor | Long.MIN_VALUE) ? r : r - divisor
524
525        // Here we subtract divisor if (r - divisor) is positive, else the result is r.
526        // This can be done by flipping the sign bit and
527        // creating a mask as -1 or 0 by signed shift.
528        return r - (divisor & (~(r - divisor) >> 63));
529    }
530
531    /**
532     * Returns the unsigned quotient of dividing the first argument by
533     * the second where each argument and the result is interpreted as
534     * an unsigned value.
535     * <p>Note that in two's complement arithmetic, the three other
536     * basic arithmetic operations of add, subtract, and multiply are
537     * bit-wise identical if the two operands are regarded as both
538     * being signed or both being unsigned. Therefore separate {@code
539     * addUnsigned}, etc. methods are not provided.</p>
540     *
541     * <p>Implementation note
542     *
543     * <p>In v1.0 this method did not use the {@code long} datatype.
544     * Modern 64-bit processors make use of the {@code long} datatype
545     * faster than an algorithm using the {@code int} datatype. This method
546     * now delegates to {@link Integer#divideUnsigned(int, int)}
547     * which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
548     *
549     * @param dividend the value to be divided
550     * @param divisor the value doing the dividing
551     * @return the unsigned quotient of the first argument divided by
552     * the second argument
553     * @see Integer#divideUnsigned(int, int)
554     */
555    public static int divideUnsigned(int dividend, int divisor) {
556        return Integer.divideUnsigned(dividend, divisor);
557    }
558
559    /**
560     * Returns the unsigned quotient of dividing the first argument by
561     * the second where each argument and the result is interpreted as
562     * an unsigned value.
563     * <p>Note that in two's complement arithmetic, the three other
564     * basic arithmetic operations of add, subtract, and multiply are
565     * bit-wise identical if the two operands are regarded as both
566     * being signed or both being unsigned. Therefore separate {@code
567     * addUnsigned}, etc. methods are not provided.</p>
568     *
569     * <p>Implementation note
570     *
571     * <p>This method does not use the {@code BigInteger} datatype.
572     * The JDK implementation of {@link Long#divideUnsigned(long, long)}
573     * uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
574     * From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
575     * even faster due to use of an intrinsic method.
576     *
577     * @param dividend the value to be divided
578     * @param divisor the value doing the dividing
579     * @return the unsigned quotient of the first argument divided by
580     * the second argument.
581     * @see Long#divideUnsigned(long, long)
582     */
583    public static long divideUnsigned(long dividend, long divisor) {
584        // The implementation is a Java port of algorithm described in the book
585        // "Hacker's Delight 2.0" (section 9.3 "Unsigned short division from signed division").
586        // Adapts 6-line predicate expressions program with (u >=) an unsigned compare
587        // using the provided branchless variants.
588        if (divisor < 0) {
589            // line 1 branchless:
590            // q <- (dividend (u >=) divisor)
591            return (dividend & ~(dividend - divisor)) >>> 63;
592        }
593        final long q = ((dividend >>> 1) / divisor) << 1;
594        final long r = dividend - q * divisor;
595        // line 5 branchless:
596        // q <- q + (r (u >=) divisor)
597        return q + ((r | ~(r - divisor)) >>> 63);
598    }
599
600    /**
601     * Exception.
602     */
603    private static class NumbersArithmeticException extends ArithmeticException {
604        /** Serializable version Id. */
605        private static final long serialVersionUID = 20180130L;
606
607        /**
608         * Create an exception where the message is constructed by applying
609         * {@link String#format(String, Object...)}.
610         *
611         * @param message Exception message format string
612         * @param args Arguments for formatting the message
613         */
614        NumbersArithmeticException(String message, Object... args) {
615            super(String.format(message, args));
616        }
617    }
618}