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