1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.numbers.core; 18 19 import java.math.BigInteger; 20 import java.text.MessageFormat; 21 22 /** 23 * Some useful, arithmetics related, additions to the built-in functions in 24 * {@link Math}. 25 * 26 */ 27 public final class ArithmeticUtils { 28 29 /** Overflow gcd exception message for 2^63. */ 30 private static final String OVERFLOW_GCD_MESSAGE_2_POWER_63 = "overflow: gcd({0}, {1}) is 2^63"; 31 32 /** Negative exponent exception message part 1. */ 33 private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({"; 34 /** Negative exponent exception message part 2. */ 35 private static final String NEGATIVE_EXPONENT_2 = "})"; 36 37 /** Private constructor. */ 38 private ArithmeticUtils() { 39 // intentionally empty. 40 } 41 42 /** 43 * Computes the greatest common divisor of the absolute value of two 44 * numbers, using a modified version of the "binary gcd" method. 45 * See Knuth 4.5.2 algorithm B. 46 * The algorithm is due to Josef Stein (1961). 47 * <br> 48 * Special cases: 49 * <ul> 50 * <li>The invocations 51 * {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)}, 52 * {@code gcd(Integer.MIN_VALUE, 0)} and 53 * {@code gcd(0, Integer.MIN_VALUE)} throw an 54 * {@code ArithmeticException}, because the result would be 2^31, which 55 * is too large for an int value.</li> 56 * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and 57 * {@code gcd(x, 0)} is the absolute value of {@code x}, except 58 * for the special cases above.</li> 59 * <li>The invocation {@code gcd(0, 0)} is the only one which returns 60 * {@code 0}.</li> 61 * </ul> 62 * 63 * @param p Number. 64 * @param q Number. 65 * @return the greatest common divisor (never negative). 66 * @throws ArithmeticException if the result cannot be represented as 67 * a non-negative {@code int} value. 68 */ 69 public static int gcd(int p, int q) { 70 // Perform the gcd algorithm on negative numbers, so that -2^31 does not 71 // need to be handled separately 72 int a = p > 0 ? -p : p; 73 int b = q > 0 ? -q : q; 74 75 int negatedGcd; 76 if (a == 0) { 77 negatedGcd = b; 78 } else if (b == 0) { 79 negatedGcd = a; 80 } else { 81 // Make "a" and "b" odd, keeping track of common power of 2. 82 final int aTwos = Integer.numberOfTrailingZeros(a); 83 final int bTwos = Integer.numberOfTrailingZeros(b); 84 a >>= aTwos; 85 b >>= bTwos; 86 final int shift = Math.min(aTwos, bTwos); 87 88 // "a" and "b" are negative and odd. 89 // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)". 90 // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)". 91 // Hence, in the successive iterations: 92 // "a" becomes the negative absolute difference of the current values, 93 // "b" becomes that value of the two that is closer to zero. 94 while (a != b) { 95 final int delta = a - b; 96 b = Math.max(a, b); 97 a = delta > 0 ? -delta : delta; 98 99 // 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 } 110 return -negatedGcd; 111 } 112 113 /** 114 * <p> 115 * Gets the greatest common divisor of the absolute value of two numbers, 116 * using the "binary gcd" method which avoids division and modulo 117 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef 118 * Stein (1961). 119 * </p> 120 * Special cases: 121 * <ul> 122 * <li>The invocations 123 * {@code gcd(Long.MIN_VALUE, Long.MIN_VALUE)}, 124 * {@code gcd(Long.MIN_VALUE, 0L)} and 125 * {@code gcd(0L, Long.MIN_VALUE)} throw an 126 * {@code ArithmeticException}, because the result would be 2^63, which 127 * is too large for a long value.</li> 128 * <li>The result of {@code gcd(x, x)}, {@code gcd(0L, x)} and 129 * {@code gcd(x, 0L)} is the absolute value of {@code x}, except 130 * for the special cases above. 131 * <li>The invocation {@code gcd(0L, 0L)} is the only one which returns 132 * {@code 0L}.</li> 133 * </ul> 134 * 135 * @param p Number. 136 * @param q Number. 137 * @return the greatest common divisor, never negative. 138 * @throws ArithmeticException if the result cannot be represented as 139 * a non-negative {@code long} value. 140 */ 141 public static long gcd(final long p, final long q) { 142 long u = p; 143 long v = q; 144 if (u == 0 || v == 0) { 145 if (u == Long.MIN_VALUE || v == Long.MIN_VALUE) { 146 throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63, 147 p, q); 148 } 149 return Math.abs(u) + Math.abs(v); 150 } 151 // keep u and v negative, as negative integers range down to 152 // -2^63, while positive numbers can only be as large as 2^63-1 153 // (i.e. we can't necessarily negate a negative number without 154 // overflow) 155 /* assert u!=0 && v!=0; */ 156 if (u > 0) { 157 u = -u; 158 } // make u negative 159 if (v > 0) { 160 v = -v; 161 } // make v negative 162 // B1. [Find power of 2] 163 int k = 0; 164 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are 165 // both even... 166 u /= 2; 167 v /= 2; 168 k++; // cast out twos. 169 } 170 if (k == 63) { 171 throw new NumbersArithmeticException(OVERFLOW_GCD_MESSAGE_2_POWER_63, 172 p, q); 173 } 174 // B2. Initialize: u and v have been divided by 2^k and at least 175 // one is odd. 176 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */; 177 // t negative: u was odd, v may be even (t replaces v) 178 // t positive: u was even, v is odd (t replaces u) 179 do { 180 /* assert u<0 && v<0; */ 181 // B4/B3: cast out twos from t. 182 while ((t & 1) == 0) { // while t is even.. 183 t /= 2; // cast out twos 184 } 185 // B5 [reset max(u,v)] 186 if (t > 0) { 187 u = -t; 188 } else { 189 v = t; 190 } 191 // B6/B3. at this point both u and v should be odd. 192 t = (v - u) / 2; 193 // |u| larger: t positive (replace u) 194 // |v| larger: t negative (replace v) 195 } while (t != 0); 196 return -u * (1L << k); // gcd is u*2^k 197 } 198 199 /** 200 * <p> 201 * Returns the least common multiple of the absolute value of two numbers, 202 * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}. 203 * </p> 204 * Special cases: 205 * <ul> 206 * <li>The invocations {@code lcm(Integer.MIN_VALUE, n)} and 207 * {@code lcm(n, Integer.MIN_VALUE)}, where {@code abs(n)} is a 208 * power of 2, throw an {@code ArithmeticException}, because the result 209 * would be 2^31, which is too large for an int value.</li> 210 * <li>The result of {@code lcm(0, x)} and {@code lcm(x, 0)} is 211 * {@code 0} for any {@code x}. 212 * </ul> 213 * 214 * @param a Number. 215 * @param b Number. 216 * @return the least common multiple, never negative. 217 * @throws ArithmeticException if the result cannot be represented as 218 * a non-negative {@code int} value. 219 */ 220 public static int lcm(int a, int b) { 221 if (a == 0 || b == 0) { 222 return 0; 223 } 224 final int lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b)); 225 if (lcm == Integer.MIN_VALUE) { 226 throw new NumbersArithmeticException("overflow: lcm({0}, {1}) is 2^31", 227 a, b); 228 } 229 return lcm; 230 } 231 232 /** 233 * <p> 234 * Returns the least common multiple of the absolute value of two numbers, 235 * using the formula {@code lcm(a,b) = (a / gcd(a,b)) * b}. 236 * </p> 237 * Special cases: 238 * <ul> 239 * <li>The invocations {@code lcm(Long.MIN_VALUE, n)} and 240 * {@code lcm(n, Long.MIN_VALUE)}, where {@code abs(n)} is a 241 * power of 2, throw an {@code ArithmeticException}, because the result 242 * would be 2^63, which is too large for an int value.</li> 243 * <li>The result of {@code lcm(0L, x)} and {@code lcm(x, 0L)} is 244 * {@code 0L} for any {@code x}. 245 * </ul> 246 * 247 * @param a Number. 248 * @param b Number. 249 * @return the least common multiple, never negative. 250 * @throws ArithmeticException if the result cannot be represented 251 * as a non-negative {@code long} value. 252 */ 253 public static long lcm(long a, long b) { 254 if (a == 0 || b == 0) { 255 return 0; 256 } 257 final long lcm = Math.abs(Math.multiplyExact(a / gcd(a, b), b)); 258 if (lcm == Long.MIN_VALUE) { 259 throw new NumbersArithmeticException("overflow: lcm({0}, {1}) is 2^63", 260 a, b); 261 } 262 return lcm; 263 } 264 265 /** 266 * Raise an int to an int power. 267 * 268 * <p>Special cases:</p> 269 * <ul> 270 * <li>{@code k^0} returns {@code 1} (including {@code k=0}) 271 * <li>{@code k^1} returns {@code k} (including {@code k=0}) 272 * <li>{@code 0^0} returns {@code 1} 273 * <li>{@code 0^e} returns {@code 0} 274 * <li>{@code 1^e} returns {@code 1} 275 * <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even 276 * </ul> 277 * 278 * @param k Number to raise. 279 * @param e Exponent (must be positive or zero). 280 * @return \( k^e \) 281 * @throws IllegalArgumentException if {@code e < 0}. 282 * @throws ArithmeticException if the result would overflow. 283 */ 284 public static int pow(final int k, 285 final int e) { 286 if (e < 0) { 287 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 288 } 289 290 if (k == 0) { 291 return e == 0 ? 1 : 0; 292 } 293 294 if (k == 1) { 295 return 1; 296 } 297 298 if (k == -1) { 299 return (e & 1) == 0 ? 1 : -1; 300 } 301 302 if (e >= 31) { 303 throw new ArithmeticException("integer overflow"); 304 } 305 306 int exp = e; 307 int result = 1; 308 int k2p = k; 309 while (true) { 310 if ((exp & 0x1) != 0) { 311 result = Math.multiplyExact(result, k2p); 312 } 313 314 exp >>= 1; 315 if (exp == 0) { 316 break; 317 } 318 319 k2p = Math.multiplyExact(k2p, k2p); 320 } 321 322 return result; 323 } 324 325 /** 326 * Raise a long to an int power. 327 * 328 * <p>Special cases:</p> 329 * <ul> 330 * <li>{@code k^0} returns {@code 1} (including {@code k=0}) 331 * <li>{@code k^1} returns {@code k} (including {@code k=0}) 332 * <li>{@code 0^0} returns {@code 1} 333 * <li>{@code 0^e} returns {@code 0} 334 * <li>{@code 1^e} returns {@code 1} 335 * <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even 336 * </ul> 337 * 338 * @param k Number to raise. 339 * @param e Exponent (must be positive or zero). 340 * @return \( k^e \) 341 * @throws IllegalArgumentException if {@code e < 0}. 342 * @throws ArithmeticException if the result would overflow. 343 */ 344 public static long pow(final long k, 345 final int e) { 346 if (e < 0) { 347 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 348 } 349 350 if (k == 0L) { 351 return e == 0 ? 1L : 0L; 352 } 353 354 if (k == 1L) { 355 return 1L; 356 } 357 358 if (k == -1L) { 359 return (e & 1) == 0 ? 1L : -1L; 360 } 361 362 if (e >= 63) { 363 throw new ArithmeticException("long overflow"); 364 } 365 366 int exp = e; 367 long result = 1; 368 long k2p = k; 369 while (true) { 370 if ((exp & 0x1) != 0) { 371 result = Math.multiplyExact(result, k2p); 372 } 373 374 exp >>= 1; 375 if (exp == 0) { 376 break; 377 } 378 379 k2p = Math.multiplyExact(k2p, k2p); 380 } 381 382 return result; 383 } 384 385 /** 386 * Raise a BigInteger to an int power. 387 * 388 * @param k Number to raise. 389 * @param e Exponent (must be positive or zero). 390 * @return k<sup>e</sup> 391 * @throws IllegalArgumentException if {@code e < 0}. 392 */ 393 public static BigInteger pow(final BigInteger k, int e) { 394 if (e < 0) { 395 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 396 } 397 398 return k.pow(e); 399 } 400 401 /** 402 * Raise a BigInteger to a long power. 403 * 404 * @param k Number to raise. 405 * @param e Exponent (must be positive or zero). 406 * @return k<sup>e</sup> 407 * @throws IllegalArgumentException if {@code e < 0}. 408 */ 409 public static BigInteger pow(final BigInteger k, final long e) { 410 if (e < 0) { 411 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 412 } 413 414 long exp = e; 415 BigInteger result = BigInteger.ONE; 416 BigInteger k2p = k; 417 while (exp != 0) { 418 if ((exp & 0x1) != 0) { 419 result = result.multiply(k2p); 420 } 421 k2p = k2p.multiply(k2p); 422 exp >>= 1; 423 } 424 425 return result; 426 427 } 428 429 /** 430 * Raise a BigInteger to a BigInteger power. 431 * 432 * @param k Number to raise. 433 * @param e Exponent (must be positive or zero). 434 * @return k<sup>e</sup> 435 * @throws IllegalArgumentException if {@code e < 0}. 436 */ 437 public static BigInteger pow(final BigInteger k, final BigInteger e) { 438 if (e.compareTo(BigInteger.ZERO) < 0) { 439 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 440 } 441 442 BigInteger exp = e; 443 BigInteger result = BigInteger.ONE; 444 BigInteger k2p = k; 445 while (!BigInteger.ZERO.equals(exp)) { 446 if (exp.testBit(0)) { 447 result = result.multiply(k2p); 448 } 449 k2p = k2p.multiply(k2p); 450 exp = exp.shiftRight(1); 451 } 452 453 return result; 454 } 455 456 /** 457 * Returns true if the argument is a power of two. 458 * 459 * @param n the number to test 460 * @return true if the argument is a power of two 461 */ 462 public static boolean isPowerOfTwo(long n) { 463 return n > 0 && (n & (n - 1)) == 0; 464 } 465 466 /** 467 * Returns the unsigned remainder from dividing the first argument 468 * by the second where each argument and the result is interpreted 469 * as an unsigned value. 470 * <p>This method does not use the {@code long} datatype.</p> 471 * 472 * @param dividend the value to be divided 473 * @param divisor the value doing the dividing 474 * @return the unsigned remainder of the first argument divided by 475 * the second argument. 476 */ 477 public static int remainderUnsigned(int dividend, int divisor) { 478 if (divisor >= 0) { 479 if (dividend >= 0) { 480 return dividend % divisor; 481 } 482 // The implementation is a Java port of algorithm described in the book 483 // "Hacker's Delight" (section "Unsigned short division from signed division"). 484 final int q = ((dividend >>> 1) / divisor) << 1; 485 dividend -= q * divisor; 486 if (dividend < 0 || dividend >= divisor) { 487 dividend -= divisor; 488 } 489 return dividend; 490 } 491 return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor; 492 } 493 494 /** 495 * Returns the unsigned remainder from dividing the first argument 496 * by the second where each argument and the result is interpreted 497 * as an unsigned value. 498 * <p>This method does not use the {@code BigInteger} datatype.</p> 499 * 500 * @param dividend the value to be divided 501 * @param divisor the value doing the dividing 502 * @return the unsigned remainder of the first argument divided by 503 * the second argument. 504 */ 505 public static long remainderUnsigned(long dividend, long divisor) { 506 if (divisor >= 0L) { 507 if (dividend >= 0L) { 508 return dividend % divisor; 509 } 510 // The implementation is a Java port of algorithm described in the book 511 // "Hacker's Delight" (section "Unsigned short division from signed division"). 512 final long q = ((dividend >>> 1) / divisor) << 1; 513 dividend -= q * divisor; 514 if (dividend < 0L || dividend >= divisor) { 515 dividend -= divisor; 516 } 517 return dividend; 518 } 519 return dividend >= 0L || dividend < divisor ? dividend : dividend - divisor; 520 } 521 522 /** 523 * Returns the unsigned quotient of dividing the first argument by 524 * the second where each argument and the result is interpreted as 525 * an unsigned value. 526 * <p>Note that in two's complement arithmetic, the three other 527 * basic arithmetic operations of add, subtract, and multiply are 528 * bit-wise identical if the two operands are regarded as both 529 * being signed or both being unsigned. Therefore separate {@code 530 * addUnsigned}, etc. methods are not provided.</p> 531 * <p>This method does not use the {@code long} datatype.</p> 532 * 533 * @param dividend the value to be divided 534 * @param divisor the value doing the dividing 535 * @return the unsigned quotient of the first argument divided by 536 * the second argument 537 */ 538 public static int divideUnsigned(int dividend, int divisor) { 539 if (divisor >= 0) { 540 if (dividend >= 0) { 541 return dividend / divisor; 542 } 543 // The implementation is a Java port of algorithm described in the book 544 // "Hacker's Delight" (section "Unsigned short division from signed division"). 545 int q = ((dividend >>> 1) / divisor) << 1; 546 dividend -= q * divisor; 547 if (dividend < 0L || dividend >= divisor) { 548 q++; 549 } 550 return q; 551 } 552 return dividend >= 0 || dividend < divisor ? 0 : 1; 553 } 554 555 /** 556 * Returns the unsigned quotient of dividing the first argument by 557 * the second where each argument and the result is interpreted as 558 * an unsigned value. 559 * <p>Note that in two's complement arithmetic, the three other 560 * basic arithmetic operations of add, subtract, and multiply are 561 * bit-wise identical if the two operands are regarded as both 562 * being signed or both being unsigned. Therefore separate {@code 563 * addUnsigned}, etc. methods are not provided.</p> 564 * <p>This method does not use the {@code BigInteger} datatype.</p> 565 * 566 * @param dividend the value to be divided 567 * @param divisor the value doing the dividing 568 * @return the unsigned quotient of the first argument divided by 569 * the second argument. 570 */ 571 public static long divideUnsigned(long dividend, long divisor) { 572 if (divisor >= 0L) { 573 if (dividend >= 0L) { 574 return dividend / divisor; 575 } 576 // The implementation is a Java port of algorithm described in the book 577 // "Hacker's Delight" (section "Unsigned short division from signed division"). 578 long q = ((dividend >>> 1) / divisor) << 1; 579 dividend -= q * divisor; 580 if (dividend < 0L || dividend >= divisor) { 581 q++; 582 } 583 return q; 584 } 585 return dividend >= 0L || dividend < divisor ? 0L : 1L; 586 } 587 588 /** 589 * Exception. 590 */ 591 private static class NumbersArithmeticException extends ArithmeticException { 592 /** Serializable version Id. */ 593 private static final long serialVersionUID = 20180130L; 594 595 /** 596 * Constructor with a specific message. 597 * 598 * @param message Message pattern providing the specific context of 599 * the error. 600 * @param args Arguments. 601 */ 602 NumbersArithmeticException(String message, Object... args) { 603 super(MessageFormat.format(message, args)); 604 } 605 } 606 }