1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.util;
18
19 import java.math.BigDecimal;
20
21 /**
22 * Some useful additions to the built-in functions in {@link Math}.
23 * @version $Revision: 321510 $ $Date: 2005-10-15 15:33:14 -0700 (Sat, 15 Oct 2005) $
24 */
25 public final class MathUtils {
26
27 /** -1.0 cast as a byte. */
28 private static final byte NB = (byte)-1;
29
30 /** -1.0 cast as a short. */
31 private static final short NS = (short)-1;
32
33 /** 1.0 cast as a byte. */
34 private static final byte PB = (byte)1;
35
36 /** 1.0 cast as a short. */
37 private static final short PS = (short)1;
38
39 /** 0.0 cast as a byte. */
40 private static final byte ZB = (byte)0;
41
42 /** 0.0 cast as a short. */
43 private static final short ZS = (short)0;
44
45 /**
46 * Private Constructor
47 */
48 private MathUtils() {
49 super();
50 }
51
52 /**
53 * Add two integers, checking for overflow.
54 *
55 * @param x an addend
56 * @param y an addend
57 * @return the sum <code>x+y</code>
58 * @throws ArithmeticException if the result can not be represented as an
59 * int
60 * @since 1.1
61 */
62 public static int addAndCheck(int x, int y) {
63 long s = (long)x + (long)y;
64 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
65 throw new ArithmeticException("overflow: add");
66 }
67 return (int)s;
68 }
69
70 /**
71 * Returns an exact representation of the <a
72 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
73 * Coefficient</a>, "<code>n choose k</code>", the number of
74 * <code>k</code>-element subsets that can be selected from an
75 * <code>n</code>-element set.
76 * <p>
77 * <Strong>Preconditions</strong>:
78 * <ul>
79 * <li> <code>0 <= k <= n </code> (otherwise
80 * <code>IllegalArgumentException</code> is thrown)</li>
81 * <li> The result is small enough to fit into a <code>long</code>. The
82 * largest value of <code>n</code> for which all coefficients are
83 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
84 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException
85 * </code> is
86 * thrown.</li>
87 * </ul>
88 *
89 * @param n the size of the set
90 * @param k the size of the subsets to be counted
91 * @return <code>n choose k</code>
92 * @throws IllegalArgumentException if preconditions are not met.
93 * @throws ArithmeticException if the result is too large to be represented
94 * by a long integer.
95 */
96 public static long binomialCoefficient(final int n, final int k) {
97 if (n < k) {
98 throw new IllegalArgumentException(
99 "must have n >= k for binomial coefficient (n,k)");
100 }
101 if (n < 0) {
102 throw new IllegalArgumentException(
103 "must have n >= 0 for binomial coefficient (n,k)");
104 }
105 if ((n == k) || (k == 0)) {
106 return 1;
107 }
108 if ((k == 1) || (k == n - 1)) {
109 return n;
110 }
111
112 long result = Math.round(binomialCoefficientDouble(n, k));
113 if (result == Long.MAX_VALUE) {
114 throw new ArithmeticException(
115 "result too large to represent in a long integer");
116 }
117 return result;
118 }
119
120 /**
121 * Returns a <code>double</code> representation of the <a
122 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
123 * Coefficient</a>, "<code>n choose k</code>", the number of
124 * <code>k</code>-element subsets that can be selected from an
125 * <code>n</code>-element set.
126 * <p>
127 * <Strong>Preconditions</strong>:
128 * <ul>
129 * <li> <code>0 <= k <= n </code> (otherwise
130 * <code>IllegalArgumentException</code> is thrown)</li>
131 * <li> The result is small enough to fit into a <code>double</code>. The
132 * largest value of <code>n</code> for which all coefficients are <
133 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
134 * Double.POSITIVE_INFINITY is returned</li>
135 * </ul>
136 *
137 * @param n the size of the set
138 * @param k the size of the subsets to be counted
139 * @return <code>n choose k</code>
140 * @throws IllegalArgumentException if preconditions are not met.
141 */
142 public static double binomialCoefficientDouble(final int n, final int k) {
143 return Math.floor(Math.exp(binomialCoefficientLog(n, k)) + 0.5);
144 }
145
146 /**
147 * Returns the natural <code>log</code> of the <a
148 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
149 * Coefficient</a>, "<code>n choose k</code>", the number of
150 * <code>k</code>-element subsets that can be selected from an
151 * <code>n</code>-element set.
152 * <p>
153 * <Strong>Preconditions</strong>:
154 * <ul>
155 * <li> <code>0 <= k <= n </code> (otherwise
156 * <code>IllegalArgumentException</code> is thrown)</li>
157 * </ul>
158 *
159 * @param n the size of the set
160 * @param k the size of the subsets to be counted
161 * @return <code>n choose k</code>
162 * @throws IllegalArgumentException if preconditions are not met.
163 */
164 public static double binomialCoefficientLog(final int n, final int k) {
165 if (n < k) {
166 throw new IllegalArgumentException(
167 "must have n >= k for binomial coefficient (n,k)");
168 }
169 if (n < 0) {
170 throw new IllegalArgumentException(
171 "must have n >= 0 for binomial coefficient (n,k)");
172 }
173 if ((n == k) || (k == 0)) {
174 return 0;
175 }
176 if ((k == 1) || (k == n - 1)) {
177 return Math.log((double)n);
178 }
179 double logSum = 0;
180
181
182 for (int i = k + 1; i <= n; i++) {
183 logSum += Math.log((double)i);
184 }
185
186
187 for (int i = 2; i <= n - k; i++) {
188 logSum -= Math.log((double)i);
189 }
190
191 return logSum;
192 }
193
194 /**
195 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
196 * hyperbolic cosine</a> of x.
197 *
198 * @param x double value for which to find the hyperbolic cosine
199 * @return hyperbolic cosine of x
200 */
201 public static double cosh(double x) {
202 return (Math.exp(x) + Math.exp(-x)) / 2.0;
203 }
204
205 /**
206 * Returns true iff both arguments are NaN or neither is NaN and they are
207 * equal
208 *
209 * @param x first value
210 * @param y second value
211 * @return true if the values are equal or both are NaN
212 */
213 public static boolean equals(double x, double y) {
214 return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
215 }
216
217 /**
218 * Returns n!. Shorthand for <code>n</code> <a
219 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
220 * product of the numbers <code>1,...,n</code>.
221 * <p>
222 * <Strong>Preconditions</strong>:
223 * <ul>
224 * <li> <code>n >= 0</code> (otherwise
225 * <code>IllegalArgumentException</code> is thrown)</li>
226 * <li> The result is small enough to fit into a <code>long</code>. The
227 * largest value of <code>n</code> for which <code>n!</code> <
228 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
229 * an <code>ArithMeticException </code> is thrown.</li>
230 * </ul>
231 * </p>
232 *
233 * @param n argument
234 * @return <code>n!</code>
235 * @throws ArithmeticException if the result is too large to be represented
236 * by a long integer.
237 * @throws IllegalArgumentException if n < 0
238 */
239 public static long factorial(final int n) {
240 long result = Math.round(factorialDouble(n));
241 if (result == Long.MAX_VALUE) {
242 throw new ArithmeticException(
243 "result too large to represent in a long integer");
244 }
245 return result;
246 }
247
248 /**
249 * Returns n!. Shorthand for <code>n</code> <a
250 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
251 * product of the numbers <code>1,...,n</code> as a <code>double</code>.
252 * <p>
253 * <Strong>Preconditions</strong>:
254 * <ul>
255 * <li> <code>n >= 0</code> (otherwise
256 * <code>IllegalArgumentException</code> is thrown)</li>
257 * <li> The result is small enough to fit into a <code>double</code>. The
258 * largest value of <code>n</code> for which <code>n!</code> <
259 * Double.MAX_VALUE</code> is 170. If the computed value exceeds
260 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
261 * </ul>
262 * </p>
263 *
264 * @param n argument
265 * @return <code>n!</code>
266 * @throws IllegalArgumentException if n < 0
267 */
268 public static double factorialDouble(final int n) {
269 if (n < 0) {
270 throw new IllegalArgumentException("must have n >= 0 for n!");
271 }
272 return Math.floor(Math.exp(factorialLog(n)) + 0.5);
273 }
274
275 /**
276 * Returns the natural logarithm of n!.
277 * <p>
278 * <Strong>Preconditions</strong>:
279 * <ul>
280 * <li> <code>n >= 0</code> (otherwise
281 * <code>IllegalArgumentException</code> is thrown)</li>
282 * </ul>
283 *
284 * @param n argument
285 * @return <code>n!</code>
286 * @throws IllegalArgumentException if preconditions are not met.
287 */
288 public static double factorialLog(final int n) {
289 if (n < 0) {
290 throw new IllegalArgumentException("must have n > 0 for n!");
291 }
292 double logSum = 0;
293 for (int i = 2; i <= n; i++) {
294 logSum += Math.log((double)i);
295 }
296 return logSum;
297 }
298
299 /**
300 * <p>
301 * Gets the greatest common divisor of the absolute value of two numbers,
302 * using the "binary gcd" method which avoids division and modulo
303 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
304 * Stein (1961).
305 * </p>
306 *
307 * @param u a non-zero number
308 * @param v a non-zero number
309 * @return the greatest common divisor, never zero
310 * @since 1.1
311 */
312 public static int gcd(int u, int v) {
313 if (u * v == 0) {
314 return (Math.abs(u) + Math.abs(v));
315 }
316
317
318
319
320
321 if (u > 0) {
322 u = -u;
323 }
324 if (v > 0) {
325 v = -v;
326 }
327
328 int k = 0;
329 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) {
330
331 u /= 2;
332 v /= 2;
333 k++;
334 }
335 if (k == 31) {
336 throw new ArithmeticException("overflow: gcd is 2^31");
337 }
338
339
340 int t = ((u & 1) == 1) ? v : -(u / 2)
341
342
343 do {
344
345
346 while ((t & 1) == 0) {
347 t /= 2;
348 }
349
350 if (t > 0) {
351 u = -t;
352 } else {
353 v = t;
354 }
355
356 t = (v - u) / 2;
357
358
359 } while (t != 0);
360 return -u * (1 << k);
361 }
362
363 /**
364 * Returns an integer hash code representing the given double value.
365 *
366 * @param value the value to be hashed
367 * @return the hash code
368 */
369 public static int hash(double value) {
370 long bits = Double.doubleToLongBits(value);
371 return (int)(bits ^ (bits >>> 32));
372 }
373
374 /**
375 * For a byte value x, this method returns (byte)(+1) if x >= 0 and
376 * (byte)(-1) if x < 0.
377 *
378 * @param x the value, a byte
379 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
380 */
381 public static byte indicator(final byte x) {
382 return (x >= ZB) ? PB : NB;
383 }
384
385 /**
386 * For a double precision value x, this method returns +1.0 if x >= 0 and
387 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
388 * <code>NaN</code>.
389 *
390 * @param x the value, a double
391 * @return +1.0 or -1.0, depending on the sign of x
392 */
393 public static double indicator(final double x) {
394 if (Double.isNaN(x)) {
395 return Double.NaN;
396 }
397 return (x >= 0.0) ? 1.0 : -1.0;
398 }
399
400 /**
401 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
402 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
403 *
404 * @param x the value, a float
405 * @return +1.0F or -1.0F, depending on the sign of x
406 */
407 public static float indicator(final float x) {
408 if (Float.isNaN(x)) {
409 return Float.NaN;
410 }
411 return (x >= 0.0F) ? 1.0F : -1.0F;
412 }
413
414 /**
415 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
416 *
417 * @param x the value, an int
418 * @return +1 or -1, depending on the sign of x
419 */
420 public static int indicator(final int x) {
421 return (x >= 0) ? 1 : -1;
422 }
423
424 /**
425 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
426 *
427 * @param x the value, a long
428 * @return +1L or -1L, depending on the sign of x
429 */
430 public static long indicator(final long x) {
431 return (x >= 0L) ? 1L : -1L;
432 }
433
434 /**
435 * For a short value x, this method returns (short)(+1) if x >= 0 and
436 * (short)(-1) if x < 0.
437 *
438 * @param x the value, a short
439 * @return (short)(+1) or (short)(-1), depending on the sign of x
440 */
441 public static short indicator(final short x) {
442 return (x >= ZS) ? PS : NS;
443 }
444
445 /**
446 * Returns the least common multiple between two integer values.
447 *
448 * @param a the first integer value.
449 * @param b the second integer value.
450 * @return the least common multiple between a and b.
451 * @throws ArithmeticException if the lcm is too large to store as an int
452 * @since 1.1
453 */
454 public static int lcm(int a, int b) {
455 return Math.abs(mulAndCheck(a / gcd(a, b), b));
456 }
457
458 /**
459 * Multiply two integers, checking for overflow.
460 *
461 * @param x a factor
462 * @param y a factor
463 * @return the product <code>x*y</code>
464 * @throws ArithmeticException if the result can not be represented as an
465 * int
466 * @since 1.1
467 */
468 public static int mulAndCheck(int x, int y) {
469 long m = ((long)x) * ((long)y);
470 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
471 throw new ArithmeticException("overflow: mul");
472 }
473 return (int)m;
474 }
475
476 /**
477 * Round the given value to the specified number of decimal places. The
478 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
479 *
480 * @param x the value to round.
481 * @param scale the number of digits to the right of the decimal point.
482 * @return the rounded value.
483 * @since 1.1
484 */
485 public static double round(double x, int scale) {
486 return round(x, scale, BigDecimal.ROUND_HALF_UP);
487 }
488
489 /**
490 * Round the given value to the specified number of decimal places. The
491 * value is rounded using the given method which is any method defined in
492 * {@link BigDecimal}.
493 *
494 * @param x the value to round.
495 * @param scale the number of digits to the right of the decimal point.
496 * @param roundingMethod the rounding method as defined in
497 * {@link BigDecimal}.
498 * @return the rounded value.
499 * @since 1.1
500 */
501 public static double round(double x, int scale, int roundingMethod) {
502 double sign = indicator(x);
503 double factor = Math.pow(10.0, scale) * sign;
504 return roundUnscaled(x * factor, sign, roundingMethod) / factor;
505 }
506
507 /**
508 * Round the given value to the specified number of decimal places. The
509 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
510 *
511 * @param x the value to round.
512 * @param scale the number of digits to the right of the decimal point.
513 * @return the rounded value.
514 * @since 1.1
515 */
516 public static float round(float x, int scale) {
517 return round(x, scale, BigDecimal.ROUND_HALF_UP);
518 }
519
520 /**
521 * Round the given value to the specified number of decimal places. The
522 * value is rounded using the given method which is any method defined in
523 * {@link BigDecimal}.
524 *
525 * @param x the value to round.
526 * @param scale the number of digits to the right of the decimal point.
527 * @param roundingMethod the rounding method as defined in
528 * {@link BigDecimal}.
529 * @return the rounded value.
530 * @since 1.1
531 */
532 public static float round(float x, int scale, int roundingMethod) {
533 float sign = indicator(x);
534 float factor = (float)Math.pow(10.0f, scale) * sign;
535 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
536 }
537
538 /**
539 * Round the given non-negative, value to the "nearest" integer. Nearest is
540 * determined by the rounding method specified. Rounding methods are defined
541 * in {@link BigDecimal}.
542 *
543 * @param unscaled the value to round.
544 * @param sign the sign of the original, scaled value.
545 * @param roundingMethod the rounding method as defined in
546 * {@link BigDecimal}.
547 * @return the rounded value.
548 * @since 1.1
549 */
550 private static double roundUnscaled(double unscaled, double sign,
551 int roundingMethod) {
552 switch (roundingMethod) {
553 case BigDecimal.ROUND_CEILING :
554 if (sign == -1) {
555 unscaled = Math.floor(unscaled);
556 } else {
557 unscaled = Math.ceil(unscaled);
558 }
559 break;
560 case BigDecimal.ROUND_DOWN :
561 unscaled = Math.floor(unscaled);
562 break;
563 case BigDecimal.ROUND_FLOOR :
564 if (sign == -1) {
565 unscaled = Math.ceil(unscaled);
566 } else {
567 unscaled = Math.floor(unscaled);
568 }
569 break;
570 case BigDecimal.ROUND_HALF_DOWN : {
571 double fraction = Math.abs(unscaled - Math.floor(unscaled));
572 if (fraction > 0.5) {
573 unscaled = Math.ceil(unscaled);
574 } else {
575 unscaled = Math.floor(unscaled);
576 }
577 break;
578 }
579 case BigDecimal.ROUND_HALF_EVEN : {
580 double fraction = Math.abs(unscaled - Math.floor(unscaled));
581 if (fraction > 0.5) {
582 unscaled = Math.ceil(unscaled);
583 } else if (fraction < 0.5) {
584 unscaled = Math.floor(unscaled);
585 } else {
586 if (Math.floor(unscaled) / 2.0 == Math.floor(Math
587 .floor(unscaled) / 2.0)) {
588 unscaled = Math.floor(unscaled);
589 } else {
590 unscaled = Math.ceil(unscaled);
591 }
592 }
593 break;
594 }
595 case BigDecimal.ROUND_HALF_UP : {
596 double fraction = Math.abs(unscaled - Math.floor(unscaled));
597 if (fraction >= 0.5) {
598 unscaled = Math.ceil(unscaled);
599 } else {
600 unscaled = Math.floor(unscaled);
601 }
602 break;
603 }
604 case BigDecimal.ROUND_UNNECESSARY :
605 if (unscaled != Math.floor(unscaled)) {
606 throw new ArithmeticException("Inexact result from rounding");
607 }
608 break;
609 case BigDecimal.ROUND_UP :
610 unscaled = Math.ceil(unscaled);
611 break;
612 default :
613 throw new IllegalArgumentException("Invalid rounding method.");
614 }
615 return unscaled;
616 }
617
618 /**
619 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
620 * for byte value <code>x</code>.
621 * <p>
622 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
623 * x = 0, and (byte)(-1) if x < 0.
624 *
625 * @param x the value, a byte
626 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
627 */
628 public static byte sign(final byte x) {
629 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
630 }
631
632 /**
633 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
634 * for double precision <code>x</code>.
635 * <p>
636 * For a double value <code>x</code>, this method returns
637 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
638 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
639 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
640 *
641 * @param x the value, a double
642 * @return +1.0, 0.0, or -1.0, depending on the sign of x
643 */
644 public static double sign(final double x) {
645 if (Double.isNaN(x)) {
646 return Double.NaN;
647 }
648 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
649 }
650
651 /**
652 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
653 * for float value <code>x</code>.
654 * <p>
655 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
656 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
657 * is <code>NaN</code>.
658 *
659 * @param x the value, a float
660 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
661 */
662 public static float sign(final float x) {
663 if (Float.isNaN(x)) {
664 return Float.NaN;
665 }
666 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
667 }
668
669 /**
670 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
671 * for int value <code>x</code>.
672 * <p>
673 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
674 * if x < 0.
675 *
676 * @param x the value, an int
677 * @return +1, 0, or -1, depending on the sign of x
678 */
679 public static int sign(final int x) {
680 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
681 }
682
683 /**
684 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
685 * for long value <code>x</code>.
686 * <p>
687 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
688 * -1L if x < 0.
689 *
690 * @param x the value, a long
691 * @return +1L, 0L, or -1L, depending on the sign of x
692 */
693 public static long sign(final long x) {
694 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
695 }
696
697 /**
698 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
699 * for short value <code>x</code>.
700 * <p>
701 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
702 * if x = 0, and (short)(-1) if x < 0.
703 *
704 * @param x the value, a short
705 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
706 * x
707 */
708 public static short sign(final short x) {
709 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
710 }
711
712 /**
713 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
714 * hyperbolic sine</a> of x.
715 *
716 * @param x double value for which to find the hyperbolic sine
717 * @return hyperbolic sine of x
718 */
719 public static double sinh(double x) {
720 return (Math.exp(x) - Math.exp(-x)) / 2.0;
721 }
722
723 /**
724 * Subtract two integers, checking for overflow.
725 *
726 * @param x the minuend
727 * @param y the subtrahend
728 * @return the difference <code>x-y</code>
729 * @throws ArithmeticException if the result can not be represented as an
730 * int
731 * @since 1.1
732 */
733 public static int subAndCheck(int x, int y) {
734 long s = (long)x - (long)y;
735 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
736 throw new ArithmeticException("overflow: subtract");
737 }
738 return (int)s;
739 }
740 }