commons-numbers-core/src/main/java/org/apache/commons/numbers/core/ArithmeticUtils.java - commons-numbers - Git at Google

 /*
  * 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);
         } else {
             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));
         }
     }
 }
	/*
	* 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);
	} else {
	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));
	}
	}
	}