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 import org.apache.commons.math.exception.util.Localizable;
026 import org.apache.commons.math.exception.util.LocalizedFormats;
027 import org.apache.commons.math.exception.NonMonotonousSequenceException;
028
029 /**
030 * Some useful additions to the built-in functions in {@link Math}.
031 * @version $Revision: 1073472 $ $Date: 2011-02-22 20:49:07 +0100 (mar. 22 f??vr. 2011) $
032 */
033 public final class MathUtils {
034
035 /** Smallest positive number such that 1 - EPSILON is not numerically equal to 1. */
036 public static final double EPSILON = 0x1.0p-53;
037
038 /** Safe minimum, such that 1 / SAFE_MIN does not overflow.
039 * <p>In IEEE 754 arithmetic, this is also the smallest normalized
040 * number 2<sup>-1022</sup>.</p>
041 */
042 public static final double SAFE_MIN = 0x1.0p-1022;
043
044 /**
045 * 2 π.
046 * @since 2.1
047 */
048 public static final double TWO_PI = 2 * FastMath.PI;
049
050 /** -1.0 cast as a byte. */
051 private static final byte NB = (byte)-1;
052
053 /** -1.0 cast as a short. */
054 private static final short NS = (short)-1;
055
056 /** 1.0 cast as a byte. */
057 private static final byte PB = (byte)1;
058
059 /** 1.0 cast as a short. */
060 private static final short PS = (short)1;
061
062 /** 0.0 cast as a byte. */
063 private static final byte ZB = (byte)0;
064
065 /** 0.0 cast as a short. */
066 private static final short ZS = (short)0;
067
068 /** Gap between NaN and regular numbers. */
069 private static final int NAN_GAP = 4 * 1024 * 1024;
070
071 /** Offset to order signed double numbers lexicographically. */
072 private static final long SGN_MASK = 0x8000000000000000L;
073
074 /** Offset to order signed double numbers lexicographically. */
075 private static final int SGN_MASK_FLOAT = 0x80000000;
076
077 /** All long-representable factorials */
078 private static final long[] FACTORIALS = new long[] {
079 1l, 1l, 2l,
080 6l, 24l, 120l,
081 720l, 5040l, 40320l,
082 362880l, 3628800l, 39916800l,
083 479001600l, 6227020800l, 87178291200l,
084 1307674368000l, 20922789888000l, 355687428096000l,
085 6402373705728000l, 121645100408832000l, 2432902008176640000l };
086
087 /**
088 * Private Constructor
089 */
090 private MathUtils() {
091 super();
092 }
093
094 /**
095 * Add two integers, checking for overflow.
096 *
097 * @param x an addend
098 * @param y an addend
099 * @return the sum <code>x+y</code>
100 * @throws ArithmeticException if the result can not be represented as an
101 * int
102 * @since 1.1
103 */
104 public static int addAndCheck(int x, int y) {
105 long s = (long)x + (long)y;
106 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
107 throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, x, y);
108 }
109 return (int)s;
110 }
111
112 /**
113 * Add two long integers, checking for overflow.
114 *
115 * @param a an addend
116 * @param b an addend
117 * @return the sum <code>a+b</code>
118 * @throws ArithmeticException if the result can not be represented as an
119 * long
120 * @since 1.2
121 */
122 public static long addAndCheck(long a, long b) {
123 return addAndCheck(a, b, LocalizedFormats.OVERFLOW_IN_ADDITION);
124 }
125
126 /**
127 * Add two long integers, checking for overflow.
128 *
129 * @param a an addend
130 * @param b an addend
131 * @param pattern the pattern to use for any thrown exception.
132 * @return the sum <code>a+b</code>
133 * @throws ArithmeticException if the result can not be represented as an
134 * long
135 * @since 1.2
136 */
137 private static long addAndCheck(long a, long b, Localizable pattern) {
138 long ret;
139 if (a > b) {
140 // use symmetry to reduce boundary cases
141 ret = addAndCheck(b, a, pattern);
142 } else {
143 // assert a <= b
144
145 if (a < 0) {
146 if (b < 0) {
147 // check for negative overflow
148 if (Long.MIN_VALUE - b <= a) {
149 ret = a + b;
150 } else {
151 throw MathRuntimeException.createArithmeticException(pattern, a, b);
152 }
153 } else {
154 // opposite sign addition is always safe
155 ret = a + b;
156 }
157 } else {
158 // assert a >= 0
159 // assert b >= 0
160
161 // check for positive overflow
162 if (a <= Long.MAX_VALUE - b) {
163 ret = a + b;
164 } else {
165 throw MathRuntimeException.createArithmeticException(pattern, a, b);
166 }
167 }
168 }
169 return ret;
170 }
171
172 /**
173 * Returns an exact representation of the <a
174 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
175 * Coefficient</a>, "<code>n choose k</code>", the number of
176 * <code>k</code>-element subsets that can be selected from an
177 * <code>n</code>-element set.
178 * <p>
179 * <Strong>Preconditions</strong>:
180 * <ul>
181 * <li> <code>0 <= k <= n </code> (otherwise
182 * <code>IllegalArgumentException</code> is thrown)</li>
183 * <li> The result is small enough to fit into a <code>long</code>. The
184 * largest value of <code>n</code> for which all coefficients are
185 * <code> < Long.MAX_VALUE</code> is 66. If the computed value exceeds
186 * <code>Long.MAX_VALUE</code> an <code>ArithMeticException</code> is
187 * thrown.</li>
188 * </ul></p>
189 *
190 * @param n the size of the set
191 * @param k the size of the subsets to be counted
192 * @return <code>n choose k</code>
193 * @throws IllegalArgumentException if preconditions are not met.
194 * @throws ArithmeticException if the result is too large to be represented
195 * by a long integer.
196 */
197 public static long binomialCoefficient(final int n, final int k) {
198 checkBinomial(n, k);
199 if ((n == k) || (k == 0)) {
200 return 1;
201 }
202 if ((k == 1) || (k == n - 1)) {
203 return n;
204 }
205 // Use symmetry for large k
206 if (k > n / 2)
207 return binomialCoefficient(n, n - k);
208
209 // We use the formula
210 // (n choose k) = n! / (n-k)! / k!
211 // (n choose k) == ((n-k+1)*...*n) / (1*...*k)
212 // which could be written
213 // (n choose k) == (n-1 choose k-1) * n / k
214 long result = 1;
215 if (n <= 61) {
216 // For n <= 61, the naive implementation cannot overflow.
217 int i = n - k + 1;
218 for (int j = 1; j <= k; j++) {
219 result = result * i / j;
220 i++;
221 }
222 } else if (n <= 66) {
223 // For n > 61 but n <= 66, the result cannot overflow,
224 // but we must take care not to overflow intermediate values.
225 int i = n - k + 1;
226 for (int j = 1; j <= k; j++) {
227 // We know that (result * i) is divisible by j,
228 // but (result * i) may overflow, so we split j:
229 // Filter out the gcd, d, so j/d and i/d are integer.
230 // result is divisible by (j/d) because (j/d)
231 // is relative prime to (i/d) and is a divisor of
232 // result * (i/d).
233 final long d = gcd(i, j);
234 result = (result / (j / d)) * (i / d);
235 i++;
236 }
237 } else {
238 // For n > 66, a result overflow might occur, so we check
239 // the multiplication, taking care to not overflow
240 // unnecessary.
241 int i = n - k + 1;
242 for (int j = 1; j <= k; j++) {
243 final long d = gcd(i, j);
244 result = mulAndCheck(result / (j / d), i / d);
245 i++;
246 }
247 }
248 return result;
249 }
250
251 /**
252 * Returns a <code>double</code> representation of the <a
253 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
254 * Coefficient</a>, "<code>n choose k</code>", the number of
255 * <code>k</code>-element subsets that can be selected from an
256 * <code>n</code>-element set.
257 * <p>
258 * <Strong>Preconditions</strong>:
259 * <ul>
260 * <li> <code>0 <= k <= n </code> (otherwise
261 * <code>IllegalArgumentException</code> is thrown)</li>
262 * <li> The result is small enough to fit into a <code>double</code>. The
263 * largest value of <code>n</code> for which all coefficients are <
264 * Double.MAX_VALUE is 1029. If the computed value exceeds Double.MAX_VALUE,
265 * Double.POSITIVE_INFINITY is returned</li>
266 * </ul></p>
267 *
268 * @param n the size of the set
269 * @param k the size of the subsets to be counted
270 * @return <code>n choose k</code>
271 * @throws IllegalArgumentException if preconditions are not met.
272 */
273 public static double binomialCoefficientDouble(final int n, final int k) {
274 checkBinomial(n, k);
275 if ((n == k) || (k == 0)) {
276 return 1d;
277 }
278 if ((k == 1) || (k == n - 1)) {
279 return n;
280 }
281 if (k > n/2) {
282 return binomialCoefficientDouble(n, n - k);
283 }
284 if (n < 67) {
285 return binomialCoefficient(n,k);
286 }
287
288 double result = 1d;
289 for (int i = 1; i <= k; i++) {
290 result *= (double)(n - k + i) / (double)i;
291 }
292
293 return FastMath.floor(result + 0.5);
294 }
295
296 /**
297 * Returns the natural <code>log</code> of the <a
298 * href="http://mathworld.wolfram.com/BinomialCoefficient.html"> Binomial
299 * Coefficient</a>, "<code>n choose k</code>", the number of
300 * <code>k</code>-element subsets that can be selected from an
301 * <code>n</code>-element set.
302 * <p>
303 * <Strong>Preconditions</strong>:
304 * <ul>
305 * <li> <code>0 <= k <= n </code> (otherwise
306 * <code>IllegalArgumentException</code> is thrown)</li>
307 * </ul></p>
308 *
309 * @param n the size of the set
310 * @param k the size of the subsets to be counted
311 * @return <code>n choose k</code>
312 * @throws IllegalArgumentException if preconditions are not met.
313 */
314 public static double binomialCoefficientLog(final int n, final int k) {
315 checkBinomial(n, k);
316 if ((n == k) || (k == 0)) {
317 return 0;
318 }
319 if ((k == 1) || (k == n - 1)) {
320 return FastMath.log(n);
321 }
322
323 /*
324 * For values small enough to do exact integer computation,
325 * return the log of the exact value
326 */
327 if (n < 67) {
328 return FastMath.log(binomialCoefficient(n,k));
329 }
330
331 /*
332 * Return the log of binomialCoefficientDouble for values that will not
333 * overflow binomialCoefficientDouble
334 */
335 if (n < 1030) {
336 return FastMath.log(binomialCoefficientDouble(n, k));
337 }
338
339 if (k > n / 2) {
340 return binomialCoefficientLog(n, n - k);
341 }
342
343 /*
344 * Sum logs for values that could overflow
345 */
346 double logSum = 0;
347
348 // n!/(n-k)!
349 for (int i = n - k + 1; i <= n; i++) {
350 logSum += FastMath.log(i);
351 }
352
353 // divide by k!
354 for (int i = 2; i <= k; i++) {
355 logSum -= FastMath.log(i);
356 }
357
358 return logSum;
359 }
360
361 /**
362 * Check binomial preconditions.
363 * @param n the size of the set
364 * @param k the size of the subsets to be counted
365 * @exception IllegalArgumentException if preconditions are not met.
366 */
367 private static void checkBinomial(final int n, final int k)
368 throws IllegalArgumentException {
369 if (n < k) {
370 throw MathRuntimeException.createIllegalArgumentException(
371 LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
372 n, k);
373 }
374 if (n < 0) {
375 throw MathRuntimeException.createIllegalArgumentException(
376 LocalizedFormats.BINOMIAL_NEGATIVE_PARAMETER,
377 n);
378 }
379 }
380
381 /**
382 * Compares two numbers given some amount of allowed error.
383 *
384 * @param x the first number
385 * @param y the second number
386 * @param eps the amount of error to allow when checking for equality
387 * @return <ul><li>0 if {@link #equals(double, double, double) equals(x, y, eps)}</li>
388 * <li>< 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x < y</li>
389 * <li>> 0 if !{@link #equals(double, double, double) equals(x, y, eps)} && x > y</li></ul>
390 */
391 public static int compareTo(double x, double y, double eps) {
392 if (equals(x, y, eps)) {
393 return 0;
394 } else if (x < y) {
395 return -1;
396 }
397 return 1;
398 }
399
400 /**
401 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicCosine.html">
402 * hyperbolic cosine</a> of x.
403 *
404 * @param x double value for which to find the hyperbolic cosine
405 * @return hyperbolic cosine of x
406 */
407 public static double cosh(double x) {
408 return (FastMath.exp(x) + FastMath.exp(-x)) / 2.0;
409 }
410
411 /**
412 * Returns true iff they are strictly equal.
413 *
414 * @param x first value
415 * @param y second value
416 * @return {@code true} if the values are equal.
417 * @deprecated as of 2.2 his method considers that {@code NaN == NaN}. In release
418 * 3.0, the semantics will change in order to comply with IEEE754 where it
419 * is specified that {@code NaN != NaN}.
420 * New methods have been added for those cases wher the old semantics is
421 * useful (see e.g. {@link #equalsIncludingNaN(float,float)
422 * equalsIncludingNaN}.
423 */
424 @Deprecated
425 public static boolean equals(float x, float y) {
426 return (Float.isNaN(x) && Float.isNaN(y)) || x == y;
427 }
428
429 /**
430 * Returns true if both arguments are NaN or neither is NaN and they are
431 * equal as defined by {@link #equals(float,float,int) equals(x, y, 1)}.
432 *
433 * @param x first value
434 * @param y second value
435 * @return {@code true} if the values are equal or both are NaN.
436 * @since 2.2
437 */
438 public static boolean equalsIncludingNaN(float x, float y) {
439 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, 1);
440 }
441
442 /**
443 * Returns true if both arguments are equal or within the range of allowed
444 * error (inclusive).
445 *
446 * @param x first value
447 * @param y second value
448 * @param eps the amount of absolute error to allow.
449 * @return {@code true} if the values are equal or within range of each other.
450 * @since 2.2
451 */
452 public static boolean equals(float x, float y, float eps) {
453 return equals(x, y, 1) || FastMath.abs(y - x) <= eps;
454 }
455
456 /**
457 * Returns true if both arguments are NaN or are equal or within the range
458 * of allowed error (inclusive).
459 *
460 * @param x first value
461 * @param y second value
462 * @param eps the amount of absolute error to allow.
463 * @return {@code true} if the values are equal or within range of each other,
464 * or both are NaN.
465 * @since 2.2
466 */
467 public static boolean equalsIncludingNaN(float x, float y, float eps) {
468 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
469 }
470
471 /**
472 * Returns true if both arguments are equal or within the range of allowed
473 * error (inclusive).
474 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
475 * (or fewer) floating point numbers between them, i.e. two adjacent floating
476 * point numbers are considered equal.
477 * Adapted from <a
478 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
479 * Bruce Dawson</a>
480 *
481 * @param x first value
482 * @param y second value
483 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
484 * values between {@code x} and {@code y}.
485 * @return {@code true} if there are fewer than {@code maxUlps} floating
486 * point values between {@code x} and {@code y}.
487 * @since 2.2
488 */
489 public static boolean equals(float x, float y, int maxUlps) {
490 // Check that "maxUlps" is non-negative and small enough so that
491 // NaN won't compare as equal to anything (except another NaN).
492 assert maxUlps > 0 && maxUlps < NAN_GAP;
493
494 int xInt = Float.floatToIntBits(x);
495 int yInt = Float.floatToIntBits(y);
496
497 // Make lexicographically ordered as a two's-complement integer.
498 if (xInt < 0) {
499 xInt = SGN_MASK_FLOAT - xInt;
500 }
501 if (yInt < 0) {
502 yInt = SGN_MASK_FLOAT - yInt;
503 }
504
505 final boolean isEqual = FastMath.abs(xInt - yInt) <= maxUlps;
506
507 return isEqual && !Float.isNaN(x) && !Float.isNaN(y);
508 }
509
510 /**
511 * Returns true if both arguments are NaN or if they are equal as defined
512 * by {@link #equals(float,float,int) equals(x, y, maxUlps)}.
513 *
514 * @param x first value
515 * @param y second value
516 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
517 * values between {@code x} and {@code y}.
518 * @return {@code true} if both arguments are NaN or if there are less than
519 * {@code maxUlps} floating point values between {@code x} and {@code y}.
520 * @since 2.2
521 */
522 public static boolean equalsIncludingNaN(float x, float y, int maxUlps) {
523 return (Float.isNaN(x) && Float.isNaN(y)) || equals(x, y, maxUlps);
524 }
525
526 /**
527 * Returns true iff both arguments are null or have same dimensions and all
528 * their elements are equal as defined by
529 * {@link #equals(float,float)}.
530 *
531 * @param x first array
532 * @param y second array
533 * @return true if the values are both null or have same dimension
534 * and equal elements.
535 * @deprecated as of 2.2 this method considers that {@code NaN == NaN}. In release
536 * 3.0, the semantics will change in order to comply with IEEE754 where it
537 * is specified that {@code NaN != NaN}.
538 * New methods have been added for those cases where the old semantics is
539 * useful (see e.g. {@link #equalsIncludingNaN(float[],float[])
540 * equalsIncludingNaN}.
541 */
542 @Deprecated
543 public static boolean equals(float[] x, float[] y) {
544 if ((x == null) || (y == null)) {
545 return !((x == null) ^ (y == null));
546 }
547 if (x.length != y.length) {
548 return false;
549 }
550 for (int i = 0; i < x.length; ++i) {
551 if (!equals(x[i], y[i])) {
552 return false;
553 }
554 }
555 return true;
556 }
557
558 /**
559 * Returns true iff both arguments are null or have same dimensions and all
560 * their elements are equal as defined by
561 * {@link #equalsIncludingNaN(float,float)}.
562 *
563 * @param x first array
564 * @param y second array
565 * @return true if the values are both null or have same dimension and
566 * equal elements
567 * @since 2.2
568 */
569 public static boolean equalsIncludingNaN(float[] x, float[] y) {
570 if ((x == null) || (y == null)) {
571 return !((x == null) ^ (y == null));
572 }
573 if (x.length != y.length) {
574 return false;
575 }
576 for (int i = 0; i < x.length; ++i) {
577 if (!equalsIncludingNaN(x[i], y[i])) {
578 return false;
579 }
580 }
581 return true;
582 }
583
584 /**
585 * Returns true iff both arguments are NaN or neither is NaN and they are
586 * equal
587 *
588 * <p>This method considers that {@code NaN == NaN}. In release
589 * 3.0, the semantics will change in order to comply with IEEE754 where it
590 * is specified that {@code NaN != NaN}.
591 * New methods have been added for those cases where the old semantics
592 * (w.r.t. NaN) is useful (see e.g.
593 * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}.
594 * </p>
595 *
596 * @param x first value
597 * @param y second value
598 * @return {@code true} if the values are equal.
599 */
600 public static boolean equals(double x, double y) {
601 return (Double.isNaN(x) && Double.isNaN(y)) || x == y;
602 }
603
604 /**
605 * Returns true if both arguments are NaN or neither is NaN and they are
606 * equal as defined by {@link #equals(double,double,int) equals(x, y, 1)}.
607 *
608 * @param x first value
609 * @param y second value
610 * @return {@code true} if the values are equal or both are NaN.
611 * @since 2.2
612 */
613 public static boolean equalsIncludingNaN(double x, double y) {
614 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, 1);
615 }
616
617 /**
618 * Returns true if both arguments are equal or within the range of allowed
619 * error (inclusive).
620 * <p>
621 * Two NaNs are considered equals, as are two infinities with same sign.
622 * </p>
623 * <p>This method considers that {@code NaN == NaN}. In release
624 * 3.0, the semantics will change in order to comply with IEEE754 where it
625 * is specified that {@code NaN != NaN}.
626 * New methods have been added for those cases where the old semantics
627 * (w.r.t. NaN) is useful (see e.g.
628 * {@link #equalsIncludingNaN(double,double, double) equalsIncludingNaN}.
629 * </p>
630 * @param x first value
631 * @param y second value
632 * @param eps the amount of absolute error to allow.
633 * @return {@code true} if the values are equal or within range of each other.
634 */
635 public static boolean equals(double x, double y, double eps) {
636 return equals(x, y) || FastMath.abs(y - x) <= eps;
637 }
638
639 /**
640 * Returns true if both arguments are NaN or are equal or within the range
641 * of allowed error (inclusive).
642 *
643 * @param x first value
644 * @param y second value
645 * @param eps the amount of absolute error to allow.
646 * @return {@code true} if the values are equal or within range of each other,
647 * or both are NaN.
648 * @since 2.2
649 */
650 public static boolean equalsIncludingNaN(double x, double y, double eps) {
651 return equalsIncludingNaN(x, y) || (FastMath.abs(y - x) <= eps);
652 }
653
654 /**
655 * Returns true if both arguments are equal or within the range of allowed
656 * error (inclusive).
657 * Two float numbers are considered equal if there are {@code (maxUlps - 1)}
658 * (or fewer) floating point numbers between them, i.e. two adjacent floating
659 * point numbers are considered equal.
660 * Adapted from <a
661 * href="http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm">
662 * Bruce Dawson</a>
663 *
664 * <p>This method considers that {@code NaN == NaN}. In release
665 * 3.0, the semantics will change in order to comply with IEEE754 where it
666 * is specified that {@code NaN != NaN}.
667 * New methods have been added for those cases where the old semantics
668 * (w.r.t. NaN) is useful (see e.g.
669 * {@link #equalsIncludingNaN(double,double, int) equalsIncludingNaN}.
670 * </p>
671 *
672 * @param x first value
673 * @param y second value
674 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
675 * values between {@code x} and {@code y}.
676 * @return {@code true} if there are fewer than {@code maxUlps} floating
677 * point values between {@code x} and {@code y}.
678 */
679 public static boolean equals(double x, double y, int maxUlps) {
680 // Check that "maxUlps" is non-negative and small enough so that
681 // NaN won't compare as equal to anything (except another NaN).
682 assert maxUlps > 0 && maxUlps < NAN_GAP;
683
684 long xInt = Double.doubleToLongBits(x);
685 long yInt = Double.doubleToLongBits(y);
686
687 // Make lexicographically ordered as a two's-complement integer.
688 if (xInt < 0) {
689 xInt = SGN_MASK - xInt;
690 }
691 if (yInt < 0) {
692 yInt = SGN_MASK - yInt;
693 }
694
695 return FastMath.abs(xInt - yInt) <= maxUlps;
696 }
697
698 /**
699 * Returns true if both arguments are NaN or if they are equal as defined
700 * by {@link #equals(double,double,int) equals(x, y, maxUlps}.
701 *
702 * @param x first value
703 * @param y second value
704 * @param maxUlps {@code (maxUlps - 1)} is the number of floating point
705 * values between {@code x} and {@code y}.
706 * @return {@code true} if both arguments are NaN or if there are less than
707 * {@code maxUlps} floating point values between {@code x} and {@code y}.
708 * @since 2.2
709 */
710 public static boolean equalsIncludingNaN(double x, double y, int maxUlps) {
711 return (Double.isNaN(x) && Double.isNaN(y)) || equals(x, y, maxUlps);
712 }
713
714 /**
715 * Returns true iff both arguments are null or have same dimensions and all
716 * their elements are equal as defined by
717 * {@link #equals(double,double)}.
718 *
719 * <p>This method considers that {@code NaN == NaN}. In release
720 * 3.0, the semantics will change in order to comply with IEEE754 where it
721 * is specified that {@code NaN != NaN}.
722 * New methods have been added for those cases wher the old semantics is
723 * useful (see e.g. {@link #equalsIncludingNaN(double[],double[])
724 * equalsIncludingNaN}.
725 * </p>
726 * @param x first array
727 * @param y second array
728 * @return true if the values are both null or have same dimension
729 * and equal elements.
730 */
731 public static boolean equals(double[] x, double[] y) {
732 if ((x == null) || (y == null)) {
733 return !((x == null) ^ (y == null));
734 }
735 if (x.length != y.length) {
736 return false;
737 }
738 for (int i = 0; i < x.length; ++i) {
739 if (!equals(x[i], y[i])) {
740 return false;
741 }
742 }
743 return true;
744 }
745
746 /**
747 * Returns true iff both arguments are null or have same dimensions and all
748 * their elements are equal as defined by
749 * {@link #equalsIncludingNaN(double,double)}.
750 *
751 * @param x first array
752 * @param y second array
753 * @return true if the values are both null or have same dimension and
754 * equal elements
755 * @since 2.2
756 */
757 public static boolean equalsIncludingNaN(double[] x, double[] y) {
758 if ((x == null) || (y == null)) {
759 return !((x == null) ^ (y == null));
760 }
761 if (x.length != y.length) {
762 return false;
763 }
764 for (int i = 0; i < x.length; ++i) {
765 if (!equalsIncludingNaN(x[i], y[i])) {
766 return false;
767 }
768 }
769 return true;
770 }
771
772 /**
773 * Returns n!. Shorthand for <code>n</code> <a
774 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
775 * product of the numbers <code>1,...,n</code>.
776 * <p>
777 * <Strong>Preconditions</strong>:
778 * <ul>
779 * <li> <code>n >= 0</code> (otherwise
780 * <code>IllegalArgumentException</code> is thrown)</li>
781 * <li> The result is small enough to fit into a <code>long</code>. The
782 * largest value of <code>n</code> for which <code>n!</code> <
783 * Long.MAX_VALUE</code> is 20. If the computed value exceeds <code>Long.MAX_VALUE</code>
784 * an <code>ArithMeticException </code> is thrown.</li>
785 * </ul>
786 * </p>
787 *
788 * @param n argument
789 * @return <code>n!</code>
790 * @throws ArithmeticException if the result is too large to be represented
791 * by a long integer.
792 * @throws IllegalArgumentException if n < 0
793 */
794 public static long factorial(final int n) {
795 if (n < 0) {
796 throw MathRuntimeException.createIllegalArgumentException(
797 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
798 n);
799 }
800 if (n > 20) {
801 throw new ArithmeticException(
802 "factorial value is too large to fit in a long");
803 }
804 return FACTORIALS[n];
805 }
806
807 /**
808 * Returns n!. Shorthand for <code>n</code> <a
809 * href="http://mathworld.wolfram.com/Factorial.html"> Factorial</a>, the
810 * product of the numbers <code>1,...,n</code> as a <code>double</code>.
811 * <p>
812 * <Strong>Preconditions</strong>:
813 * <ul>
814 * <li> <code>n >= 0</code> (otherwise
815 * <code>IllegalArgumentException</code> is thrown)</li>
816 * <li> The result is small enough to fit into a <code>double</code>. The
817 * largest value of <code>n</code> for which <code>n!</code> <
818 * Double.MAX_VALUE</code> is 170. If the computed value exceeds
819 * Double.MAX_VALUE, Double.POSITIVE_INFINITY is returned</li>
820 * </ul>
821 * </p>
822 *
823 * @param n argument
824 * @return <code>n!</code>
825 * @throws IllegalArgumentException if n < 0
826 */
827 public static double factorialDouble(final int n) {
828 if (n < 0) {
829 throw MathRuntimeException.createIllegalArgumentException(
830 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
831 n);
832 }
833 if (n < 21) {
834 return factorial(n);
835 }
836 return FastMath.floor(FastMath.exp(factorialLog(n)) + 0.5);
837 }
838
839 /**
840 * Returns the natural logarithm of n!.
841 * <p>
842 * <Strong>Preconditions</strong>:
843 * <ul>
844 * <li> <code>n >= 0</code> (otherwise
845 * <code>IllegalArgumentException</code> is thrown)</li>
846 * </ul></p>
847 *
848 * @param n argument
849 * @return <code>n!</code>
850 * @throws IllegalArgumentException if preconditions are not met.
851 */
852 public static double factorialLog(final int n) {
853 if (n < 0) {
854 throw MathRuntimeException.createIllegalArgumentException(
855 LocalizedFormats.FACTORIAL_NEGATIVE_PARAMETER,
856 n);
857 }
858 if (n < 21) {
859 return FastMath.log(factorial(n));
860 }
861 double logSum = 0;
862 for (int i = 2; i <= n; i++) {
863 logSum += FastMath.log(i);
864 }
865 return logSum;
866 }
867
868 /**
869 * <p>
870 * Gets the greatest common divisor of the absolute value of two numbers,
871 * using the "binary gcd" method which avoids division and modulo
872 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
873 * Stein (1961).
874 * </p>
875 * Special cases:
876 * <ul>
877 * <li>The invocations
878 * <code>gcd(Integer.MIN_VALUE, Integer.MIN_VALUE)</code>,
879 * <code>gcd(Integer.MIN_VALUE, 0)</code> and
880 * <code>gcd(0, Integer.MIN_VALUE)</code> throw an
881 * <code>ArithmeticException</code>, because the result would be 2^31, which
882 * is too large for an int value.</li>
883 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0, x)</code> and
884 * <code>gcd(x, 0)</code> is the absolute value of <code>x</code>, except
885 * for the special cases above.
886 * <li>The invocation <code>gcd(0, 0)</code> is the only one which returns
887 * <code>0</code>.</li>
888 * </ul>
889 *
890 * @param p any number
891 * @param q any number
892 * @return the greatest common divisor, never negative
893 * @throws ArithmeticException if the result cannot be represented as a
894 * nonnegative int value
895 * @since 1.1
896 */
897 public static int gcd(final int p, final int q) {
898 int u = p;
899 int v = q;
900 if ((u == 0) || (v == 0)) {
901 if ((u == Integer.MIN_VALUE) || (v == Integer.MIN_VALUE)) {
902 throw MathRuntimeException.createArithmeticException(
903 LocalizedFormats.GCD_OVERFLOW_32_BITS,
904 p, q);
905 }
906 return FastMath.abs(u) + FastMath.abs(v);
907 }
908 // keep u and v negative, as negative integers range down to
909 // -2^31, while positive numbers can only be as large as 2^31-1
910 // (i.e. we can't necessarily negate a negative number without
911 // overflow)
912 /* assert u!=0 && v!=0; */
913 if (u > 0) {
914 u = -u;
915 } // make u negative
916 if (v > 0) {
917 v = -v;
918 } // make v negative
919 // B1. [Find power of 2]
920 int k = 0;
921 while ((u & 1) == 0 && (v & 1) == 0 && k < 31) { // while u and v are
922 // both even...
923 u /= 2;
924 v /= 2;
925 k++; // cast out twos.
926 }
927 if (k == 31) {
928 throw MathRuntimeException.createArithmeticException(
929 LocalizedFormats.GCD_OVERFLOW_32_BITS,
930 p, q);
931 }
932 // B2. Initialize: u and v have been divided by 2^k and at least
933 // one is odd.
934 int t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
935 // t negative: u was odd, v may be even (t replaces v)
936 // t positive: u was even, v is odd (t replaces u)
937 do {
938 /* assert u<0 && v<0; */
939 // B4/B3: cast out twos from t.
940 while ((t & 1) == 0) { // while t is even..
941 t /= 2; // cast out twos
942 }
943 // B5 [reset max(u,v)]
944 if (t > 0) {
945 u = -t;
946 } else {
947 v = t;
948 }
949 // B6/B3. at this point both u and v should be odd.
950 t = (v - u) / 2;
951 // |u| larger: t positive (replace u)
952 // |v| larger: t negative (replace v)
953 } while (t != 0);
954 return -u * (1 << k); // gcd is u*2^k
955 }
956
957 /**
958 * <p>
959 * Gets the greatest common divisor of the absolute value of two numbers,
960 * using the "binary gcd" method which avoids division and modulo
961 * operations. See Knuth 4.5.2 algorithm B. This algorithm is due to Josef
962 * Stein (1961).
963 * </p>
964 * Special cases:
965 * <ul>
966 * <li>The invocations
967 * <code>gcd(Long.MIN_VALUE, Long.MIN_VALUE)</code>,
968 * <code>gcd(Long.MIN_VALUE, 0L)</code> and
969 * <code>gcd(0L, Long.MIN_VALUE)</code> throw an
970 * <code>ArithmeticException</code>, because the result would be 2^63, which
971 * is too large for a long value.</li>
972 * <li>The result of <code>gcd(x, x)</code>, <code>gcd(0L, x)</code> and
973 * <code>gcd(x, 0L)</code> is the absolute value of <code>x</code>, except
974 * for the special cases above.
975 * <li>The invocation <code>gcd(0L, 0L)</code> is the only one which returns
976 * <code>0L</code>.</li>
977 * </ul>
978 *
979 * @param p any number
980 * @param q any number
981 * @return the greatest common divisor, never negative
982 * @throws ArithmeticException if the result cannot be represented as a nonnegative long
983 * value
984 * @since 2.1
985 */
986 public static long gcd(final long p, final long q) {
987 long u = p;
988 long v = q;
989 if ((u == 0) || (v == 0)) {
990 if ((u == Long.MIN_VALUE) || (v == Long.MIN_VALUE)){
991 throw MathRuntimeException.createArithmeticException(
992 LocalizedFormats.GCD_OVERFLOW_64_BITS,
993 p, q);
994 }
995 return FastMath.abs(u) + FastMath.abs(v);
996 }
997 // keep u and v negative, as negative integers range down to
998 // -2^63, while positive numbers can only be as large as 2^63-1
999 // (i.e. we can't necessarily negate a negative number without
1000 // overflow)
1001 /* assert u!=0 && v!=0; */
1002 if (u > 0) {
1003 u = -u;
1004 } // make u negative
1005 if (v > 0) {
1006 v = -v;
1007 } // make v negative
1008 // B1. [Find power of 2]
1009 int k = 0;
1010 while ((u & 1) == 0 && (v & 1) == 0 && k < 63) { // while u and v are
1011 // both even...
1012 u /= 2;
1013 v /= 2;
1014 k++; // cast out twos.
1015 }
1016 if (k == 63) {
1017 throw MathRuntimeException.createArithmeticException(
1018 LocalizedFormats.GCD_OVERFLOW_64_BITS,
1019 p, q);
1020 }
1021 // B2. Initialize: u and v have been divided by 2^k and at least
1022 // one is odd.
1023 long t = ((u & 1) == 1) ? v : -(u / 2)/* B3 */;
1024 // t negative: u was odd, v may be even (t replaces v)
1025 // t positive: u was even, v is odd (t replaces u)
1026 do {
1027 /* assert u<0 && v<0; */
1028 // B4/B3: cast out twos from t.
1029 while ((t & 1) == 0) { // while t is even..
1030 t /= 2; // cast out twos
1031 }
1032 // B5 [reset max(u,v)]
1033 if (t > 0) {
1034 u = -t;
1035 } else {
1036 v = t;
1037 }
1038 // B6/B3. at this point both u and v should be odd.
1039 t = (v - u) / 2;
1040 // |u| larger: t positive (replace u)
1041 // |v| larger: t negative (replace v)
1042 } while (t != 0);
1043 return -u * (1L << k); // gcd is u*2^k
1044 }
1045
1046 /**
1047 * Returns an integer hash code representing the given double value.
1048 *
1049 * @param value the value to be hashed
1050 * @return the hash code
1051 */
1052 public static int hash(double value) {
1053 return new Double(value).hashCode();
1054 }
1055
1056 /**
1057 * Returns an integer hash code representing the given double array.
1058 *
1059 * @param value the value to be hashed (may be null)
1060 * @return the hash code
1061 * @since 1.2
1062 */
1063 public static int hash(double[] value) {
1064 return Arrays.hashCode(value);
1065 }
1066
1067 /**
1068 * For a byte value x, this method returns (byte)(+1) if x >= 0 and
1069 * (byte)(-1) if x < 0.
1070 *
1071 * @param x the value, a byte
1072 * @return (byte)(+1) or (byte)(-1), depending on the sign of x
1073 */
1074 public static byte indicator(final byte x) {
1075 return (x >= ZB) ? PB : NB;
1076 }
1077
1078 /**
1079 * For a double precision value x, this method returns +1.0 if x >= 0 and
1080 * -1.0 if x < 0. Returns <code>NaN</code> if <code>x</code> is
1081 * <code>NaN</code>.
1082 *
1083 * @param x the value, a double
1084 * @return +1.0 or -1.0, depending on the sign of x
1085 */
1086 public static double indicator(final double x) {
1087 if (Double.isNaN(x)) {
1088 return Double.NaN;
1089 }
1090 return (x >= 0.0) ? 1.0 : -1.0;
1091 }
1092
1093 /**
1094 * For a float value x, this method returns +1.0F if x >= 0 and -1.0F if x <
1095 * 0. Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.
1096 *
1097 * @param x the value, a float
1098 * @return +1.0F or -1.0F, depending on the sign of x
1099 */
1100 public static float indicator(final float x) {
1101 if (Float.isNaN(x)) {
1102 return Float.NaN;
1103 }
1104 return (x >= 0.0F) ? 1.0F : -1.0F;
1105 }
1106
1107 /**
1108 * For an int value x, this method returns +1 if x >= 0 and -1 if x < 0.
1109 *
1110 * @param x the value, an int
1111 * @return +1 or -1, depending on the sign of x
1112 */
1113 public static int indicator(final int x) {
1114 return (x >= 0) ? 1 : -1;
1115 }
1116
1117 /**
1118 * For a long value x, this method returns +1L if x >= 0 and -1L if x < 0.
1119 *
1120 * @param x the value, a long
1121 * @return +1L or -1L, depending on the sign of x
1122 */
1123 public static long indicator(final long x) {
1124 return (x >= 0L) ? 1L : -1L;
1125 }
1126
1127 /**
1128 * For a short value x, this method returns (short)(+1) if x >= 0 and
1129 * (short)(-1) if x < 0.
1130 *
1131 * @param x the value, a short
1132 * @return (short)(+1) or (short)(-1), depending on the sign of x
1133 */
1134 public static short indicator(final short x) {
1135 return (x >= ZS) ? PS : NS;
1136 }
1137
1138 /**
1139 * <p>
1140 * Returns the least common multiple of the absolute value of two numbers,
1141 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
1142 * </p>
1143 * Special cases:
1144 * <ul>
1145 * <li>The invocations <code>lcm(Integer.MIN_VALUE, n)</code> and
1146 * <code>lcm(n, Integer.MIN_VALUE)</code>, where <code>abs(n)</code> is a
1147 * power of 2, throw an <code>ArithmeticException</code>, because the result
1148 * would be 2^31, which is too large for an int value.</li>
1149 * <li>The result of <code>lcm(0, x)</code> and <code>lcm(x, 0)</code> is
1150 * <code>0</code> for any <code>x</code>.
1151 * </ul>
1152 *
1153 * @param a any number
1154 * @param b any number
1155 * @return the least common multiple, never negative
1156 * @throws ArithmeticException
1157 * if the result cannot be represented as a nonnegative int
1158 * value
1159 * @since 1.1
1160 */
1161 public static int lcm(int a, int b) {
1162 if (a==0 || b==0){
1163 return 0;
1164 }
1165 int lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b));
1166 if (lcm == Integer.MIN_VALUE) {
1167 throw MathRuntimeException.createArithmeticException(
1168 LocalizedFormats.LCM_OVERFLOW_32_BITS,
1169 a, b);
1170 }
1171 return lcm;
1172 }
1173
1174 /**
1175 * <p>
1176 * Returns the least common multiple of the absolute value of two numbers,
1177 * using the formula <code>lcm(a,b) = (a / gcd(a,b)) * b</code>.
1178 * </p>
1179 * Special cases:
1180 * <ul>
1181 * <li>The invocations <code>lcm(Long.MIN_VALUE, n)</code> and
1182 * <code>lcm(n, Long.MIN_VALUE)</code>, where <code>abs(n)</code> is a
1183 * power of 2, throw an <code>ArithmeticException</code>, because the result
1184 * would be 2^63, which is too large for an int value.</li>
1185 * <li>The result of <code>lcm(0L, x)</code> and <code>lcm(x, 0L)</code> is
1186 * <code>0L</code> for any <code>x</code>.
1187 * </ul>
1188 *
1189 * @param a any number
1190 * @param b any number
1191 * @return the least common multiple, never negative
1192 * @throws ArithmeticException if the result cannot be represented as a nonnegative long
1193 * value
1194 * @since 2.1
1195 */
1196 public static long lcm(long a, long b) {
1197 if (a==0 || b==0){
1198 return 0;
1199 }
1200 long lcm = FastMath.abs(mulAndCheck(a / gcd(a, b), b));
1201 if (lcm == Long.MIN_VALUE){
1202 throw MathRuntimeException.createArithmeticException(
1203 LocalizedFormats.LCM_OVERFLOW_64_BITS,
1204 a, b);
1205 }
1206 return lcm;
1207 }
1208
1209 /**
1210 * <p>Returns the
1211 * <a href="http://mathworld.wolfram.com/Logarithm.html">logarithm</a>
1212 * for base <code>b</code> of <code>x</code>.
1213 * </p>
1214 * <p>Returns <code>NaN<code> if either argument is negative. If
1215 * <code>base</code> is 0 and <code>x</code> is positive, 0 is returned.
1216 * If <code>base</code> is positive and <code>x</code> is 0,
1217 * <code>Double.NEGATIVE_INFINITY</code> is returned. If both arguments
1218 * are 0, the result is <code>NaN</code>.</p>
1219 *
1220 * @param base the base of the logarithm, must be greater than 0
1221 * @param x argument, must be greater than 0
1222 * @return the value of the logarithm - the number y such that base^y = x.
1223 * @since 1.2
1224 */
1225 public static double log(double base, double x) {
1226 return FastMath.log(x)/FastMath.log(base);
1227 }
1228
1229 /**
1230 * Multiply two integers, checking for overflow.
1231 *
1232 * @param x a factor
1233 * @param y a factor
1234 * @return the product <code>x*y</code>
1235 * @throws ArithmeticException if the result can not be represented as an
1236 * int
1237 * @since 1.1
1238 */
1239 public static int mulAndCheck(int x, int y) {
1240 long m = ((long)x) * ((long)y);
1241 if (m < Integer.MIN_VALUE || m > Integer.MAX_VALUE) {
1242 throw new ArithmeticException("overflow: mul");
1243 }
1244 return (int)m;
1245 }
1246
1247 /**
1248 * Multiply two long integers, checking for overflow.
1249 *
1250 * @param a first value
1251 * @param b second value
1252 * @return the product <code>a * b</code>
1253 * @throws ArithmeticException if the result can not be represented as an
1254 * long
1255 * @since 1.2
1256 */
1257 public static long mulAndCheck(long a, long b) {
1258 long ret;
1259 String msg = "overflow: multiply";
1260 if (a > b) {
1261 // use symmetry to reduce boundary cases
1262 ret = mulAndCheck(b, a);
1263 } else {
1264 if (a < 0) {
1265 if (b < 0) {
1266 // check for positive overflow with negative a, negative b
1267 if (a >= Long.MAX_VALUE / b) {
1268 ret = a * b;
1269 } else {
1270 throw new ArithmeticException(msg);
1271 }
1272 } else if (b > 0) {
1273 // check for negative overflow with negative a, positive b
1274 if (Long.MIN_VALUE / b <= a) {
1275 ret = a * b;
1276 } else {
1277 throw new ArithmeticException(msg);
1278
1279 }
1280 } else {
1281 // assert b == 0
1282 ret = 0;
1283 }
1284 } else if (a > 0) {
1285 // assert a > 0
1286 // assert b > 0
1287
1288 // check for positive overflow with positive a, positive b
1289 if (a <= Long.MAX_VALUE / b) {
1290 ret = a * b;
1291 } else {
1292 throw new ArithmeticException(msg);
1293 }
1294 } else {
1295 // assert a == 0
1296 ret = 0;
1297 }
1298 }
1299 return ret;
1300 }
1301
1302 /**
1303 * Get the next machine representable number after a number, moving
1304 * in the direction of another number.
1305 * <p>
1306 * If <code>direction</code> is greater than or equal to<code>d</code>,
1307 * the smallest machine representable number strictly greater than
1308 * <code>d</code> is returned; otherwise the largest representable number
1309 * strictly less than <code>d</code> is returned.</p>
1310 * <p>
1311 * If <code>d</code> is NaN or Infinite, it is returned unchanged.</p>
1312 *
1313 * @param d base number
1314 * @param direction (the only important thing is whether
1315 * direction is greater or smaller than d)
1316 * @return the next machine representable number in the specified direction
1317 * @since 1.2
1318 * @deprecated as of 2.2, replaced by {@link FastMath#nextAfter(double, double)}
1319 * which handles Infinities differently, and returns direction if d and direction compare equal.
1320 */
1321 @Deprecated
1322 public static double nextAfter(double d, double direction) {
1323
1324 // handling of some important special cases
1325 if (Double.isNaN(d) || Double.isInfinite(d)) {
1326 return d;
1327 } else if (d == 0) {
1328 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
1329 }
1330 // special cases MAX_VALUE to infinity and MIN_VALUE to 0
1331 // are handled just as normal numbers
1332
1333 // split the double in raw components
1334 long bits = Double.doubleToLongBits(d);
1335 long sign = bits & 0x8000000000000000L;
1336 long exponent = bits & 0x7ff0000000000000L;
1337 long mantissa = bits & 0x000fffffffffffffL;
1338
1339 if (d * (direction - d) >= 0) {
1340 // we should increase the mantissa
1341 if (mantissa == 0x000fffffffffffffL) {
1342 return Double.longBitsToDouble(sign |
1343 (exponent + 0x0010000000000000L));
1344 } else {
1345 return Double.longBitsToDouble(sign |
1346 exponent | (mantissa + 1));
1347 }
1348 } else {
1349 // we should decrease the mantissa
1350 if (mantissa == 0L) {
1351 return Double.longBitsToDouble(sign |
1352 (exponent - 0x0010000000000000L) |
1353 0x000fffffffffffffL);
1354 } else {
1355 return Double.longBitsToDouble(sign |
1356 exponent | (mantissa - 1));
1357 }
1358 }
1359 }
1360
1361 /**
1362 * Scale a number by 2<sup>scaleFactor</sup>.
1363 * <p>If <code>d</code> is 0 or NaN or Infinite, it is returned unchanged.</p>
1364 *
1365 * @param d base number
1366 * @param scaleFactor power of two by which d should be multiplied
1367 * @return d × 2<sup>scaleFactor</sup>
1368 * @since 2.0
1369 * @deprecated as of 2.2, replaced by {@link FastMath#scalb(double, int)}
1370 */
1371 @Deprecated
1372 public static double scalb(final double d, final int scaleFactor) {
1373 return FastMath.scalb(d, scaleFactor);
1374 }
1375
1376 /**
1377 * Normalize an angle in a 2&pi wide interval around a center value.
1378 * <p>This method has three main uses:</p>
1379 * <ul>
1380 * <li>normalize an angle between 0 and 2π:<br/>
1381 * <code>a = MathUtils.normalizeAngle(a, FastMath.PI);</code></li>
1382 * <li>normalize an angle between -π and +π<br/>
1383 * <code>a = MathUtils.normalizeAngle(a, 0.0);</code></li>
1384 * <li>compute the angle between two defining angular positions:<br>
1385 * <code>angle = MathUtils.normalizeAngle(end, start) - start;</code></li>
1386 * </ul>
1387 * <p>Note that due to numerical accuracy and since π cannot be represented
1388 * exactly, the result interval is <em>closed</em>, it cannot be half-closed
1389 * as would be more satisfactory in a purely mathematical view.</p>
1390 * @param a angle to normalize
1391 * @param center center of the desired 2π interval for the result
1392 * @return a-2kπ with integer k and center-π <= a-2kπ <= center+π
1393 * @since 1.2
1394 */
1395 public static double normalizeAngle(double a, double center) {
1396 return a - TWO_PI * FastMath.floor((a + FastMath.PI - center) / TWO_PI);
1397 }
1398
1399 /**
1400 * <p>Normalizes an array to make it sum to a specified value.
1401 * Returns the result of the transformation <pre>
1402 * x |-> x * normalizedSum / sum
1403 * </pre>
1404 * applied to each non-NaN element x of the input array, where sum is the
1405 * sum of the non-NaN entries in the input array.</p>
1406 *
1407 * <p>Throws IllegalArgumentException if <code>normalizedSum</code> is infinite
1408 * or NaN and ArithmeticException if the input array contains any infinite elements
1409 * or sums to 0</p>
1410 *
1411 * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
1412 *
1413 * @param values input array to be normalized
1414 * @param normalizedSum target sum for the normalized array
1415 * @return normalized array
1416 * @throws ArithmeticException if the input array contains infinite elements or sums to zero
1417 * @throws IllegalArgumentException if the target sum is infinite or NaN
1418 * @since 2.1
1419 */
1420 public static double[] normalizeArray(double[] values, double normalizedSum)
1421 throws ArithmeticException, IllegalArgumentException {
1422 if (Double.isInfinite(normalizedSum)) {
1423 throw MathRuntimeException.createIllegalArgumentException(
1424 LocalizedFormats.NORMALIZE_INFINITE);
1425 }
1426 if (Double.isNaN(normalizedSum)) {
1427 throw MathRuntimeException.createIllegalArgumentException(
1428 LocalizedFormats.NORMALIZE_NAN);
1429 }
1430 double sum = 0d;
1431 final int len = values.length;
1432 double[] out = new double[len];
1433 for (int i = 0; i < len; i++) {
1434 if (Double.isInfinite(values[i])) {
1435 throw MathRuntimeException.createArithmeticException(
1436 LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
1437 }
1438 if (!Double.isNaN(values[i])) {
1439 sum += values[i];
1440 }
1441 }
1442 if (sum == 0) {
1443 throw MathRuntimeException.createArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
1444 }
1445 for (int i = 0; i < len; i++) {
1446 if (Double.isNaN(values[i])) {
1447 out[i] = Double.NaN;
1448 } else {
1449 out[i] = values[i] * normalizedSum / sum;
1450 }
1451 }
1452 return out;
1453 }
1454
1455 /**
1456 * Round the given value to the specified number of decimal places. The
1457 * value is rounded using the {@link BigDecimal#ROUND_HALF_UP} method.
1458 *
1459 * @param x the value to round.
1460 * @param scale the number of digits to the right of the decimal point.
1461 * @return the rounded value.
1462 * @since 1.1
1463 */
1464 public static double round(double x, int scale) {
1465 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1466 }
1467
1468 /**
1469 * Round the given value to the specified number of decimal places. The
1470 * value is rounded using the given method which is any method defined in
1471 * {@link BigDecimal}.
1472 *
1473 * @param x the value to round.
1474 * @param scale the number of digits to the right of the decimal point.
1475 * @param roundingMethod the rounding method as defined in
1476 * {@link BigDecimal}.
1477 * @return the rounded value.
1478 * @since 1.1
1479 */
1480 public static double round(double x, int scale, int roundingMethod) {
1481 try {
1482 return (new BigDecimal
1483 (Double.toString(x))
1484 .setScale(scale, roundingMethod))
1485 .doubleValue();
1486 } catch (NumberFormatException ex) {
1487 if (Double.isInfinite(x)) {
1488 return x;
1489 } else {
1490 return Double.NaN;
1491 }
1492 }
1493 }
1494
1495 /**
1496 * Round the given value to the specified number of decimal places. The
1497 * value is rounding using the {@link BigDecimal#ROUND_HALF_UP} method.
1498 *
1499 * @param x the value to round.
1500 * @param scale the number of digits to the right of the decimal point.
1501 * @return the rounded value.
1502 * @since 1.1
1503 */
1504 public static float round(float x, int scale) {
1505 return round(x, scale, BigDecimal.ROUND_HALF_UP);
1506 }
1507
1508 /**
1509 * Round the given value to the specified number of decimal places. The
1510 * value is rounded using the given method which is any method defined in
1511 * {@link BigDecimal}.
1512 *
1513 * @param x the value to round.
1514 * @param scale the number of digits to the right of the decimal point.
1515 * @param roundingMethod the rounding method as defined in
1516 * {@link BigDecimal}.
1517 * @return the rounded value.
1518 * @since 1.1
1519 */
1520 public static float round(float x, int scale, int roundingMethod) {
1521 float sign = indicator(x);
1522 float factor = (float)FastMath.pow(10.0f, scale) * sign;
1523 return (float)roundUnscaled(x * factor, sign, roundingMethod) / factor;
1524 }
1525
1526 /**
1527 * Round the given non-negative, value to the "nearest" integer. Nearest is
1528 * determined by the rounding method specified. Rounding methods are defined
1529 * in {@link BigDecimal}.
1530 *
1531 * @param unscaled the value to round.
1532 * @param sign the sign of the original, scaled value.
1533 * @param roundingMethod the rounding method as defined in
1534 * {@link BigDecimal}.
1535 * @return the rounded value.
1536 * @since 1.1
1537 */
1538 private static double roundUnscaled(double unscaled, double sign,
1539 int roundingMethod) {
1540 switch (roundingMethod) {
1541 case BigDecimal.ROUND_CEILING :
1542 if (sign == -1) {
1543 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1544 } else {
1545 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1546 }
1547 break;
1548 case BigDecimal.ROUND_DOWN :
1549 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1550 break;
1551 case BigDecimal.ROUND_FLOOR :
1552 if (sign == -1) {
1553 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1554 } else {
1555 unscaled = FastMath.floor(nextAfter(unscaled, Double.NEGATIVE_INFINITY));
1556 }
1557 break;
1558 case BigDecimal.ROUND_HALF_DOWN : {
1559 unscaled = nextAfter(unscaled, Double.NEGATIVE_INFINITY);
1560 double fraction = unscaled - FastMath.floor(unscaled);
1561 if (fraction > 0.5) {
1562 unscaled = FastMath.ceil(unscaled);
1563 } else {
1564 unscaled = FastMath.floor(unscaled);
1565 }
1566 break;
1567 }
1568 case BigDecimal.ROUND_HALF_EVEN : {
1569 double fraction = unscaled - FastMath.floor(unscaled);
1570 if (fraction > 0.5) {
1571 unscaled = FastMath.ceil(unscaled);
1572 } else if (fraction < 0.5) {
1573 unscaled = FastMath.floor(unscaled);
1574 } else {
1575 // The following equality test is intentional and needed for rounding purposes
1576 if (FastMath.floor(unscaled) / 2.0 == FastMath.floor(Math
1577 .floor(unscaled) / 2.0)) { // even
1578 unscaled = FastMath.floor(unscaled);
1579 } else { // odd
1580 unscaled = FastMath.ceil(unscaled);
1581 }
1582 }
1583 break;
1584 }
1585 case BigDecimal.ROUND_HALF_UP : {
1586 unscaled = nextAfter(unscaled, Double.POSITIVE_INFINITY);
1587 double fraction = unscaled - FastMath.floor(unscaled);
1588 if (fraction >= 0.5) {
1589 unscaled = FastMath.ceil(unscaled);
1590 } else {
1591 unscaled = FastMath.floor(unscaled);
1592 }
1593 break;
1594 }
1595 case BigDecimal.ROUND_UNNECESSARY :
1596 if (unscaled != FastMath.floor(unscaled)) {
1597 throw new ArithmeticException("Inexact result from rounding");
1598 }
1599 break;
1600 case BigDecimal.ROUND_UP :
1601 unscaled = FastMath.ceil(nextAfter(unscaled, Double.POSITIVE_INFINITY));
1602 break;
1603 default :
1604 throw MathRuntimeException.createIllegalArgumentException(
1605 LocalizedFormats.INVALID_ROUNDING_METHOD,
1606 roundingMethod,
1607 "ROUND_CEILING", BigDecimal.ROUND_CEILING,
1608 "ROUND_DOWN", BigDecimal.ROUND_DOWN,
1609 "ROUND_FLOOR", BigDecimal.ROUND_FLOOR,
1610 "ROUND_HALF_DOWN", BigDecimal.ROUND_HALF_DOWN,
1611 "ROUND_HALF_EVEN", BigDecimal.ROUND_HALF_EVEN,
1612 "ROUND_HALF_UP", BigDecimal.ROUND_HALF_UP,
1613 "ROUND_UNNECESSARY", BigDecimal.ROUND_UNNECESSARY,
1614 "ROUND_UP", BigDecimal.ROUND_UP);
1615 }
1616 return unscaled;
1617 }
1618
1619 /**
1620 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1621 * for byte value <code>x</code>.
1622 * <p>
1623 * For a byte value x, this method returns (byte)(+1) if x > 0, (byte)(0) if
1624 * x = 0, and (byte)(-1) if x < 0.</p>
1625 *
1626 * @param x the value, a byte
1627 * @return (byte)(+1), (byte)(0), or (byte)(-1), depending on the sign of x
1628 */
1629 public static byte sign(final byte x) {
1630 return (x == ZB) ? ZB : (x > ZB) ? PB : NB;
1631 }
1632
1633 /**
1634 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1635 * for double precision <code>x</code>.
1636 * <p>
1637 * For a double value <code>x</code>, this method returns
1638 * <code>+1.0</code> if <code>x > 0</code>, <code>0.0</code> if
1639 * <code>x = 0.0</code>, and <code>-1.0</code> if <code>x < 0</code>.
1640 * Returns <code>NaN</code> if <code>x</code> is <code>NaN</code>.</p>
1641 *
1642 * @param x the value, a double
1643 * @return +1.0, 0.0, or -1.0, depending on the sign of x
1644 */
1645 public static double sign(final double x) {
1646 if (Double.isNaN(x)) {
1647 return Double.NaN;
1648 }
1649 return (x == 0.0) ? 0.0 : (x > 0.0) ? 1.0 : -1.0;
1650 }
1651
1652 /**
1653 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1654 * for float value <code>x</code>.
1655 * <p>
1656 * For a float value x, this method returns +1.0F if x > 0, 0.0F if x =
1657 * 0.0F, and -1.0F if x < 0. Returns <code>NaN</code> if <code>x</code>
1658 * is <code>NaN</code>.</p>
1659 *
1660 * @param x the value, a float
1661 * @return +1.0F, 0.0F, or -1.0F, depending on the sign of x
1662 */
1663 public static float sign(final float x) {
1664 if (Float.isNaN(x)) {
1665 return Float.NaN;
1666 }
1667 return (x == 0.0F) ? 0.0F : (x > 0.0F) ? 1.0F : -1.0F;
1668 }
1669
1670 /**
1671 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1672 * for int value <code>x</code>.
1673 * <p>
1674 * For an int value x, this method returns +1 if x > 0, 0 if x = 0, and -1
1675 * if x < 0.</p>
1676 *
1677 * @param x the value, an int
1678 * @return +1, 0, or -1, depending on the sign of x
1679 */
1680 public static int sign(final int x) {
1681 return (x == 0) ? 0 : (x > 0) ? 1 : -1;
1682 }
1683
1684 /**
1685 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1686 * for long value <code>x</code>.
1687 * <p>
1688 * For a long value x, this method returns +1L if x > 0, 0L if x = 0, and
1689 * -1L if x < 0.</p>
1690 *
1691 * @param x the value, a long
1692 * @return +1L, 0L, or -1L, depending on the sign of x
1693 */
1694 public static long sign(final long x) {
1695 return (x == 0L) ? 0L : (x > 0L) ? 1L : -1L;
1696 }
1697
1698 /**
1699 * Returns the <a href="http://mathworld.wolfram.com/Sign.html"> sign</a>
1700 * for short value <code>x</code>.
1701 * <p>
1702 * For a short value x, this method returns (short)(+1) if x > 0, (short)(0)
1703 * if x = 0, and (short)(-1) if x < 0.</p>
1704 *
1705 * @param x the value, a short
1706 * @return (short)(+1), (short)(0), or (short)(-1), depending on the sign of
1707 * x
1708 */
1709 public static short sign(final short x) {
1710 return (x == ZS) ? ZS : (x > ZS) ? PS : NS;
1711 }
1712
1713 /**
1714 * Returns the <a href="http://mathworld.wolfram.com/HyperbolicSine.html">
1715 * hyperbolic sine</a> of x.
1716 *
1717 * @param x double value for which to find the hyperbolic sine
1718 * @return hyperbolic sine of x
1719 */
1720 public static double sinh(double x) {
1721 return (FastMath.exp(x) - FastMath.exp(-x)) / 2.0;
1722 }
1723
1724 /**
1725 * Subtract two integers, checking for overflow.
1726 *
1727 * @param x the minuend
1728 * @param y the subtrahend
1729 * @return the difference <code>x-y</code>
1730 * @throws ArithmeticException if the result can not be represented as an
1731 * int
1732 * @since 1.1
1733 */
1734 public static int subAndCheck(int x, int y) {
1735 long s = (long)x - (long)y;
1736 if (s < Integer.MIN_VALUE || s > Integer.MAX_VALUE) {
1737 throw MathRuntimeException.createArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, x, y);
1738 }
1739 return (int)s;
1740 }
1741
1742 /**
1743 * Subtract two long integers, checking for overflow.
1744 *
1745 * @param a first value
1746 * @param b second value
1747 * @return the difference <code>a-b</code>
1748 * @throws ArithmeticException if the result can not be represented as an
1749 * long
1750 * @since 1.2
1751 */
1752 public static long subAndCheck(long a, long b) {
1753 long ret;
1754 String msg = "overflow: subtract";
1755 if (b == Long.MIN_VALUE) {
1756 if (a < 0) {
1757 ret = a - b;
1758 } else {
1759 throw new ArithmeticException(msg);
1760 }
1761 } else {
1762 // use additive inverse
1763 ret = addAndCheck(a, -b, LocalizedFormats.OVERFLOW_IN_ADDITION);
1764 }
1765 return ret;
1766 }
1767
1768 /**
1769 * Raise an int to an int power.
1770 * @param k number to raise
1771 * @param e exponent (must be positive or null)
1772 * @return k<sup>e</sup>
1773 * @exception IllegalArgumentException if e is negative
1774 */
1775 public static int pow(final int k, int e)
1776 throws IllegalArgumentException {
1777
1778 if (e < 0) {
1779 throw MathRuntimeException.createIllegalArgumentException(
1780 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1781 k, e);
1782 }
1783
1784 int result = 1;
1785 int k2p = k;
1786 while (e != 0) {
1787 if ((e & 0x1) != 0) {
1788 result *= k2p;
1789 }
1790 k2p *= k2p;
1791 e = e >> 1;
1792 }
1793
1794 return result;
1795
1796 }
1797
1798 /**
1799 * Raise an int to a long power.
1800 * @param k number to raise
1801 * @param e exponent (must be positive or null)
1802 * @return k<sup>e</sup>
1803 * @exception IllegalArgumentException if e is negative
1804 */
1805 public static int pow(final int k, long e)
1806 throws IllegalArgumentException {
1807
1808 if (e < 0) {
1809 throw MathRuntimeException.createIllegalArgumentException(
1810 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1811 k, e);
1812 }
1813
1814 int result = 1;
1815 int k2p = k;
1816 while (e != 0) {
1817 if ((e & 0x1) != 0) {
1818 result *= k2p;
1819 }
1820 k2p *= k2p;
1821 e = e >> 1;
1822 }
1823
1824 return result;
1825
1826 }
1827
1828 /**
1829 * Raise a long to an int power.
1830 * @param k number to raise
1831 * @param e exponent (must be positive or null)
1832 * @return k<sup>e</sup>
1833 * @exception IllegalArgumentException if e is negative
1834 */
1835 public static long pow(final long k, int e)
1836 throws IllegalArgumentException {
1837
1838 if (e < 0) {
1839 throw MathRuntimeException.createIllegalArgumentException(
1840 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1841 k, e);
1842 }
1843
1844 long result = 1l;
1845 long k2p = k;
1846 while (e != 0) {
1847 if ((e & 0x1) != 0) {
1848 result *= k2p;
1849 }
1850 k2p *= k2p;
1851 e = e >> 1;
1852 }
1853
1854 return result;
1855
1856 }
1857
1858 /**
1859 * Raise a long to a long power.
1860 * @param k number to raise
1861 * @param e exponent (must be positive or null)
1862 * @return k<sup>e</sup>
1863 * @exception IllegalArgumentException if e is negative
1864 */
1865 public static long pow(final long k, long e)
1866 throws IllegalArgumentException {
1867
1868 if (e < 0) {
1869 throw MathRuntimeException.createIllegalArgumentException(
1870 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1871 k, e);
1872 }
1873
1874 long result = 1l;
1875 long k2p = k;
1876 while (e != 0) {
1877 if ((e & 0x1) != 0) {
1878 result *= k2p;
1879 }
1880 k2p *= k2p;
1881 e = e >> 1;
1882 }
1883
1884 return result;
1885
1886 }
1887
1888 /**
1889 * Raise a BigInteger to an int power.
1890 * @param k number to raise
1891 * @param e exponent (must be positive or null)
1892 * @return k<sup>e</sup>
1893 * @exception IllegalArgumentException if e is negative
1894 */
1895 public static BigInteger pow(final BigInteger k, int e)
1896 throws IllegalArgumentException {
1897
1898 if (e < 0) {
1899 throw MathRuntimeException.createIllegalArgumentException(
1900 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1901 k, e);
1902 }
1903
1904 return k.pow(e);
1905
1906 }
1907
1908 /**
1909 * Raise a BigInteger to a long power.
1910 * @param k number to raise
1911 * @param e exponent (must be positive or null)
1912 * @return k<sup>e</sup>
1913 * @exception IllegalArgumentException if e is negative
1914 */
1915 public static BigInteger pow(final BigInteger k, long e)
1916 throws IllegalArgumentException {
1917
1918 if (e < 0) {
1919 throw MathRuntimeException.createIllegalArgumentException(
1920 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1921 k, e);
1922 }
1923
1924 BigInteger result = BigInteger.ONE;
1925 BigInteger k2p = k;
1926 while (e != 0) {
1927 if ((e & 0x1) != 0) {
1928 result = result.multiply(k2p);
1929 }
1930 k2p = k2p.multiply(k2p);
1931 e = e >> 1;
1932 }
1933
1934 return result;
1935
1936 }
1937
1938 /**
1939 * Raise a BigInteger to a BigInteger power.
1940 * @param k number to raise
1941 * @param e exponent (must be positive or null)
1942 * @return k<sup>e</sup>
1943 * @exception IllegalArgumentException if e is negative
1944 */
1945 public static BigInteger pow(final BigInteger k, BigInteger e)
1946 throws IllegalArgumentException {
1947
1948 if (e.compareTo(BigInteger.ZERO) < 0) {
1949 throw MathRuntimeException.createIllegalArgumentException(
1950 LocalizedFormats.POWER_NEGATIVE_PARAMETERS,
1951 k, e);
1952 }
1953
1954 BigInteger result = BigInteger.ONE;
1955 BigInteger k2p = k;
1956 while (!BigInteger.ZERO.equals(e)) {
1957 if (e.testBit(0)) {
1958 result = result.multiply(k2p);
1959 }
1960 k2p = k2p.multiply(k2p);
1961 e = e.shiftRight(1);
1962 }
1963
1964 return result;
1965
1966 }
1967
1968 /**
1969 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1970 *
1971 * @param p1 the first point
1972 * @param p2 the second point
1973 * @return the L<sub>1</sub> distance between the two points
1974 */
1975 public static double distance1(double[] p1, double[] p2) {
1976 double sum = 0;
1977 for (int i = 0; i < p1.length; i++) {
1978 sum += FastMath.abs(p1[i] - p2[i]);
1979 }
1980 return sum;
1981 }
1982
1983 /**
1984 * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
1985 *
1986 * @param p1 the first point
1987 * @param p2 the second point
1988 * @return the L<sub>1</sub> distance between the two points
1989 */
1990 public static int distance1(int[] p1, int[] p2) {
1991 int sum = 0;
1992 for (int i = 0; i < p1.length; i++) {
1993 sum += FastMath.abs(p1[i] - p2[i]);
1994 }
1995 return sum;
1996 }
1997
1998 /**
1999 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
2000 *
2001 * @param p1 the first point
2002 * @param p2 the second point
2003 * @return the L<sub>2</sub> distance between the two points
2004 */
2005 public static double distance(double[] p1, double[] p2) {
2006 double sum = 0;
2007 for (int i = 0; i < p1.length; i++) {
2008 final double dp = p1[i] - p2[i];
2009 sum += dp * dp;
2010 }
2011 return FastMath.sqrt(sum);
2012 }
2013
2014 /**
2015 * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
2016 *
2017 * @param p1 the first point
2018 * @param p2 the second point
2019 * @return the L<sub>2</sub> distance between the two points
2020 */
2021 public static double distance(int[] p1, int[] p2) {
2022 double sum = 0;
2023 for (int i = 0; i < p1.length; i++) {
2024 final double dp = p1[i] - p2[i];
2025 sum += dp * dp;
2026 }
2027 return FastMath.sqrt(sum);
2028 }
2029
2030 /**
2031 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
2032 *
2033 * @param p1 the first point
2034 * @param p2 the second point
2035 * @return the L<sub>∞</sub> distance between the two points
2036 */
2037 public static double distanceInf(double[] p1, double[] p2) {
2038 double max = 0;
2039 for (int i = 0; i < p1.length; i++) {
2040 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
2041 }
2042 return max;
2043 }
2044
2045 /**
2046 * Calculates the L<sub>∞</sub> (max of abs) distance between two points.
2047 *
2048 * @param p1 the first point
2049 * @param p2 the second point
2050 * @return the L<sub>∞</sub> distance between the two points
2051 */
2052 public static int distanceInf(int[] p1, int[] p2) {
2053 int max = 0;
2054 for (int i = 0; i < p1.length; i++) {
2055 max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
2056 }
2057 return max;
2058 }
2059
2060 /**
2061 * Specification of ordering direction.
2062 */
2063 public static enum OrderDirection {
2064 /** Constant for increasing direction. */
2065 INCREASING,
2066 /** Constant for decreasing direction. */
2067 DECREASING
2068 }
2069
2070 /**
2071 * Checks that the given array is sorted.
2072 *
2073 * @param val Values.
2074 * @param dir Ordering direction.
2075 * @param strict Whether the order should be strict.
2076 * @throws NonMonotonousSequenceException if the array is not sorted.
2077 * @since 2.2
2078 */
2079 public static void checkOrder(double[] val, OrderDirection dir, boolean strict) {
2080 double previous = val[0];
2081 boolean ok = true;
2082
2083 int max = val.length;
2084 for (int i = 1; i < max; i++) {
2085 switch (dir) {
2086 case INCREASING:
2087 if (strict) {
2088 if (val[i] <= previous) {
2089 ok = false;
2090 }
2091 } else {
2092 if (val[i] < previous) {
2093 ok = false;
2094 }
2095 }
2096 break;
2097 case DECREASING:
2098 if (strict) {
2099 if (val[i] >= previous) {
2100 ok = false;
2101 }
2102 } else {
2103 if (val[i] > previous) {
2104 ok = false;
2105 }
2106 }
2107 break;
2108 default:
2109 // Should never happen.
2110 throw new IllegalArgumentException();
2111 }
2112
2113 if (!ok) {
2114 throw new NonMonotonousSequenceException(val[i], previous, i, dir, strict);
2115 }
2116 previous = val[i];
2117 }
2118 }
2119
2120 /**
2121 * Checks that the given array is sorted in strictly increasing order.
2122 *
2123 * @param val Values.
2124 * @throws NonMonotonousSequenceException if the array is not sorted.
2125 * @since 2.2
2126 */
2127 public static void checkOrder(double[] val) {
2128 checkOrder(val, OrderDirection.INCREASING, true);
2129 }
2130
2131 /**
2132 * Checks that the given array is sorted.
2133 *
2134 * @param val Values
2135 * @param dir Order direction (-1 for decreasing, 1 for increasing)
2136 * @param strict Whether the order should be strict
2137 * @throws NonMonotonousSequenceException if the array is not sorted.
2138 * @deprecated as of 2.2 (please use the new {@link #checkOrder(double[],OrderDirection,boolean)
2139 * checkOrder} method). To be removed in 3.0.
2140 */
2141 @Deprecated
2142 public static void checkOrder(double[] val, int dir, boolean strict) {
2143 if (dir > 0) {
2144 checkOrder(val, OrderDirection.INCREASING, strict);
2145 } else {
2146 checkOrder(val, OrderDirection.DECREASING, strict);
2147 }
2148 }
2149
2150 /**
2151 * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
2152 * Translation of the minpack enorm subroutine.
2153 *
2154 * The redistribution policy for MINPACK is available <a
2155 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it
2156 * is reproduced below.</p>
2157 *
2158 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
2159 * <tr><td>
2160 * Minpack Copyright Notice (1999) University of Chicago.
2161 * All rights reserved
2162 * </td></tr>
2163 * <tr><td>
2164 * Redistribution and use in source and binary forms, with or without
2165 * modification, are permitted provided that the following conditions
2166 * are met:
2167 * <ol>
2168 * <li>Redistributions of source code must retain the above copyright
2169 * notice, this list of conditions and the following disclaimer.</li>
2170 * <li>Redistributions in binary form must reproduce the above
2171 * copyright notice, this list of conditions and the following
2172 * disclaimer in the documentation and/or other materials provided
2173 * with the distribution.</li>
2174 * <li>The end-user documentation included with the redistribution, if any,
2175 * must include the following acknowledgment:
2176 * <code>This product includes software developed by the University of
2177 * Chicago, as Operator of Argonne National Laboratory.</code>
2178 * Alternately, this acknowledgment may appear in the software itself,
2179 * if and wherever such third-party acknowledgments normally appear.</li>
2180 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
2181 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
2182 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
2183 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
2184 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
2185 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
2186 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
2187 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
2188 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
2189 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
2190 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
2191 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
2192 * BE CORRECTED.</strong></li>
2193 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
2194 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
2195 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
2196 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
2197 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
2198 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
2199 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
2200 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
2201 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
2202 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
2203 * <ol></td></tr>
2204 * </table>
2205 *
2206 * @param v vector of doubles
2207 * @return the 2-norm of the vector
2208 * @since 2.2
2209 */
2210 public static double safeNorm(double[] v) {
2211 double rdwarf = 3.834e-20;
2212 double rgiant = 1.304e+19;
2213 double s1=0.0;
2214 double s2=0.0;
2215 double s3=0.0;
2216 double x1max = 0.0;
2217 double x3max = 0.0;
2218 double floatn = (double)v.length;
2219 double agiant = rgiant/floatn;
2220 for (int i=0;i<v.length;i++) {
2221 double xabs = Math.abs(v[i]);
2222 if (xabs<rdwarf || xabs>agiant) {
2223 if (xabs>rdwarf) {
2224 if (xabs>x1max) {
2225 double r=x1max/xabs;
2226 s1=1.0+s1*r*r;
2227 x1max=xabs;
2228 } else {
2229 double r=xabs/x1max;
2230 s1+=r*r;
2231 }
2232 } else {
2233 if (xabs>x3max) {
2234 double r=x3max/xabs;
2235 s3=1.0+s3*r*r;
2236 x3max=xabs;
2237 } else {
2238 if (xabs!=0.0) {
2239 double r=xabs/x3max;
2240 s3+=r*r;
2241 }
2242 }
2243 }
2244 } else {
2245 s2+=xabs*xabs;
2246 }
2247 }
2248 double norm;
2249 if (s1!=0.0) {
2250 norm = x1max*Math.sqrt(s1+(s2/x1max)/x1max);
2251 } else {
2252 if (s2==0.0) {
2253 norm = x3max*Math.sqrt(s3);
2254 } else {
2255 if (s2>=x3max) {
2256 norm = Math.sqrt(s2*(1.0+(x3max/s2)*(x3max*s3)));
2257 } else {
2258 norm = Math.sqrt(x3max*((s2/x3max)+(x3max*s3)));
2259 }
2260 }
2261 }
2262 return norm;
2263 }
2264
2265 }