ArithmeticUtils.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.numbers.core;
import java.math.BigInteger;
/**
* Some useful, arithmetics related, additions to the built-in functions in
* {@link Math}.
*
*/
public final class ArithmeticUtils {
/** Negative exponent exception message part 1. */
private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({";
/** Negative exponent exception message part 2. */
private static final String NEGATIVE_EXPONENT_2 = "})";
/** Private constructor. */
private ArithmeticUtils() {
// intentionally empty.
}
/**
* Computes the greatest common divisor of the absolute value of two
* numbers, using a modified version of the "binary gcd" method.
* See Knuth 4.5.2 algorithm B.
* The algorithm is due to Josef Stein (1961).
* <br>
* Special cases:
* <ul>
* <li>The invocations
* {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
* {@code gcd(Integer.MIN_VALUE, 0)} and
* {@code gcd(0, Integer.MIN_VALUE)} throw an
* {@code ArithmeticException}, because the result would be 2^31, which
* is too large for an int value.</li>
* <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
* {@code gcd(x, 0)} is the absolute value of {@code x}, except
* for the special cases above.</li>
* <li>The invocation {@code gcd(0, 0)} is the only one which returns
* {@code 0}.</li>
* </ul>
*
* <p>Two numbers are relatively prime, or coprime, if their gcd is 1.</p>
*
* @param p Number.
* @param q Number.
* @return the greatest common divisor (never negative).
* @throws ArithmeticException if the result cannot be represented as
* a non-negative {@code int} value.
*/
public static int gcd(int p, int q) {
// Perform the gcd algorithm on negative numbers, so that -2^31 does not
// need to be handled separately
int a = p > 0 ? -p : p;
int b = q > 0 ? -q : q;
int negatedGcd;
if (a == 0) {
negatedGcd = b;
} else if (b == 0) {
negatedGcd = a;
} else {
// Make "a" and "b" odd, keeping track of common power of 2.
final int aTwos = Integer.numberOfTrailingZeros(a);
final int bTwos = Integer.numberOfTrailingZeros(b);
a >>= aTwos;
b >>= bTwos;
final int shift = Math.min(aTwos, bTwos);
// "a" and "b" are negative and odd.
// If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
// If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
// Hence, in the successive iterations:
// "a" becomes the negative absolute difference of the current values,
// "b" becomes that value of the two that is closer to zero.
while (a != b) {
final int delta = a - b;
b = Math.max(a, b);
a = delta > 0 ? -delta : delta;
// Remove any power of 2 in "a" ("b" is guaranteed to be odd).
a >>= Integer.numberOfTrailingZeros(a);
}
// Recover the common power of 2.
negatedGcd = a << shift;
}
if (negatedGcd == Integer.MIN_VALUE) {
throw new NumbersArithmeticException("overflow: gcd(%d, %d) is 2^31",
p, q);
}
return -negatedGcd;
}
/**
* <p>
* Gets the greatest common divisor of the absolute value of two numbers,
* using the "binary gcd" method which avoids division and modulo
* operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
* Stein (1961).
* </p>
* Special cases:
* <ul>
* <li>The invocations
* {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)},
* {@code gcd(Long.MIN_VALUE, 0L)} and
* {@code gcd(0L, Long.MIN_VALUE)} throw an
* {@code ArithmeticException}, because the result would be 2^63, which
* is too large for a long value.</li>
* <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and
* {@code gcd(x, 0L)} is the absolute value of {@code x}, except
* for the special cases above.
* <li>The invocation {@code gcd(0L, 0L)} is the only one which returns
* {@code 0L}.</li>
* </ul>
*
* <p>Two numbers are relatively prime, or coprime, if their gcd is 1.</p>
*
* @param p Number.
* @param q Number.
* @return the greatest common divisor, never negative.
* @throws ArithmeticException if the result cannot be represented as
* a non-negative {@code long} value.
*/
public static long gcd(long p, long q) {
// Perform the gcd algorithm on negative numbers, so that -2^63 does not
// need to be handled separately
long a = p > 0 ? -p : p;
long b = q > 0 ? -q : q;
long negatedGcd;
if (a == 0) {
negatedGcd = b;
} else if (b == 0) {
negatedGcd = a;
} else {
// Make "a" and "b" odd, keeping track of common power of 2.
final int aTwos = Long.numberOfTrailingZeros(a);
final int bTwos = Long.numberOfTrailingZeros(b);
a >>= aTwos;
b >>= bTwos;
final int shift = Math.min(aTwos, bTwos);
// "a" and "b" are negative and odd.
// If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
// If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
// Hence, in the successive iterations:
// "a" becomes the negative absolute difference of the current values,
// "b" becomes that value of the two that is closer to zero.
while (true) {
final long delta = a - b;
if (delta == 0) {
// This way of terminating the loop is intentionally different from the int gcd implementation.
// Benchmarking shows that testing for long inequality (a != b) is slow compared to
// testing the delta against zero. The same change on the int gcd reduces performance there,
// hence we have two variants of this loop.
break;
}
b = Math.max(a, b);
a = delta > 0 ? -delta : delta;
// Remove any power of 2 in "a" ("b" is guaranteed to be odd).
a >>= Long.numberOfTrailingZeros(a);
}
// Recover the common power of 2.
negatedGcd = a << shift;
}
if (negatedGcd == Long.MIN_VALUE) {
throw new NumbersArithmeticException("overflow: gcd(%d, %d) is 2^63",
p, q);
}
return -negatedGcd;
}
/**
* <p>
* Returns the least common multiple of the absolute value of two numbers,
* using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
* </p>
* Special cases:
* <ul>
* <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and
* {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a
* power of 2, throw an {@code ArithmeticException}, because the result
* would be 2^31, which is too large for an int value.</li>
* <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is
* {@code 0} for any {@code x}.
* </ul>
*
* @param a Number.
* @param b Number.
* @return the least common multiple, never negative.
* @throws ArithmeticException if the result cannot be represented as
* a non-negative {@code int} value.
*/
public static int lcm(int a, int b) {
if (a == 0 || b == 0) {
return 0;
}
final int lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
if (lcm == Integer.MIN_VALUE) {
throw new NumbersArithmeticException("overflow: lcm(%d, %d) is 2^31",
a, b);
}
return lcm;
}
/**
* <p>
* Returns the least common multiple of the absolute value of two numbers,
* using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}.
* </p>
* Special cases:
* <ul>
* <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and
* {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a
* power of 2, throw an {@code ArithmeticException}, because the result
* would be 2^63, which is too large for an int value.</li>
* <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is
* {@code 0L} for any {@code x}.
* </ul>
*
* @param a Number.
* @param b Number.
* @return the least common multiple, never negative.
* @throws ArithmeticException if the result cannot be represented
* as a non-negative {@code long} value.
*/
public static long lcm(long a, long b) {
if (a == 0 || b == 0) {
return 0;
}
final long lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b));
if (lcm == Long.MIN_VALUE) {
throw new NumbersArithmeticException("overflow: lcm(%d, %d) is 2^63",
a, b);
}
return lcm;
}
/**
* Raise an int to an int power.
*
* <p>Special cases:</p>
* <ul>
* <li>{@code k^0} returns {@code 1} (including {@code k=0})
* <li>{@code k^1} returns {@code k} (including {@code k=0})
* <li>{@code 0^0} returns {@code 1}
* <li>{@code 0^e} returns {@code 0}
* <li>{@code 1^e} returns {@code 1}
* <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
* </ul>
*
* @param k Number to raise.
* @param e Exponent (must be positive or zero).
* @return \( k^e \)
* @throws IllegalArgumentException if {@code e < 0}.
* @throws ArithmeticException if the result would overflow.
*/
public static int pow(final int k,
final int e) {
if (e < 0) {
throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
}
if (k == 0) {
return e == 0 ? 1 : 0;
}
if (k == 1) {
return 1;
}
if (k == -1) {
return (e & 1) == 0 ? 1 : -1;
}
if (e >= 31) {
throw new ArithmeticException("integer overflow");
}
int exp = e;
int result = 1;
int k2p = k;
while (true) {
if ((exp & 0x1) != 0) {
result = Math.multiplyExact(result, k2p);
}
exp >>= 1;
if (exp == 0) {
break;
}
k2p = Math.multiplyExact(k2p, k2p);
}
return result;
}
/**
* Raise a long to an int power.
*
* <p>Special cases:</p>
* <ul>
* <li>{@code k^0} returns {@code 1} (including {@code k=0})
* <li>{@code k^1} returns {@code k} (including {@code k=0})
* <li>{@code 0^0} returns {@code 1}
* <li>{@code 0^e} returns {@code 0}
* <li>{@code 1^e} returns {@code 1}
* <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even
* </ul>
*
* @param k Number to raise.
* @param e Exponent (must be positive or zero).
* @return \( k^e \)
* @throws IllegalArgumentException if {@code e < 0}.
* @throws ArithmeticException if the result would overflow.
*/
public static long pow(final long k,
final int e) {
if (e < 0) {
throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
}
if (k == 0L) {
return e == 0 ? 1L : 0L;
}
if (k == 1L) {
return 1L;
}
if (k == -1L) {
return (e & 1) == 0 ? 1L : -1L;
}
if (e >= 63) {
throw new ArithmeticException("long overflow");
}
int exp = e;
long result = 1;
long k2p = k;
while (true) {
if ((exp & 0x1) != 0) {
result = Math.multiplyExact(result, k2p);
}
exp >>= 1;
if (exp == 0) {
break;
}
k2p = Math.multiplyExact(k2p, k2p);
}
return result;
}
/**
* Raise a BigInteger to an int power.
*
* @param k Number to raise.
* @param e Exponent (must be positive or zero).
* @return k<sup>e</sup>
* @throws IllegalArgumentException if {@code e < 0}.
*/
public static BigInteger pow(final BigInteger k, int e) {
if (e < 0) {
throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
}
return k.pow(e);
}
/**
* Raise a BigInteger to a long power.
*
* @param k Number to raise.
* @param e Exponent (must be positive or zero).
* @return k<sup>e</sup>
* @throws IllegalArgumentException if {@code e < 0}.
*/
public static BigInteger pow(final BigInteger k, final long e) {
if (e < 0) {
throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
}
long exp = e;
BigInteger result = BigInteger.ONE;
BigInteger k2p = k;
while (exp != 0) {
if ((exp & 0x1) != 0) {
result = result.multiply(k2p);
}
k2p = k2p.multiply(k2p);
exp >>= 1;
}
return result;
}
/**
* Raise a BigInteger to a BigInteger power.
*
* @param k Number to raise.
* @param e Exponent (must be positive or zero).
* @return k<sup>e</sup>
* @throws IllegalArgumentException if {@code e < 0}.
*/
public static BigInteger pow(final BigInteger k, final BigInteger e) {
if (e.compareTo(BigInteger.ZERO) < 0) {
throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2);
}
BigInteger exp = e;
BigInteger result = BigInteger.ONE;
BigInteger k2p = k;
while (!BigInteger.ZERO.equals(exp)) {
if (exp.testBit(0)) {
result = result.multiply(k2p);
}
k2p = k2p.multiply(k2p);
exp = exp.shiftRight(1);
}
return result;
}
/**
* Returns true if the argument is a power of two.
*
* @param n the number to test
* @return true if the argument is a power of two
*/
public static boolean isPowerOfTwo(long n) {
return n > 0 && (n & (n - 1)) == 0;
}
/**
* Returns the unsigned remainder from dividing the first argument
* by the second where each argument and the result is interpreted
* as an unsigned value.
*
* <p>Implementation note
*
* <p>In v1.0 this method did not use the {@code long} datatype.
* Modern 64-bit processors make use of the {@code long} datatype
* faster than an algorithm using the {@code int} datatype. This method
* now delegates to {@link Integer#remainderUnsigned(int, int)}
* which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
*
* @param dividend the value to be divided
* @param divisor the value doing the dividing
* @return the unsigned remainder of the first argument divided by
* the second argument.
* @see Integer#remainderUnsigned(int, int)
*/
public static int remainderUnsigned(int dividend, int divisor) {
return Integer.remainderUnsigned(dividend, divisor);
}
/**
* Returns the unsigned remainder from dividing the first argument
* by the second where each argument and the result is interpreted
* as an unsigned value.
*
* <p>Implementation note
*
* <p>This method does not use the {@code BigInteger} datatype.
* The JDK implementation of {@link Long#remainderUnsigned(long, long)}
* uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
* From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
* even faster due to use of an intrinsic method.
*
* @param dividend the value to be divided
* @param divisor the value doing the dividing
* @return the unsigned remainder of the first argument divided by
* the second argument.
* @see Long#remainderUnsigned(long, long)
*/
public static long remainderUnsigned(long dividend, long divisor) {
// Adapts the divideUnsigned method to compute the remainder.
if (divisor < 0) {
// Using unsigned compare:
// if dividend < divisor: return dividend
// else: return dividend - divisor
// Subtracting divisor using masking is more complex in this case
// and we use a condition
return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor;
}
// From Hacker's Delight 2.0, section 9.3
final long q = ((dividend >>> 1) / divisor) << 1;
final long r = dividend - q * divisor;
// unsigned r: 0 <= r < 2 * divisor
// if (r < divisor): r
// else: r - divisor
// The compare of unsigned r can be done using:
// return (r + Long.MIN_VALUE) < (divisor | Long.MIN_VALUE) ? r : r - divisor
// Here we subtract divisor if (r - divisor) is positive, else the result is r.
// This can be done by flipping the sign bit and
// creating a mask as -1 or 0 by signed shift.
return r - (divisor & (~(r - divisor) >> 63));
}
/**
* Returns the unsigned quotient of dividing the first argument by
* the second where each argument and the result is interpreted as
* an unsigned value.
* <p>Note that in two's complement arithmetic, the three other
* basic arithmetic operations of add, subtract, and multiply are
* bit-wise identical if the two operands are regarded as both
* being signed or both being unsigned. Therefore separate {@code
* addUnsigned}, etc. methods are not provided.</p>
*
* <p>Implementation note
*
* <p>In v1.0 this method did not use the {@code long} datatype.
* Modern 64-bit processors make use of the {@code long} datatype
* faster than an algorithm using the {@code int} datatype. This method
* now delegates to {@link Integer#divideUnsigned(int, int)}
* which uses {@code long} arithmetic; or from JDK 19 an intrinsic method.
*
* @param dividend the value to be divided
* @param divisor the value doing the dividing
* @return the unsigned quotient of the first argument divided by
* the second argument
* @see Integer#divideUnsigned(int, int)
*/
public static int divideUnsigned(int dividend, int divisor) {
return Integer.divideUnsigned(dividend, divisor);
}
/**
* Returns the unsigned quotient of dividing the first argument by
* the second where each argument and the result is interpreted as
* an unsigned value.
* <p>Note that in two's complement arithmetic, the three other
* basic arithmetic operations of add, subtract, and multiply are
* bit-wise identical if the two operands are regarded as both
* being signed or both being unsigned. Therefore separate {@code
* addUnsigned}, etc. methods are not provided.</p>
*
* <p>Implementation note
*
* <p>This method does not use the {@code BigInteger} datatype.
* The JDK implementation of {@link Long#divideUnsigned(long, long)}
* uses {@code BigInteger} prior to JDK 17 and this method is 15-25x faster.
* From JDK 17 onwards the JDK implementation is as fast; or from JDK 19
* even faster due to use of an intrinsic method.
*
* @param dividend the value to be divided
* @param divisor the value doing the dividing
* @return the unsigned quotient of the first argument divided by
* the second argument.
* @see Long#divideUnsigned(long, long)
*/
public static long divideUnsigned(long dividend, long divisor) {
// The implementation is a Java port of algorithm described in the book
// "Hacker's Delight 2.0" (section 9.3 "Unsigned short division from signed division").
// Adapts 6-line predicate expressions program with (u >=) an unsigned compare
// using the provided branchless variants.
if (divisor < 0) {
// line 1 branchless:
// q <- (dividend (u >=) divisor)
return (dividend & ~(dividend - divisor)) >>> 63;
}
final long q = ((dividend >>> 1) / divisor) << 1;
final long r = dividend - q * divisor;
// line 5 branchless:
// q <- q + (r (u >=) divisor)
return q + ((r | ~(r - divisor)) >>> 63);
}
/**
* Exception.
*/
private static class NumbersArithmeticException extends ArithmeticException {
/** Serializable version Id. */
private static final long serialVersionUID = 20180130L;
/**
* Create an exception where the message is constructed by applying
* {@link String#format(String, Object...)}.
*
* @param message Exception message format string
* @param args Arguments for formatting the message
*/
NumbersArithmeticException(String message, Object... args) {
super(String.format(message, args));
}
}
}