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