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 */
017
018 package org.apache.commons.math.util;
019
020 import java.math.BigDecimal;
021 import java.math.BigInteger;
022 import java.util.Arrays;
023
024 import org.apache.commons.math.MathRuntimeException;
025
026 /**
027 * Some useful additions to the built-in functions in {@link Math}.
028 * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
029 */
030 public final class MathUtils {
031
032 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
033 public static final double EPSILON = 0x1.0p-53;
034
035 /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
036 * <p>In IEEE 754 arithmetic, this is also the smallest normalized
037 * number 2<sup>-1022</sup>.</p>
038 */
039 public static final double SAFE_MIN = 0x1.0p-1022;
040
041 /** -1.0 cast as a byte. */
042 private static final byte NB = (byte)-1;
043
044 /** -1.0 cast as a short. */
045 private static final short NS = (short)-1;
046
047 /** 1.0 cast as a byte. */
048 private static final byte PB = (byte)1;
049
050 /** 1.0 cast as a short. */
051 private static final short PS = (short)1;
052
053 /** 0.0 cast as a byte. */
054 private static final byte ZB = (byte)0;
055
056 /** 0.0 cast as a short. */
057 private static final short ZS = (short)0;
058
059 /** 2 π. */
060 private static final double TWO_PI = 2 * Math.PI;
061
062 /** Gap between NaN and regular numbers. */
063 private static final int NAN_GAP = 4 * 1024 * 1024;
064
065 /** Offset to order signed double numbers lexicographically. */
066 private static final long SGN_MASK = 0x8000000000000000L;
067
068 /**
069 * Private Constructor
070 */
071 private MathUtils() {
072 super();
073 }
074
075 /**
076 * Add two integers, checking for overflow.
077 *
078 * @param x an addend
079 * @param y an addend
080 * @return the sum <code>x+y</code>
081 * @throws ArithmeticException if the result can not be represented as an
082 * int
083 * @since 1.1
084 */
085 public static int addAndCheck(int x, int y) {
086 long s = (long)x + (long)y;
087 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
088 throw new ArithmeticException("overflow: add");
089 }
090 return (int)s;
091 }
092
093 /**
094 * Add two long integers, checking for overflow.
095 *
096 * @param a an addend
097 * @param b an addend
098 * @return the sum <code>a+b</code>
099 * @throws ArithmeticException if the result can not be represented as an
100 * long
101 * @since 1.2
102 */
103 public static long addAndCheck(long a, long b) {
104 return addAndCheck(a, b, "overflow: add");
105 }
106
107 /**
108 * Add two long integers, checking for overflow.
109 *
110 * @param a an addend
111 * @param b an addend
112 * @param msg the message to use for any thrown exception.
113 * @return the sum <code>a+b</code>
114 * @throws ArithmeticException if the result can not be represented as an
115 * long
116 * @since 1.2
117 */
118 private static long addAndCheck(long a, long b, String msg) {
119 long ret;
120 if (a > b) {
121 // use symmetry to reduce boundary cases
122 ret = addAndCheck(b, a, msg);
123 } else {
124 // assert a <= b
125
126 if (a < 0) {
127 if (b < 0) {
128 // check for negative overflow
129 if (Long.MIN_VALUE - b <= a) {
130 ret = a + b;
131 } else {
132 throw new ArithmeticException(msg);
133 }
134 } else {
135 // opposite sign addition is always safe
136 ret = a + b;
137 }
138 } else {
139 // assert a >= 0
140 // assert b >= 0
141
142 // check for positive overflow
143 if (a <= Long.MAX_VALUE - b) {
144 ret = a + b;
145 } else {
146 throw new ArithmeticException(msg);
147 }
148 }
149 }
150 return ret;
151 }
152
153 /**
154 * Returns an exact representation of the <a
155 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
156 * Coefficient</a>, "<code>n choose k</code>", the number of
157 * <code>k</code>-element subsets that can be selected from an
158 * <code>n</code>-element set.
159 * <p>
160 * <Strong>Preconditions</strong>:
161 * <ul>
162 * <li> <code>0 <= k <= n </code> (otherwise
163 * <code>IllegalArgumentException</code> is thrown)</li>
164 * <li> The result is small enough to fit into a <code>long</code>. The
165 * largest value of <code>n</code> for which all coefficients are
166 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
167 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is
168 * thrown.</li>
169 * </ul></p>
170 *
171 * @param n the size of the set
172 * @param k the size of the subsets to be counted
173 * @return <code>n choose k</code>
174 * @throws IllegalArgumentException if preconditions are not met.
175 * @throws ArithmeticException if the result is too large to be represented
176 * by a long integer.
177 */
178 public static long binomialCoefficient(final int n, final int k) {
179 checkBinomial(n, k);
180 if ((n == k) || (k == 0)) {
181 return 1;
182 }
183 if ((k == 1) || (k == n - 1)) {
184 return n;
185 }
186 // Use symmetry for large k
187 if (k > n / 2)
188 return binomialCoefficient(n, n - k);
189
190 // We use the formula
191 // (n choose k) = n! / (n-k)! / k!
192 // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
193 // which could be written
194 // (n choose k) == (n-1 choose k-1) * n / k
195 long result = 1;
196 if (n <= 61) {
197 // For n <= 61, the naive implementation cannot overflow.
198 for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
199 result = result * i / j;
200 }
201 } else if (n <= 66) {
202 // For n > 61 but n <= 66, the result cannot overflow,
203 // but we must take care not to overflow intermediate values.
204 for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
205 // We know that (result * i) is divisible by j,
206 // but (result * i) may overflow, so we split j:
207 // Filter out the gcd, d, so j/d and i/d are integer.
208 // result is divisible by (j/d) because (j/d)
209 // is relative prime to (i/d) and is a divisor of
210 // result * (i/d).
211 long d = gcd(i, j);
212 result = (result / (j / d)) * (i / d);
213 }
214 } else {
215 // For n > 66, a result overflow might occur, so we check
216 // the multiplication, taking care to not overflow
217 // unnecessary.
218 for (int j = 1, i = n - k + 1; j <= k; i++, j++) {
219 long d = gcd(i, j);
220 result = mulAndCheck((result / (j / d)), (i / d));
221 }
222 }
223 return result;
224 }
225
226 /**
227 * Returns a <code>double</code> representation of the <a
228 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
229 * Coefficient</a>, "<code>n choose k</code>", the number of
230 * <code>k</code>-element subsets that can be selected from an
231 * <code>n</code>-element set.
232 * <p>
233 * <Strong>Preconditions</strong>:
234 * <ul>
235 * <li> <code>0 <= k <= n </code> (otherwise
236 * <code>IllegalArgumentException</code> is thrown)</li>
237 * <li> The result is small enough to fit into a <code>double</code>. The
238 * largest value of <code>n</code> for which all coefficients are <
239 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
240 * Double.POSITIVE_INFINITY is returned</li>
241 * </ul></p>
242 *
243 * @param n the size of the set
244 * @param k the size of the subsets to be counted
245 * @return <code>n choose k</code>
246 * @throws IllegalArgumentException if preconditions are not met.
247 */
248 public static double binomialCoefficientDouble(final int n, final int k) {
249 checkBinomial(n, k);
250 if ((n == k) || (k == 0)) {
251 return 1d;
252 }
253 if ((k == 1) || (k == n - 1)) {
254 return n;
255 }
256 if (k > n/2) {
257 return binomialCoefficientDouble(n, n - k);
258 }
259 if (n < 67) {
260 return binomialCoefficient(n,k);
261 }
262
263 double result = 1d;
264 for (int i = 1; i <= k; i++) {
265 result *= (double)(n - k + i) / (double)i;
266 }
267
268 return Math.floor(result + 0.5);
269 }
270
271 /**
272 * Returns the natural <code>log</code> of the <a
273 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
274 * Coefficient</a>, "<code>n choose k</code>", the number of
275 * <code>k</code>-element subsets that can be selected from an
276 * <code>n</code>-element set.
277 * <p>
278 * <Strong>Preconditions</strong>:
279 * <ul>
280 * <li> <code>0 <= k <= n </code> (otherwise
281 * <code>IllegalArgumentException</code> is thrown)</li>
282 * </ul></p>
283 *
284 * @param n the size of the set
285 * @param k the size of the subsets to be counted
286 * @return <code>n choose k</code>
287 * @throws IllegalArgumentException if preconditions are not met.
288 */
289 public static double binomialCoefficientLog(final int n, final int k) {
290 checkBinomial(n, k);
291 if ((n == k) || (k == 0)) {
292 return 0;
293 }
294 if ((k == 1) || (k == n - 1)) {
295 return Math.log(n);
296 }
297
298 /*
299 * For values small enough to do exact integer computation,
300 * return the log of the exact value
301 */
302 if (n < 67) {
303 return Math.log(binomialCoefficient(n,k));
304 }
305
306 /*
307 * Return the log of binomialCoefficientDouble for values that will not
308 * overflow binomialCoefficientDouble
309 */
310 if (n < 1030) {
311 return Math.log(binomialCoefficientDouble(n, k));
312 }
313
314 if (k > n / 2) {
315 return binomialCoefficientLog(n, n - k);
316 }
317
318 /*
319 * Sum logs for values that could overflow
320 */
321 double logSum = 0;
322
323 // n!/(n-k)!
324 for (int i = n - k + 1; i <= n; i++) {
325 logSum += Math.log(i);
326 }
327
328 // divide by k!
329 for (int i = 2; i <= k; i++) {
330 logSum -= Math.log(i);
331 }
332
333 return logSum;
334 }
335
336 /**
337 * Check binomial preconditions.
338 * @param n the size of the set
339 * @param k the size of the subsets to be counted
340 * @exception IllegalArgumentException if preconditions are not met.
341 */
342 private static void checkBinomial(final int n, final int k)
343 throws IllegalArgumentException {
344 if (n < k) {
345 throw MathRuntimeException.createIllegalArgumentException(
346 "must have n >= k for binomial coefficient (n,k), got n = {0}, k = {1}",
347 n, k);
348 }
349 if (n < 0) {
350 throw MathRuntimeException.createIllegalArgumentException(
351 "must have n >= 0 for binomial coefficient (n,k), got n = {0}",
352 n);
353 }
354 }
355
356 /**
357 * Compares two numbers given some amount of allowed error.
358 *
359 * @param x the first number
360 * @param y the second number
361 * @param eps the amount of error to allow when checking for equality
362 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li>
363 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li>
364 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul>
365 */
366 public static int compareTo(double x, double y, double eps) {
367 if (equals(x, y, eps)) {
368 return 0;
369 } else if (x < y) {
370 return -1;
371 }
372 return 1;
373 }
374
375 /**
376 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
377 * hyperbolic cosine</a> of x.
378 *
379 * @param x double value for which to find the hyperbolic cosine
380 * @return hyperbolic cosine of x
381 */
382 public static double cosh(double x) {
383 return (Math.exp(x) + Math.exp(-x)) / 2.0;
384 }
385
386 /**
387 * Returns true iff both arguments are NaN or neither is NaN and they are
388 * equal
389 *
390 * @param x first value
391 * @param y second value
392 * @return true if the values are equal or both are NaN
393 */
394 public static boolean equals(double x, double y) {
395 return ((Double.isNaN(x) && Double.isNaN(y)) || x == y);
396 }
397
398 /**
399 * Returns true iff both arguments are equal or within the range of allowed
400 * error (inclusive).
401 * <p>
402 * Two NaNs are considered equals, as are two infinities with same sign.
403 * </p>
404 *
405 * @param x first value
406 * @param y second value
407 * @param eps the amount of absolute error to allow
408 * @return true if the values are equal or within range of each other
409 */
410 public static boolean equals(double x, double y, double eps) {
411 return equals(x, y) || (Math.abs(y - x) <= eps);
412 }
413
414 /**
415 * Returns true iff both arguments are equal or within the range of allowed
416 * error (inclusive).
417 * Adapted from <a
418 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
419 * Bruce Dawson</a>
420 *
421 * @param x first value
422 * @param y second value
423 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
424 * values between {@code x} and {@code y}.
425 * @return {@code true} if there are less than {@code maxUlps} floating
426 * point values between {@code x} and {@code y}
427 */
428 public static boolean equals(double x, double y, int maxUlps) {
429 // Check that "maxUlps" is non-negative and small enough so that the
430 // default NAN won't compare as equal to anything.
431 assert maxUlps > 0 && maxUlps < NAN_GAP;
432
433 long xInt = Double.doubleToLongBits(x);
434 long yInt = Double.doubleToLongBits(y);
435
436 // Make lexicographically ordered as a two's-complement integer.
437 if (xInt < 0) {
438 xInt = SGN_MASK - xInt;
439 }
440 if (yInt < 0) {
441 yInt = SGN_MASK - yInt;
442 }
443
444 return Math.abs(xInt - yInt) <= maxUlps;
445 }
446
447 /**
448 * Returns true iff both arguments are null or have same dimensions
449 * and all their elements are {@link #equals(double,double) equals}
450 *
451 * @param x first array
452 * @param y second array
453 * @return true if the values are both null or have same dimension
454 * and equal elements
455 * @since 1.2
456 */
457 public static boolean equals(double[] x, double[] y) {
458 if ((x == null) || (y == null)) {
459 return !((x == null) ^ (y == null));
460 }
461 if (x.length != y.length) {
462 return false;
463 }
464 for (int i = 0; i < x.length; ++i) {
465 if (!equals(x[i], y[i])) {
466 return false;
467 }
468 }
469 return true;
470 }
471
472 /** All long-representable factorials */
473 private static final long[] factorials = new long[]
474 {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
475 479001600, 6227020800l, 87178291200l, 1307674368000l, 20922789888000l,
476 355687428096000l, 6402373705728000l, 121645100408832000l,
477 2432902008176640000l};
478
479 /**
480 * Returns n!. Shorthand for <code>n</code> <a
481 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
482 * product of the numbers <code>1,...,n</code>.
483 * <p>
484 * <Strong>Preconditions</strong>:
485 * <ul>
486 * <li> <code>n >= 0</code> (otherwise
487 * <code>IllegalArgumentException</code> is thrown)</li>
488 * <li> The result is small enough to fit into a <code>long</code>. The
489 * largest value of <code>n</code> for which <code>n!</code> <
490 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
491 * an <code>ArithMeticException </code> is thrown.</li>
492 * </ul>
493 * </p>
494 *
495 * @param n argument
496 * @return <code>n!</code>
497 * @throws ArithmeticException if the result is too large to be represented
498 * by a long integer.
499 * @throws IllegalArgumentException if n < 0
500 */
501 public static long factorial(final int n) {
502 if (n < 0) {
503 throw MathRuntimeException.createIllegalArgumentException(
504 "must have n >= 0 for n!, got n = {0}",
505 n);
506 }
507 if (n > 20) {
508 throw new ArithmeticException(
509 "factorial value is too large to fit in a long");
510 }
511 return factorials[n];
512 }
513
514 /**
515 * Returns n!. Shorthand for <code>n</code> <a
516 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
517 * product of the numbers <code>1,...,n</code> as a <code>double</code>.
518 * <p>
519 * <Strong>Preconditions</strong>:
520 * <ul>
521 * <li> <code>n >= 0</code> (otherwise
522 * <code>IllegalArgumentException</code> is thrown)</li>
523 * <li> The result is small enough to fit into a <code>double</code>. The
524 * largest value of <code>n</code> for which <code>n!</code> <
525 * Double.MAX_VALUE</code> is 170. If the computed value exceeds
526 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
527 * </ul>
528 * </p>
529 *
530 * @param n argument
531 * @return <code>n!</code>
532 * @throws IllegalArgumentException if n < 0
533 */
534 public static double factorialDouble(final int n) {
535 if (n < 0) {
536 throw MathRuntimeException.createIllegalArgumentException(
537 "must have n >= 0 for n!, got n = {0}",
538 n);
539 }
540 if (n < 21) {
541 return factorial(n);
542 }
543 return Math.floor(Math.exp(factorialLog(n)) + 0.5);
544 }
545
546 /**
547 * Returns the natural logarithm of n!.
548 * <p>
549 * <Strong>Preconditions</strong>:
550 * <ul>
551 * <li> <code>n >= 0</code> (otherwise
552 * <code>IllegalArgumentException</code> is thrown)</li>
553 * </ul></p>
554 *
555 * @param n argument
556 * @return <code>n!</code>
557 * @throws IllegalArgumentException if preconditions are not met.
558 */
559 public static double factorialLog(final int n) {
560 if (n < 0) {
561 throw MathRuntimeException.createIllegalArgumentException(
562 "must have n >= 0 for n!, got n = {0}",
563 n);
564 }
565 if (n < 21) {
566 return Math.log(factorial(n));
567 }
568 double logSum = 0;
569 for (int i = 2; i <= n; i++) {
570 logSum += Math.log(i);
571 }
572 return logSum;
573 }
574
575 /**
576 * <p>
577 * Gets the greatest common divisor of the absolute value of two numbers,
578 * using the "binary gcd" method which avoids division and modulo
579 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
580 * Stein (1961).
581 * </p>
582 * Special cases:
583 * <ul>
584 * <li>The invocations
585 * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>,
586 * <code>gcd(Integer.MIN_VALUE, 0)</code> and
587 * <code>gcd(0, Integer.MIN_VALUE)</code> throw an
588 * <code>ArithmeticException</code>, because the result would be 2^31, which
589 * is too large for an int value.</li>
590 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and
591 * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except
592 * for the special cases above.
593 * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
594 * <code>0</code>.</li>
595 * </ul>
596 *
597 * @param p any number
598 * @param q any number
599 * @return the greatest common divisor, never negative
600 * @throws ArithmeticException
601 * if the result cannot be represented as a nonnegative int
602 * value
603 * @since 1.1
604 */
605 public static int gcd(final int p, final int q) {
606 int u = p;
607 int v = q;
608 if ((u == 0) || (v == 0)) {
609 if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
610 throw MathRuntimeException.createArithmeticException(
611 "overflow: gcd({0}, {1}) is 2^31",
612 p, q);
613 }
614 return (Math.abs(u) + Math.abs(v));
615 }
616 // keep u and v negative, as negative integers range down to
617 // -2^31, while positive numbers can only be as large as 2^31-1
618 // (i.e. we can't necessarily negate a negative number without
619 // overflow)
620 /* assert u!=0 && v!=0; */
621 if (u > 0) {
622 u = -u;
623 } // make u negative
624 if (v > 0) {
625 v = -v;
626 } // make v negative
627 // B1. [Find power of 2]
628 int k = 0;
629 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are
630 // both even...
631 u /= 2;
632 v /= 2;
633 k++; // cast out twos.
634 }
635 if (k == 31) {
636 throw MathRuntimeException.createArithmeticException(
637 "overflow: gcd({0}, {1}) is 2^31",
638 p, q);
639 }
640 // B2. Initialize: u and v have been divided by 2^k and at least
641 // one is odd.
642 int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
643 // t negative: u was odd, v may be even (t replaces v)
644 // t positive: u was even, v is odd (t replaces u)
645 do {
646 /* assert u<0 && v<0; */
647 // B4/B3: cast out twos from t.
648 while ((t & 1) == 0) { // while t is even..
649 t /= 2; // cast out twos
650 }
651 // B5 [reset max(u,v)]
652 if (t > 0) {
653 u = -t;
654 } else {
655 v = t;
656 }
657 // B6/B3. at this point both u and v should be odd.
658 t = (v - u) / 2;
659 // |u| larger: t positive (replace u)
660 // |v| larger: t negative (replace v)
661 } while (t != 0);
662 return -u * (1 << k); // gcd is u*2^k
663 }
664
665 /**
666 * Returns an integer hash code representing the given double value.
667 *
668 * @param value the value to be hashed
669 * @return the hash code
670 */
671 public static int hash(double value) {
672 return new Double(value).hashCode();
673 }
674
675 /**
676 * Returns an integer hash code representing the given double array.
677 *
678 * @param value the value to be hashed (may be null)
679 * @return the hash code
680 * @since 1.2
681 */
682 public static int hash(double[] value) {
683 return Arrays.hashCode(value);
684 }
685
686 /**
687 * For a byte value x, this method returns (byte)(+1) if x >= 0 and
688 * (byte)(-1) if x < 0.
689 *
690 * @param x the value, a byte
691 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
692 */
693 public static byte indicator(final byte x) {
694 return (x >= ZB) ? PB : NB;
695 }
696
697 /**
698 * For a double precision value x, this method returns +1.0 if x >= 0 and
699 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
700 * <code>NaN</code>.
701 *
702 * @param x the value, a double
703 * @return +1.0 or -1.0, depending on the sign of x
704 */
705 public static double indicator(final double x) {
706 if (Double.isNaN(x)) {
707 return Double.NaN;
708 }
709 return (x >= 0.0) ? 1.0 : -1.0;
710 }
711
712 /**
713 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
714 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
715 *
716 * @param x the value, a float
717 * @return +1.0F or -1.0F, depending on the sign of x
718 */
719 public static float indicator(final float x) {
720 if (Float.isNaN(x)) {
721 return Float.NaN;
722 }
723 return (x >= 0.0F) ? 1.0F : -1.0F;
724 }
725
726 /**
727 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
728 *
729 * @param x the value, an int
730 * @return +1 or -1, depending on the sign of x
731 */
732 public static int indicator(final int x) {
733 return (x >= 0) ? 1 : -1;
734 }
735
736 /**
737 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
738 *
739 * @param x the value, a long
740 * @return +1L or -1L, depending on the sign of x
741 */
742 public static long indicator(final long x) {
743 return (x >= 0L) ? 1L : -1L;
744 }
745
746 /**
747 * For a short value x, this method returns (short)(+1) if x >= 0 and
748 * (short)(-1) if x < 0.
749 *
750 * @param x the value, a short
751 * @return (short)(+1) or (short)(-1), depending on the sign of x
752 */
753 public static short indicator(final short x) {
754 return (x >= ZS) ? PS : NS;
755 }
756
757 /**
758 * <p>
759 * Returns the least common multiple of the absolute value of two numbers,
760 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
761 * </p>
762 * Special cases:
763 * <ul>
764 * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and
765 * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a
766 * power of 2, throw an <code>ArithmeticException</code>, because the result
767 * would be 2^31, which is too large for an int value.</li>
768 * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is
769 * <code>0</code> for any <code>x</code>.
770 * </ul>
771 *
772 * @param a any number
773 * @param b any number
774 * @return the least common multiple, never negative
775 * @throws ArithmeticException
776 * if the result cannot be represented as a nonnegative int
777 * value
778 * @since 1.1
779 */
780 public static int lcm(int a, int b) {
781 if (a==0 || b==0){
782 return 0;
783 }
784 int lcm = Math.abs(mulAndCheck(a / gcd(a, b), b));
785 if (lcm == Integer.MIN_VALUE){
786 throw new ArithmeticException("overflow: lcm is 2^31");
787 }
788 return lcm;
789 }
790
791 /**
792 * <p>Returns the
793 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
794 * for base <code>b</code> of <code>x</code>.
795 * </p>
796 * <p>Returns <code>NaN<code> if either argument is negative. If
797 * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
798 * If <code>base</code> is positive and <code>x</code> is 0,
799 * <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments
800 * are 0, the result is <code>NaN</code>.</p>
801 *
802 * @param base the base of the logarithm, must be greater than 0
803 * @param x argument, must be greater than 0
804 * @return the value of the logarithm - the number y such that base^y = x.
805 * @since 1.2
806 */
807 public static double log(double base, double x) {
808 return Math.log(x)/Math.log(base);
809 }
810
811 /**
812 * Multiply two integers, checking for overflow.
813 *
814 * @param x a factor
815 * @param y a factor
816 * @return the product <code>x*y</code>
817 * @throws ArithmeticException if the result can not be represented as an
818 * int
819 * @since 1.1
820 */
821 public static int mulAndCheck(int x, int y) {
822 long m = ((long)x) * ((long)y);
823 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
824 throw new ArithmeticException("overflow: mul");
825 }
826 return (int)m;
827 }
828
829 /**
830 * Multiply two long integers, checking for overflow.
831 *
832 * @param a first value
833 * @param b second value
834 * @return the product <code>a * b</code>
835 * @throws ArithmeticException if the result can not be represented as an
836 * long
837 * @since 1.2
838 */
839 public static long mulAndCheck(long a, long b) {
840 long ret;
841 String msg = "overflow: multiply";
842 if (a > b) {
843 // use symmetry to reduce boundary cases
844 ret = mulAndCheck(b, a);
845 } else {
846 if (a < 0) {
847 if (b < 0) {
848 // check for positive overflow with negative a, negative b
849 if (a >= Long.MAX_VALUE / b) {
850 ret = a * b;
851 } else {
852 throw new ArithmeticException(msg);
853 }
854 } else if (b > 0) {
855 // check for negative overflow with negative a, positive b
856 if (Long.MIN_VALUE / b <= a) {
857 ret = a * b;
858 } else {
859 throw new ArithmeticException(msg);
860
861 }
862 } else {
863 // assert b == 0
864 ret = 0;
865 }
866 } else if (a > 0) {
867 // assert a > 0
868 // assert b > 0
869
870 // check for positive overflow with positive a, positive b
871 if (a <= Long.MAX_VALUE / b) {
872 ret = a * b;
873 } else {
874 throw new ArithmeticException(msg);
875 }
876 } else {
877 // assert a == 0
878 ret = 0;
879 }
880 }
881 return ret;
882 }
883
884 /**
885 * Get the next machine representable number after a number, moving
886 * in the direction of another number.
887 * <p>
888 * If <code>direction</code> is greater than or equal to<code>d</code>,
889 * the smallest machine representable number strictly greater than
890 * <code>d</code> is returned; otherwise the largest representable number
891 * strictly less than <code>d</code> is returned.</p>
892 * <p>
893 * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
894 *
895 * @param d base number
896 * @param direction (the only important thing is whether
897 * direction is greater or smaller than d)
898 * @return the next machine representable number in the specified direction
899 * @since 1.2
900 */
901 public static double nextAfter(double d, double direction) {
902
903 // handling of some important special cases
904 if (Double.isNaN(d) || Double.isInfinite(d)) {
905 return d;
906 } else if (d == 0) {
907 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
908 }
909 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
910 // are handled just as normal numbers
911
912 // split the double in raw components
913 long bits = Double.doubleToLongBits(d);
914 long sign = bits & 0x8000000000000000L;
915 long exponent = bits & 0x7ff0000000000000L;
916 long mantissa = bits & 0x000fffffffffffffL;
917
918 if (d * (direction - d) >= 0) {
919 // we should increase the mantissa
920 if (mantissa == 0x000fffffffffffffL) {
921 return Double.longBitsToDouble(sign |
922 (exponent + 0x0010000000000000L));
923 } else {
924 return Double.longBitsToDouble(sign |
925 exponent | (mantissa + 1));
926 }
927 } else {
928 // we should decrease the mantissa
929 if (mantissa == 0L) {
930 return Double.longBitsToDouble(sign |
931 (exponent - 0x0010000000000000L) |
932 0x000fffffffffffffL);
933 } else {
934 return Double.longBitsToDouble(sign |
935 exponent | (mantissa - 1));
936 }
937 }
938
939 }
940
941 /**
942 * Scale a number by 2<sup>scaleFactor</sup>.
943 * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p>
944 *
945 * @param d base number
946 * @param scaleFactor power of two by which d sould be multiplied
947 * @return d × 2<sup>scaleFactor</sup>
948 * @since 2.0
949 */
950 public static double scalb(final double d, final int scaleFactor) {
951
952 // handling of some important special cases
953 if ((d == 0) || Double.isNaN(d) || Double.isInfinite(d)) {
954 return d;
955 }
956
957 // split the double in raw components
958 final long bits = Double.doubleToLongBits(d);
959 final long exponent = bits & 0x7ff0000000000000L;
960 final long rest = bits & 0x800fffffffffffffL;
961
962 // shift the exponent
963 final long newBits = rest | (exponent + (((long) scaleFactor) << 52));
964 return Double.longBitsToDouble(newBits);
965
966 }
967
968 /**
969 * Normalize an angle in a 2&pi wide interval around a center value.
970 * <p>This method has three main uses:</p>
971 * <ul>
972 * <li>normalize an angle between 0 and 2π:<br/>
973 * <code>a = MathUtils.normalizeAngle(a, Math.PI);</code></li>
974 * <li>normalize an angle between -π and +π<br/>
975 * <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li>
976 * <li>compute the angle between two defining angular positions:<br>
977 * <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li>
978 * </ul>
979 * <p>Note that due to numerical accuracy and since π cannot be represented
980 * exactly, the result interval is <em>closed</em>, it cannot be half-closed
981 * as would be more satisfactory in a purely mathematical view.</p>
982 * @param a angle to normalize
983 * @param center center of the desired 2π interval for the result
984 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π
985 * @since 1.2
986 */
987 public static double normalizeAngle(double a, double center) {
988 return a - TWO_PI * Math.floor((a + Math.PI - center) / TWO_PI);
989 }
990
991 /**
992 * Round the given value to the specified number of decimal places. The
993 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
994 *
995 * @param x the value to round.
996 * @param scale the number of digits to the right of the decimal point.
997 * @return the rounded value.
998 * @since 1.1
999 */
1000 public static double round(double x, int scale) {
1001 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1002 }
1003
1004 /**
1005 * Round the given value to the specified number of decimal places. The
1006 * value is rounded using the given method which is any method defined in
1007 * {@link BigDecimal}.
1008 *
1009 * @param x the value to round.
1010 * @param scale the number of digits to the right of the decimal point.
1011 * @param roundingMethod the rounding method as defined in
1012 * {@link BigDecimal}.
1013 * @return the rounded value.
1014 * @since 1.1
1015 */
1016 public static double round(double x, int scale, int roundingMethod) {
1017 try {
1018 return (new BigDecimal
1019 (Double.toString(x))
1020 .setScale(scale, roundingMethod))
1021 .doubleValue();
1022 } catch (NumberFormatException ex) {
1023 if (Double.isInfinite(x)) {
1024 return x;
1025 } else {
1026 return Double.NaN;
1027 }
1028 }
1029 }
1030
1031 /**
1032 * Round the given value to the specified number of decimal places. The
1033 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
1034 *
1035 * @param x the value to round.
1036 * @param scale the number of digits to the right of the decimal point.
1037 * @return the rounded value.
1038 * @since 1.1
1039 */
1040 public static float round(float x, int scale) {
1041 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1042 }
1043
1044 /**
1045 * Round the given value to the specified number of decimal places. The
1046 * value is rounded using the given method which is any method defined in
1047 * {@link BigDecimal}.
1048 *
1049 * @param x the value to round.
1050 * @param scale the number of digits to the right of the decimal point.
1051 * @param roundingMethod the rounding method as defined in
1052 * {@link BigDecimal}.
1053 * @return the rounded value.
1054 * @since 1.1
1055 */
1056 public static float round(float x, int scale, int roundingMethod) {
1057 float sign = indicator(x);
1058 float factor = (float)Math.pow(10.0f, scale) * sign;
1059 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
1060 }
1061
1062 /**
1063 * Round the given non-negative, value to the "nearest" integer. Nearest is
1064 * determined by the rounding method specified. Rounding methods are defined
1065 * in {@link BigDecimal}.
1066 *
1067 * @param unscaled the value to round.
1068 * @param sign the sign of the original, scaled value.
1069 * @param roundingMethod the rounding method as defined in
1070 * {@link BigDecimal}.
1071 * @return the rounded value.
1072 * @since 1.1
1073 */
1074 private static double roundUnscaled(double unscaled, double sign,
1075 int roundingMethod) {
1076 switch (roundingMethod) {
1077 case BigDecimal.ROUND_CEILING :
1078 if (sign == -1) {
1079 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1080 } else {
1081 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1082 }
1083 break;
1084 case BigDecimal.ROUND_DOWN :
1085 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1086 break;
1087 case BigDecimal.ROUND_FLOOR :
1088 if (sign == -1) {
1089 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1090 } else {
1091 unscaled = Math.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1092 }
1093 break;
1094 case BigDecimal.ROUND_HALF_DOWN : {
1095 unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY);
1096 double fraction = unscaled - Math.floor(unscaled);
1097 if (fraction > 0.5) {
1098 unscaled = Math.ceil(unscaled);
1099 } else {
1100 unscaled = Math.floor(unscaled);
1101 }
1102 break;
1103 }
1104 case BigDecimal.ROUND_HALF_EVEN : {
1105 double fraction = unscaled - Math.floor(unscaled);
1106 if (fraction > 0.5) {
1107 unscaled = Math.ceil(unscaled);
1108 } else if (fraction < 0.5) {
1109 unscaled = Math.floor(unscaled);
1110 } else {
1111 // The following equality test is intentional and needed for rounding purposes
1112 if (Math.floor(unscaled) / 2.0 == Math.floor(Math
1113 .floor(unscaled) / 2.0)) { // even
1114 unscaled = Math.floor(unscaled);
1115 } else { // odd
1116 unscaled = Math.ceil(unscaled);
1117 }
1118 }
1119 break;
1120 }
1121 case BigDecimal.ROUND_HALF_UP : {
1122 unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY);
1123 double fraction = unscaled - Math.floor(unscaled);
1124 if (fraction >= 0.5) {
1125 unscaled = Math.ceil(unscaled);
1126 } else {
1127 unscaled = Math.floor(unscaled);
1128 }
1129 break;
1130 }
1131 case BigDecimal.ROUND_UNNECESSARY :
1132 if (unscaled != Math.floor(unscaled)) {
1133 throw new ArithmeticException("Inexact result from rounding");
1134 }
1135 break;
1136 case BigDecimal.ROUND_UP :
1137 unscaled = Math.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1138 break;
1139 default :
1140 throw MathRuntimeException.createIllegalArgumentException(
1141 "invalid rounding method {0}, valid methods: {1} ({2}), {3} ({4})," +
1142 " {5} ({6}), {7} ({8}), {9} ({10}), {11} ({12}), {13} ({14}), {15} ({16})",
1143 roundingMethod,
1144 "ROUND_CEILING", BigDecimal.ROUND_CEILING,
1145 "ROUND_DOWN", BigDecimal.ROUND_DOWN,
1146 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
1147 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
1148 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
1149 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
1150 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
1151 "ROUND_UP", BigDecimal.ROUND_UP);
1152 }
1153 return unscaled;
1154 }
1155
1156 /**
1157 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1158 * for byte value <code>x</code>.
1159 * <p>
1160 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
1161 * x = 0, and (byte)(-1) if x < 0.</p>
1162 *
1163 * @param x the value, a byte
1164 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
1165 */
1166 public static byte sign(final byte x) {
1167 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
1168 }
1169
1170 /**
1171 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1172 * for double precision <code>x</code>.
1173 * <p>
1174 * For a double value <code>x</code>, this method returns
1175 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
1176 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
1177 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
1178 *
1179 * @param x the value, a double
1180 * @return +1.0, 0.0, or -1.0, depending on the sign of x
1181 */
1182 public static double sign(final double x) {
1183 if (Double.isNaN(x)) {
1184 return Double.NaN;
1185 }
1186 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
1187 }
1188
1189 /**
1190 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1191 * for float value <code>x</code>.
1192 * <p>
1193 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
1194 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
1195 * is <code>NaN</code>.</p>
1196 *
1197 * @param x the value, a float
1198 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
1199 */
1200 public static float sign(final float x) {
1201 if (Float.isNaN(x)) {
1202 return Float.NaN;
1203 }
1204 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
1205 }
1206
1207 /**
1208 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1209 * for int value <code>x</code>.
1210 * <p>
1211 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
1212 * if x < 0.</p>
1213 *
1214 * @param x the value, an int
1215 * @return +1, 0, or -1, depending on the sign of x
1216 */
1217 public static int sign(final int x) {
1218 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
1219 }
1220
1221 /**
1222 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1223 * for long value <code>x</code>.
1224 * <p>
1225 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
1226 * -1L if x < 0.</p>
1227 *
1228 * @param x the value, a long
1229 * @return +1L, 0L, or -1L, depending on the sign of x
1230 */
1231 public static long sign(final long x) {
1232 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
1233 }
1234
1235 /**
1236 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1237 * for short value <code>x</code>.
1238 * <p>
1239 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
1240 * if x = 0, and (short)(-1) if x < 0.</p>
1241 *
1242 * @param x the value, a short
1243 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
1244 * x
1245 */
1246 public static short sign(final short x) {
1247 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
1248 }
1249
1250 /**
1251 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
1252 * hyperbolic sine</a> of x.
1253 *
1254 * @param x double value for which to find the hyperbolic sine
1255 * @return hyperbolic sine of x
1256 */
1257 public static double sinh(double x) {
1258 return (Math.exp(x) - Math.exp(-x)) / 2.0;
1259 }
1260
1261 /**
1262 * Subtract two integers, checking for overflow.
1263 *
1264 * @param x the minuend
1265 * @param y the subtrahend
1266 * @return the difference <code>x-y</code>
1267 * @throws ArithmeticException if the result can not be represented as an
1268 * int
1269 * @since 1.1
1270 */
1271 public static int subAndCheck(int x, int y) {
1272 long s = (long)x - (long)y;
1273 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
1274 throw new ArithmeticException("overflow: subtract");
1275 }
1276 return (int)s;
1277 }
1278
1279 /**
1280 * Subtract two long integers, checking for overflow.
1281 *
1282 * @param a first value
1283 * @param b second value
1284 * @return the difference <code>a-b</code>
1285 * @throws ArithmeticException if the result can not be represented as an
1286 * long
1287 * @since 1.2
1288 */
1289 public static long subAndCheck(long a, long b) {
1290 long ret;
1291 String msg = "overflow: subtract";
1292 if (b == Long.MIN_VALUE) {
1293 if (a < 0) {
1294 ret = a - b;
1295 } else {
1296 throw new ArithmeticException(msg);
1297 }
1298 } else {
1299 // use additive inverse
1300 ret = addAndCheck(a, -b, msg);
1301 }
1302 return ret;
1303 }
1304
1305 /**
1306 * Raise an int to an int power.
1307 * @param k number to raise
1308 * @param e exponent (must be positive or null)
1309 * @return k<sup>e</sup>
1310 * @exception IllegalArgumentException if e is negative
1311 */
1312 public static int pow(final int k, int e)
1313 throws IllegalArgumentException {
1314
1315 if (e < 0) {
1316 throw MathRuntimeException.createIllegalArgumentException(
1317 "cannot raise an integral value to a negative power ({0}^{1})",
1318 k, e);
1319 }
1320
1321 int result = 1;
1322 int k2p = k;
1323 while (e != 0) {
1324 if ((e & 0x1) != 0) {
1325 result *= k2p;
1326 }
1327 k2p *= k2p;
1328 e = e >> 1;
1329 }
1330
1331 return result;
1332
1333 }
1334
1335 /**
1336 * Raise an int to a long power.
1337 * @param k number to raise
1338 * @param e exponent (must be positive or null)
1339 * @return k<sup>e</sup>
1340 * @exception IllegalArgumentException if e is negative
1341 */
1342 public static int pow(final int k, long e)
1343 throws IllegalArgumentException {
1344
1345 if (e < 0) {
1346 throw MathRuntimeException.createIllegalArgumentException(
1347 "cannot raise an integral value to a negative power ({0}^{1})",
1348 k, e);
1349 }
1350
1351 int result = 1;
1352 int k2p = k;
1353 while (e != 0) {
1354 if ((e & 0x1) != 0) {
1355 result *= k2p;
1356 }
1357 k2p *= k2p;
1358 e = e >> 1;
1359 }
1360
1361 return result;
1362
1363 }
1364
1365 /**
1366 * Raise a long to an int power.
1367 * @param k number to raise
1368 * @param e exponent (must be positive or null)
1369 * @return k<sup>e</sup>
1370 * @exception IllegalArgumentException if e is negative
1371 */
1372 public static long pow(final long k, int e)
1373 throws IllegalArgumentException {
1374
1375 if (e < 0) {
1376 throw MathRuntimeException.createIllegalArgumentException(
1377 "cannot raise an integral value to a negative power ({0}^{1})",
1378 k, e);
1379 }
1380
1381 long result = 1l;
1382 long k2p = k;
1383 while (e != 0) {
1384 if ((e & 0x1) != 0) {
1385 result *= k2p;
1386 }
1387 k2p *= k2p;
1388 e = e >> 1;
1389 }
1390
1391 return result;
1392
1393 }
1394
1395 /**
1396 * Raise a long to a long power.
1397 * @param k number to raise
1398 * @param e exponent (must be positive or null)
1399 * @return k<sup>e</sup>
1400 * @exception IllegalArgumentException if e is negative
1401 */
1402 public static long pow(final long k, long e)
1403 throws IllegalArgumentException {
1404
1405 if (e < 0) {
1406 throw MathRuntimeException.createIllegalArgumentException(
1407 "cannot raise an integral value to a negative power ({0}^{1})",
1408 k, e);
1409 }
1410
1411 long result = 1l;
1412 long k2p = k;
1413 while (e != 0) {
1414 if ((e & 0x1) != 0) {
1415 result *= k2p;
1416 }
1417 k2p *= k2p;
1418 e = e >> 1;
1419 }
1420
1421 return result;
1422
1423 }
1424
1425 /**
1426 * Raise a BigInteger to an int power.
1427 * @param k number to raise
1428 * @param e exponent (must be positive or null)
1429 * @return k<sup>e</sup>
1430 * @exception IllegalArgumentException if e is negative
1431 */
1432 public static BigInteger pow(final BigInteger k, int e)
1433 throws IllegalArgumentException {
1434
1435 if (e < 0) {
1436 throw MathRuntimeException.createIllegalArgumentException(
1437 "cannot raise an integral value to a negative power ({0}^{1})",
1438 k, e);
1439 }
1440
1441 return k.pow(e);
1442
1443 }
1444
1445 /**
1446 * Raise a BigInteger to a long power.
1447 * @param k number to raise
1448 * @param e exponent (must be positive or null)
1449 * @return k<sup>e</sup>
1450 * @exception IllegalArgumentException if e is negative
1451 */
1452 public static BigInteger pow(final BigInteger k, long e)
1453 throws IllegalArgumentException {
1454
1455 if (e < 0) {
1456 throw MathRuntimeException.createIllegalArgumentException(
1457 "cannot raise an integral value to a negative power ({0}^{1})",
1458 k, e);
1459 }
1460
1461 BigInteger result = BigInteger.ONE;
1462 BigInteger k2p = k;
1463 while (e != 0) {
1464 if ((e & 0x1) != 0) {
1465 result = result.multiply(k2p);
1466 }
1467 k2p = k2p.multiply(k2p);
1468 e = e >> 1;
1469 }
1470
1471 return result;
1472
1473 }
1474
1475 /**
1476 * Raise a BigInteger to a BigInteger power.
1477 * @param k number to raise
1478 * @param e exponent (must be positive or null)
1479 * @return k<sup>e</sup>
1480 * @exception IllegalArgumentException if e is negative
1481 */
1482 public static BigInteger pow(final BigInteger k, BigInteger e)
1483 throws IllegalArgumentException {
1484
1485 if (e.compareTo(BigInteger.ZERO) < 0) {
1486 throw MathRuntimeException.createIllegalArgumentException(
1487 "cannot raise an integral value to a negative power ({0}^{1})",
1488 k, e);
1489 }
1490
1491 BigInteger result = BigInteger.ONE;
1492 BigInteger k2p = k;
1493 while (!BigInteger.ZERO.equals(e)) {
1494 if (e.testBit(0)) {
1495 result = result.multiply(k2p);
1496 }
1497 k2p = k2p.multiply(k2p);
1498 e = e.shiftRight(1);
1499 }
1500
1501 return result;
1502
1503 }
1504
1505 /**
1506 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1507 *
1508 * @param p1 the first point
1509 * @param p2 the second point
1510 * @return the L<sub>1</sub> distance between the two points
1511 */
1512 public static final double distance1(double[] p1, double[] p2) {
1513 double sum = 0;
1514 for (int i = 0; i < p1.length; i++) {
1515 sum += Math.abs(p1[i] - p2[i]);
1516 }
1517 return sum;
1518 }
1519
1520 /**
1521 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1522 *
1523 * @param p1 the first point
1524 * @param p2 the second point
1525 * @return the L<sub>1</sub> distance between the two points
1526 */
1527 public static final int distance1(int[] p1, int[] p2) {
1528 int sum = 0;
1529 for (int i = 0; i < p1.length; i++) {
1530 sum += Math.abs(p1[i] - p2[i]);
1531 }
1532 return sum;
1533 }
1534
1535 /**
1536 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1537 *
1538 * @param p1 the first point
1539 * @param p2 the second point
1540 * @return the L<sub>2</sub> distance between the two points
1541 */
1542 public static final double distance(double[] p1, double[] p2) {
1543 double sum = 0;
1544 for (int i = 0; i < p1.length; i++) {
1545 final double dp = p1[i] - p2[i];
1546 sum += dp * dp;
1547 }
1548 return Math.sqrt(sum);
1549 }
1550
1551 /**
1552 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
1553 *
1554 * @param p1 the first point
1555 * @param p2 the second point
1556 * @return the L<sub>2</sub> distance between the two points
1557 */
1558 public static final double distance(int[] p1, int[] p2) {
1559 int sum = 0;
1560 for (int i = 0; i < p1.length; i++) {
1561 final int dp = p1[i] - p2[i];
1562 sum += dp * dp;
1563 }
1564 return Math.sqrt(sum);
1565 }
1566
1567 /**
1568 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1569 *
1570 * @param p1 the first point
1571 * @param p2 the second point
1572 * @return the L<sub>∞</sub> distance between the two points
1573 */
1574 public static final double distanceInf(double[] p1, double[] p2) {
1575 double max = 0;
1576 for (int i = 0; i < p1.length; i++) {
1577 max = Math.max(max, Math.abs(p1[i] - p2[i]));
1578 }
1579 return max;
1580 }
1581
1582 /**
1583 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
1584 *
1585 * @param p1 the first point
1586 * @param p2 the second point
1587 * @return the L<sub>∞</sub> distance between the two points
1588 */
1589 public static final int distanceInf(int[] p1, int[] p2) {
1590 int max = 0;
1591 for (int i = 0; i < p1.length; i++) {
1592 max = Math.max(max, Math.abs(p1[i] - p2[i]));
1593 }
1594 return max;
1595 }
1596
1597
1598 }