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 * <p>Special cases:</p> 270 * <ul> 271 * <li>{@code k^0} returns {@code 1} (including {@code k=0}) 272 * <li>{@code k^1} returns {@code k} (including {@code k=0}) 273 * <li>{@code 0^0} returns {@code 1} 274 * <li>{@code 0^e} returns {@code 0} 275 * <li>{@code 1^e} returns {@code 1} 276 * <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even 277 * </ul> 278 * 279 * @param k Number to raise. 280 * @param e Exponent (must be positive or zero). 281 * @return \( k^e \) 282 * @throws IllegalArgumentException if {@code e < 0}. 283 * @throws ArithmeticException if the result would overflow. 284 */ 285 public static int pow(final int k, 286 final int e) { 287 if (e < 0) { 288 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 289 } 290 291 if (k == 0) { 292 return e == 0 ? 1 : 0; 293 } 294 295 if (k == 1) { 296 return 1; 297 } 298 299 if (k == -1) { 300 return (e & 1) == 0 ? 1 : -1; 301 } 302 303 if (e >= 31) { 304 throw new ArithmeticException("integer overflow"); 305 } 306 307 int exp = e; 308 int result = 1; 309 int k2p = k; 310 while (true) { 311 if ((exp & 0x1) != 0) { 312 result = Math.multiplyExact(result, k2p); 313 } 314 315 exp >>= 1; 316 if (exp == 0) { 317 break; 318 } 319 320 k2p = Math.multiplyExact(k2p, k2p); 321 } 322 323 return result; 324 } 325 326 /** 327 * Raise a long to an int power. 328 * 329 * <p>Special cases:</p> 330 * <ul> 331 * <li>{@code k^0} returns {@code 1} (including {@code k=0}) 332 * <li>{@code k^1} returns {@code k} (including {@code k=0}) 333 * <li>{@code 0^0} returns {@code 1} 334 * <li>{@code 0^e} returns {@code 0} 335 * <li>{@code 1^e} returns {@code 1} 336 * <li>{@code (-1)^e} returns {@code -1 or 1} if {@code e} is odd or even 337 * </ul> 338 * 339 * @param k Number to raise. 340 * @param e Exponent (must be positive or zero). 341 * @return \( k^e \) 342 * @throws IllegalArgumentException if {@code e < 0}. 343 * @throws ArithmeticException if the result would overflow. 344 */ 345 public static long pow(final long k, 346 final int e) { 347 if (e < 0) { 348 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 349 } 350 351 if (k == 0L) { 352 return e == 0 ? 1L : 0L; 353 } 354 355 if (k == 1L) { 356 return 1L; 357 } 358 359 if (k == -1L) { 360 return (e & 1) == 0 ? 1L : -1L; 361 } 362 363 if (e >= 63) { 364 throw new ArithmeticException("long overflow"); 365 } 366 367 int exp = e; 368 long result = 1; 369 long k2p = k; 370 while (true) { 371 if ((exp & 0x1) != 0) { 372 result = Math.multiplyExact(result, k2p); 373 } 374 375 exp >>= 1; 376 if (exp == 0) { 377 break; 378 } 379 380 k2p = Math.multiplyExact(k2p, k2p); 381 } 382 383 return result; 384 } 385 386 /** 387 * Raise a BigInteger to an int power. 388 * 389 * @param k Number to raise. 390 * @param e Exponent (must be positive or zero). 391 * @return k<sup>e</sup> 392 * @throws IllegalArgumentException if {@code e < 0}. 393 */ 394 public static BigInteger pow(final BigInteger k, int e) { 395 if (e < 0) { 396 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 397 } 398 399 return k.pow(e); 400 } 401 402 /** 403 * Raise a BigInteger to a long power. 404 * 405 * @param k Number to raise. 406 * @param e Exponent (must be positive or zero). 407 * @return k<sup>e</sup> 408 * @throws IllegalArgumentException if {@code e < 0}. 409 */ 410 public static BigInteger pow(final BigInteger k, final long e) { 411 if (e < 0) { 412 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 413 } 414 415 long exp = e; 416 BigInteger result = BigInteger.ONE; 417 BigInteger k2p = k; 418 while (exp != 0) { 419 if ((exp & 0x1) != 0) { 420 result = result.multiply(k2p); 421 } 422 k2p = k2p.multiply(k2p); 423 exp >>= 1; 424 } 425 426 return result; 427 428 } 429 430 /** 431 * Raise a BigInteger to a BigInteger power. 432 * 433 * @param k Number to raise. 434 * @param e Exponent (must be positive or zero). 435 * @return k<sup>e</sup> 436 * @throws IllegalArgumentException if {@code e < 0}. 437 */ 438 public static BigInteger pow(final BigInteger k, final BigInteger e) { 439 if (e.compareTo(BigInteger.ZERO) < 0) { 440 throw new IllegalArgumentException(NEGATIVE_EXPONENT_1 + e + NEGATIVE_EXPONENT_2); 441 } 442 443 BigInteger exp = e; 444 BigInteger result = BigInteger.ONE; 445 BigInteger k2p = k; 446 while (!BigInteger.ZERO.equals(exp)) { 447 if (exp.testBit(0)) { 448 result = result.multiply(k2p); 449 } 450 k2p = k2p.multiply(k2p); 451 exp = exp.shiftRight(1); 452 } 453 454 return result; 455 } 456 457 /** 458 * Returns true if the argument is a power of two. 459 * 460 * @param n the number to test 461 * @return true if the argument is a power of two 462 */ 463 public static boolean isPowerOfTwo(long n) { 464 return (n > 0) && ((n & (n - 1)) == 0); 465 } 466 467 /** 468 * Returns the unsigned remainder from dividing the first argument 469 * by the second where each argument and the result is interpreted 470 * as an unsigned value. 471 * <p>This method does not use the {@code long} datatype.</p> 472 * 473 * @param dividend the value to be divided 474 * @param divisor the value doing the dividing 475 * @return the unsigned remainder of the first argument divided by 476 * the second argument. 477 */ 478 public static int remainderUnsigned(int dividend, int divisor) { 479 if (divisor >= 0) { 480 if (dividend >= 0) { 481 return dividend % divisor; 482 } 483 // The implementation is a Java port of algorithm described in the book 484 // "Hacker's Delight" (section "Unsigned short division from signed division"). 485 final int q = ((dividend >>> 1) / divisor) << 1; 486 dividend -= q * divisor; 487 if (dividend < 0 || dividend >= divisor) { 488 dividend -= divisor; 489 } 490 return dividend; 491 } 492 return dividend >= 0 || dividend < divisor ? dividend : dividend - divisor; 493 } 494 495 /** 496 * Returns the unsigned remainder from dividing the first argument 497 * by the second where each argument and the result is interpreted 498 * as an unsigned value. 499 * <p>This method does not use the {@code BigInteger} datatype.</p> 500 * 501 * @param dividend the value to be divided 502 * @param divisor the value doing the dividing 503 * @return the unsigned remainder of the first argument divided by 504 * the second argument. 505 */ 506 public static long remainderUnsigned(long dividend, long divisor) { 507 if (divisor >= 0L) { 508 if (dividend >= 0L) { 509 return dividend % divisor; 510 } 511 // The implementation is a Java port of algorithm described in the book 512 // "Hacker's Delight" (section "Unsigned short division from signed division"). 513 final long q = ((dividend >>> 1) / divisor) << 1; 514 dividend -= q * divisor; 515 if (dividend < 0L || dividend >= divisor) { 516 dividend -= divisor; 517 } 518 return dividend; 519 } 520 return dividend >= 0L || dividend < divisor ? dividend : dividend - divisor; 521 } 522 523 /** 524 * Returns the unsigned quotient of dividing the first argument by 525 * the second where each argument and the result is interpreted as 526 * an unsigned value. 527 * <p>Note that in two's complement arithmetic, the three other 528 * basic arithmetic operations of add, subtract, and multiply are 529 * bit-wise identical if the two operands are regarded as both 530 * being signed or both being unsigned. Therefore separate {@code 531 * addUnsigned}, etc. methods are not provided.</p> 532 * <p>This method does not use the {@code long} datatype.</p> 533 * 534 * @param dividend the value to be divided 535 * @param divisor the value doing the dividing 536 * @return the unsigned quotient of the first argument divided by 537 * the second argument 538 */ 539 public static int divideUnsigned(int dividend, int divisor) { 540 if (divisor >= 0) { 541 if (dividend >= 0) { 542 return dividend / divisor; 543 } 544 // The implementation is a Java port of algorithm described in the book 545 // "Hacker's Delight" (section "Unsigned short division from signed division"). 546 int q = ((dividend >>> 1) / divisor) << 1; 547 dividend -= q * divisor; 548 if (dividend < 0L || dividend >= divisor) { 549 q++; 550 } 551 return q; 552 } 553 return dividend >= 0 || dividend < divisor ? 0 : 1; 554 } 555 556 /** 557 * Returns the unsigned quotient of dividing the first argument by 558 * the second where each argument and the result is interpreted as 559 * an unsigned value. 560 * <p>Note that in two's complement arithmetic, the three other 561 * basic arithmetic operations of add, subtract, and multiply are 562 * bit-wise identical if the two operands are regarded as both 563 * being signed or both being unsigned. Therefore separate {@code 564 * addUnsigned}, etc. methods are not provided.</p> 565 * <p>This method does not use the {@code BigInteger} datatype.</p> 566 * 567 * @param dividend the value to be divided 568 * @param divisor the value doing the dividing 569 * @return the unsigned quotient of the first argument divided by 570 * the second argument. 571 */ 572 public static long divideUnsigned(long dividend, long divisor) { 573 if (divisor >= 0L) { 574 if (dividend >= 0L) { 575 return dividend / divisor; 576 } 577 // The implementation is a Java port of algorithm described in the book 578 // "Hacker's Delight" (section "Unsigned short division from signed division"). 579 long q = ((dividend >>> 1) / divisor) << 1; 580 dividend -= q * divisor; 581 if (dividend < 0L || dividend >= divisor) { 582 q++; 583 } 584 return q; 585 } 586 return dividend >= 0L || dividend < divisor ? 0L : 1L; 587 } 588 589 /** 590 * Exception. 591 */ 592 private static class NumbersArithmeticException extends ArithmeticException { 593 /** Serializable version Id. */ 594 private static final long serialVersionUID = 20180130L; 595 596 /** 597 * Constructor with a specific message. 598 * 599 * @param message Message pattern providing the specific context of 600 * the error. 601 * @param args Arguments. 602 */ 603 NumbersArithmeticException(String message, Object... args) { 604 super(MessageFormat.format(message, args)); 605 } 606 } 607}