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;
import java.text.MessageFormat;
/**
* Some useful, arithmetics related, additions to the built-in functions in
* {@link Math}.
*
*/
public final class ArithmeticUtils {
/** Overflow gcd exception message for 2^63. */
private static final String OVERFLOW_GCD_MESSAGE_2_POWER_63 = "overflow: gcd({0}, {1}) is 2^63";
/** 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>
*
* @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({0}, {1}) 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>
*
* @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(final long p, final long q) {
long u = p;
long v = q;
if (u == 0 || v == 0) {
if (u == Long.MIN_VALUE || v == Long.MIN_VALUE) {
throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63,
p, q);
}
return Math.abs(u) + Math.abs(v);
}
// keep u and v negative, as negative integers range down to
// -2^63, while positive numbers can only be as large as 2^63-1
// (i.e. we can't necessarily negate a negative number without
// overflow)
/* assert u!=0 && v!=0; */
if (u > 0) {
u = -u;
} // make u negative
if (v > 0) {
v = -v;
} // make v negative
// B1. [Find power of 2]
int k = 0;
while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
// both even...
u /= 2;
v /= 2;
k++; // cast out twos.
}
if (k == 63) {
throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63,
p, q);
}
// B2. Initialize: u and v have been divided by 2^k and at least
// one is odd.
long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
// t negative: u was odd, v may be even (t replaces v)
// t positive: u was even, v is odd (t replaces u)
do {
/* assert u<0 && v<0; */
// B4/B3: cast out twos from t.
while ((t & 1) == 0) { // while t is even..
t /= 2; // cast out twos
}
// B5 [reset max(u,v)]
if (t > 0) {
u = -t;
} else {
v = t;
}
// B6/B3. at this point both u and v should be odd.
t = (v - u) / 2;
// |u| larger: t positive (replace u)
// |v| larger: t negative (replace v)
} while (t != 0);
return -u * (1L << k); // gcd is u*2^k
}
/**
* <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({0}, {1}) 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({0}, {1}) 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>This method does not use the {@code long} datatype.</p>
*
* @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.
*/
public static int remainderUnsigned(int dividend, int divisor) {
if (divisor >= 0) {
if (dividend >= 0) {
return dividend % divisor;
}
// The implementation is a Java port of algorithm described in the book
// "Hacker's Delight" (section "Unsigned short division from signed division").
final int q = ((dividend >>> 1) / divisor) << 1;
dividend -= q * divisor;
if (dividend < 0 || dividend >= divisor) {
dividend -= divisor;
}
return dividend;
}
return dividend >= 0 || dividend < divisor ? dividend : 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>This method does not use the {@code BigInteger} datatype.</p>
*
* @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.
*/
public static long remainderUnsigned(long dividend, long divisor) {
if (divisor >= 0L) {
if (dividend >= 0L) {
return dividend % divisor;
}
// The implementation is a Java port of algorithm described in the book
// "Hacker's Delight" (section "Unsigned short division from signed division").
final long q = ((dividend >>> 1) / divisor) << 1;
dividend -= q * divisor;
if (dividend < 0L || dividend >= divisor) {
dividend -= divisor;
}
return dividend;
}
return dividend >= 0L || dividend < divisor ? dividend : 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>This method does not use the {@code long} datatype.</p>
*
* @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
*/
public static int divideUnsigned(int dividend, int divisor) {
if (divisor >= 0) {
if (dividend >= 0) {
return dividend / divisor;
}
// The implementation is a Java port of algorithm described in the book
// "Hacker's Delight" (section "Unsigned short division from signed division").
int q = ((dividend >>> 1) / divisor) << 1;
dividend -= q * divisor;
if (dividend < 0L || dividend >= divisor) {
q++;
}
return q;
}
return dividend >= 0 || dividend < divisor ? 0 : 1;
}
/**
* 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>This method does not use the {@code BigInteger} datatype.</p>
*
* @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.
*/
public static long divideUnsigned(long dividend, long divisor) {
if (divisor >= 0L) {
if (dividend >= 0L) {
return dividend / divisor;
}
// The implementation is a Java port of algorithm described in the book
// "Hacker's Delight" (section "Unsigned short division from signed division").
long q = ((dividend >>> 1) / divisor) << 1;
dividend -= q * divisor;
if (dividend < 0L || dividend >= divisor) {
q++;
}
return q;
}
return dividend >= 0L || dividend < divisor ? 0L : 1L;
}
/**
* Exception.
*/
private static class NumbersArithmeticException extends ArithmeticException {
/** Serializable version Id. */
private static final long serialVersionUID = 20180130L;
/**
* Constructor with a specific message.
*
* @param message Message pattern providing the specific context of
* the error.
* @param args Arguments.
*/
NumbersArithmeticException(String message, Object... args) {
super(MessageFormat.format(message, args));
}
}
}