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
21 /**
22 * Some useful, arithmetics related, additions to the built-in functions in
23 * {@link Math}.
24 *
25 */
26 public final class ArithmeticUtils {
27
28 /** Negative exponent exception message part 1. */
29 private static final String NEGATIVE_EXPONENT_1 = "negative exponent ({";
30 /** Negative exponent exception message part 2. */
31 private static final String NEGATIVE_EXPONENT_2 = "})";
32
33 /** Private constructor. */
34 private ArithmeticUtils() {
35 // intentionally empty.
36 }
37
38 /**
39 * Computes the greatest common divisor of the absolute value of two
40 * numbers, using a modified version of the "binary gcd" method.
41 * See Knuth 4.5.2 algorithm B.
42 * The algorithm is due to Josef Stein (1961).
43 * <br>
44 * Special cases:
45 * <ul>
46 * <li>The invocations
47 * {@code gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)},
48 * {@code gcd(Integer.MIN_VALUE, 0)} and
49 * {@code gcd(0, Integer.MIN_VALUE)} throw an
50 * {@code ArithmeticException}, because the result would be 2^31, which
51 * is too large for an int value.</li>
52 * <li>The result of {@code gcd(x, x)}, {@code gcd(0, x)} and
53 * {@code gcd(x, 0)} is the absolute value of {@code x}, except
54 * for the special cases above.</li>
55 * <li>The invocation {@code gcd(0, 0)} is the only one which returns
56 * {@code 0}.</li>
57 * </ul>
58 *
59 * <p>Two numbers are relatively prime, or coprime, if their gcd is 1.</p>
60 *
61 * @param p Number.
62 * @param q Number.
63 * @return the greatest common divisor (never negative).
64 * @throws ArithmeticException if the result cannot be represented as
65 * a non-negative {@code int} value.
66 */
67 public static int gcd(int p, int q) {
68 // Perform the gcd algorithm on negative numbers, so that -2^31 does not
69 // need to be handled separately
70 int a = p > 0 ? -p : p;
71 int b = q > 0 ? -q : q;
72
73 int negatedGcd;
74 if (a == 0) {
75 negatedGcd = b;
76 } else if (b == 0) {
77 negatedGcd = a;
78 } else {
79 // Make "a" and "b" odd, keeping track of common power of 2.
80 final int aTwos = Integer.numberOfTrailingZeros(a);
81 final int bTwos = Integer.numberOfTrailingZeros(b);
82 a >>= aTwos;
83 b >>= bTwos;
84 final int shift = Math.min(aTwos, bTwos);
85
86 // "a" and "b" are negative and odd.
87 // If a < b then "gdc(a, b)" is equal to "gcd(a - b, b)".
88 // If a > b then "gcd(a, b)" is equal to "gcd(b - a, a)".
89 // Hence, in the successive iterations:
90 // "a" becomes the negative absolute difference of the current values,
91 // "b" becomes that value of the two that is closer to zero.
92 while (a != b) {
93 final int delta = a - b;
94 b = Math.max(a, b);
95 a = delta > 0 ? -delta : delta;
96
97 // Remove any power of 2 in "a" ("b" is guaranteed to be odd).
98 a >>= Integer.numberOfTrailingZeros(a);
99 }
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 }