View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math3.util;
19  
20  import java.lang.reflect.Array;
21  import java.util.ArrayList;
22  import java.util.Arrays;
23  import java.util.Collections;
24  import java.util.Comparator;
25  import java.util.List;
26  
27  import org.apache.commons.math3.Field;
28  import org.apache.commons.math3.random.RandomGenerator;
29  import org.apache.commons.math3.random.Well19937c;
30  import org.apache.commons.math3.distribution.UniformIntegerDistribution;
31  import org.apache.commons.math3.exception.DimensionMismatchException;
32  import org.apache.commons.math3.exception.MathArithmeticException;
33  import org.apache.commons.math3.exception.MathIllegalArgumentException;
34  import org.apache.commons.math3.exception.MathInternalError;
35  import org.apache.commons.math3.exception.NoDataException;
36  import org.apache.commons.math3.exception.NonMonotonicSequenceException;
37  import org.apache.commons.math3.exception.NotPositiveException;
38  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
39  import org.apache.commons.math3.exception.NullArgumentException;
40  import org.apache.commons.math3.exception.NumberIsTooLargeException;
41  import org.apache.commons.math3.exception.NotANumberException;
42  import org.apache.commons.math3.exception.util.LocalizedFormats;
43  
44  /**
45   * Arrays utilities.
46   *
47   * @since 3.0
48   */
49  public class MathArrays {
50      /** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
51      private static final int SPLIT_FACTOR = 0x8000001;
52  
53      /**
54       * Private constructor.
55       */
56      private MathArrays() {}
57  
58      /**
59       * Real-valued function that operate on an array or a part of it.
60       * @since 3.1
61       */
62      public interface Function {
63          /**
64           * Operates on an entire array.
65           *
66           * @param array Array to operate on.
67           * @return the result of the operation.
68           */
69          double evaluate(double[] array);
70          /**
71           * @param array Array to operate on.
72           * @param startIndex Index of the first element to take into account.
73           * @param numElements Number of elements to take into account.
74           * @return the result of the operation.
75           */
76          double evaluate(double[] array,
77                          int startIndex,
78                          int numElements);
79      }
80  
81      /**
82       * Create a copy of an array scaled by a value.
83       *
84       * @param arr Array to scale.
85       * @param val Scalar.
86       * @return scaled copy of array with each entry multiplied by val.
87       * @since 3.2
88       */
89      public static double[] scale(double val, final double[] arr) {
90          double[] newArr = new double[arr.length];
91          for (int i = 0; i < arr.length; i++) {
92              newArr[i] = arr[i] * val;
93          }
94          return newArr;
95      }
96  
97      /**
98       * <p>Multiply each element of an array by a value.</p>
99       *
100      * <p>The array is modified in place (no copy is created).</p>
101      *
102      * @param arr Array to scale
103      * @param val Scalar
104      * @since 3.2
105      */
106     public static void scaleInPlace(double val, final double[] arr) {
107         for (int i = 0; i < arr.length; i++) {
108             arr[i] *= val;
109         }
110     }
111 
112     /**
113      * Creates an array whose contents will be the element-by-element
114      * addition of the arguments.
115      *
116      * @param a First term of the addition.
117      * @param b Second term of the addition.
118      * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}.
119      * @throws DimensionMismatchException if the array lengths differ.
120      * @since 3.1
121      */
122     public static double[] ebeAdd(double[] a, double[] b)
123         throws DimensionMismatchException {
124         if (a.length != b.length) {
125             throw new DimensionMismatchException(a.length, b.length);
126         }
127 
128         final double[] result = a.clone();
129         for (int i = 0; i < a.length; i++) {
130             result[i] += b[i];
131         }
132         return result;
133     }
134     /**
135      * Creates an array whose contents will be the element-by-element
136      * subtraction of the second argument from the first.
137      *
138      * @param a First term.
139      * @param b Element to be subtracted.
140      * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}.
141      * @throws DimensionMismatchException if the array lengths differ.
142      * @since 3.1
143      */
144     public static double[] ebeSubtract(double[] a, double[] b)
145         throws DimensionMismatchException {
146         if (a.length != b.length) {
147             throw new DimensionMismatchException(a.length, b.length);
148         }
149 
150         final double[] result = a.clone();
151         for (int i = 0; i < a.length; i++) {
152             result[i] -= b[i];
153         }
154         return result;
155     }
156     /**
157      * Creates an array whose contents will be the element-by-element
158      * multiplication of the arguments.
159      *
160      * @param a First factor of the multiplication.
161      * @param b Second factor of the multiplication.
162      * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}.
163      * @throws DimensionMismatchException if the array lengths differ.
164      * @since 3.1
165      */
166     public static double[] ebeMultiply(double[] a, double[] b)
167         throws DimensionMismatchException {
168         if (a.length != b.length) {
169             throw new DimensionMismatchException(a.length, b.length);
170         }
171 
172         final double[] result = a.clone();
173         for (int i = 0; i < a.length; i++) {
174             result[i] *= b[i];
175         }
176         return result;
177     }
178     /**
179      * Creates an array whose contents will be the element-by-element
180      * division of the first argument by the second.
181      *
182      * @param a Numerator of the division.
183      * @param b Denominator of the division.
184      * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}.
185      * @throws DimensionMismatchException if the array lengths differ.
186      * @since 3.1
187      */
188     public static double[] ebeDivide(double[] a, double[] b)
189         throws DimensionMismatchException {
190         if (a.length != b.length) {
191             throw new DimensionMismatchException(a.length, b.length);
192         }
193 
194         final double[] result = a.clone();
195         for (int i = 0; i < a.length; i++) {
196             result[i] /= b[i];
197         }
198         return result;
199     }
200 
201     /**
202      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
203      *
204      * @param p1 the first point
205      * @param p2 the second point
206      * @return the L<sub>1</sub> distance between the two points
207      */
208     public static double distance1(double[] p1, double[] p2) {
209         double sum = 0;
210         for (int i = 0; i < p1.length; i++) {
211             sum += FastMath.abs(p1[i] - p2[i]);
212         }
213         return sum;
214     }
215 
216     /**
217      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
218      *
219      * @param p1 the first point
220      * @param p2 the second point
221      * @return the L<sub>1</sub> distance between the two points
222      */
223     public static int distance1(int[] p1, int[] p2) {
224       int sum = 0;
225       for (int i = 0; i < p1.length; i++) {
226           sum += FastMath.abs(p1[i] - p2[i]);
227       }
228       return sum;
229     }
230 
231     /**
232      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
233      *
234      * @param p1 the first point
235      * @param p2 the second point
236      * @return the L<sub>2</sub> distance between the two points
237      */
238     public static double distance(double[] p1, double[] p2) {
239         double sum = 0;
240         for (int i = 0; i < p1.length; i++) {
241             final double dp = p1[i] - p2[i];
242             sum += dp * dp;
243         }
244         return FastMath.sqrt(sum);
245     }
246 
247     /**
248      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
249      *
250      * @param p1 the first point
251      * @param p2 the second point
252      * @return the L<sub>2</sub> distance between the two points
253      */
254     public static double distance(int[] p1, int[] p2) {
255       double sum = 0;
256       for (int i = 0; i < p1.length; i++) {
257           final double dp = p1[i] - p2[i];
258           sum += dp * dp;
259       }
260       return FastMath.sqrt(sum);
261     }
262 
263     /**
264      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
265      *
266      * @param p1 the first point
267      * @param p2 the second point
268      * @return the L<sub>&infin;</sub> distance between the two points
269      */
270     public static double distanceInf(double[] p1, double[] p2) {
271         double max = 0;
272         for (int i = 0; i < p1.length; i++) {
273             max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
274         }
275         return max;
276     }
277 
278     /**
279      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
280      *
281      * @param p1 the first point
282      * @param p2 the second point
283      * @return the L<sub>&infin;</sub> distance between the two points
284      */
285     public static int distanceInf(int[] p1, int[] p2) {
286         int max = 0;
287         for (int i = 0; i < p1.length; i++) {
288             max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
289         }
290         return max;
291     }
292 
293     /**
294      * Specification of ordering direction.
295      */
296     public static enum OrderDirection {
297         /** Constant for increasing direction. */
298         INCREASING,
299         /** Constant for decreasing direction. */
300         DECREASING
301     }
302 
303     /**
304      * Check that an array is monotonically increasing or decreasing.
305      *
306      * @param <T> the type of the elements in the specified array
307      * @param val Values.
308      * @param dir Ordering direction.
309      * @param strict Whether the order should be strict.
310      * @return {@code true} if sorted, {@code false} otherwise.
311      */
312     public static  <T extends Comparable<? super T>> boolean isMonotonic(T[] val,
313                                       OrderDirection dir,
314                                       boolean strict) {
315         T previous = val[0];
316         final int max = val.length;
317         for (int i = 1; i < max; i++) {
318             final int comp;
319             switch (dir) {
320             case INCREASING:
321                 comp = previous.compareTo(val[i]);
322                 if (strict) {
323                     if (comp >= 0) {
324                         return false;
325                     }
326                 } else {
327                     if (comp > 0) {
328                         return false;
329                     }
330                 }
331                 break;
332             case DECREASING:
333                 comp = val[i].compareTo(previous);
334                 if (strict) {
335                     if (comp >= 0) {
336                         return false;
337                     }
338                 } else {
339                     if (comp > 0) {
340                        return false;
341                     }
342                 }
343                 break;
344             default:
345                 // Should never happen.
346                 throw new MathInternalError();
347             }
348 
349             previous = val[i];
350         }
351         return true;
352     }
353 
354     /**
355      * Check that an array is monotonically increasing or decreasing.
356      *
357      * @param val Values.
358      * @param dir Ordering direction.
359      * @param strict Whether the order should be strict.
360      * @return {@code true} if sorted, {@code false} otherwise.
361      */
362     public static boolean isMonotonic(double[] val, OrderDirection dir, boolean strict) {
363         return checkOrder(val, dir, strict, false);
364     }
365 
366     /**
367      * Check that the given array is sorted.
368      *
369      * @param val Values.
370      * @param dir Ordering direction.
371      * @param strict Whether the order should be strict.
372      * @param abort Whether to throw an exception if the check fails.
373      * @return {@code true} if the array is sorted.
374      * @throws NonMonotonicSequenceException if the array is not sorted
375      * and {@code abort} is {@code true}.
376      */
377     public static boolean checkOrder(double[] val, OrderDirection dir,
378                                      boolean strict, boolean abort)
379         throws NonMonotonicSequenceException {
380         double previous = val[0];
381         final int max = val.length;
382 
383         int index;
384         ITEM:
385         for (index = 1; index < max; index++) {
386             switch (dir) {
387             case INCREASING:
388                 if (strict) {
389                     if (val[index] <= previous) {
390                         break ITEM;
391                     }
392                 } else {
393                     if (val[index] < previous) {
394                         break ITEM;
395                     }
396                 }
397                 break;
398             case DECREASING:
399                 if (strict) {
400                     if (val[index] >= previous) {
401                         break ITEM;
402                     }
403                 } else {
404                     if (val[index] > previous) {
405                         break ITEM;
406                     }
407                 }
408                 break;
409             default:
410                 // Should never happen.
411                 throw new MathInternalError();
412             }
413 
414             previous = val[index];
415         }
416 
417         if (index == max) {
418             // Loop completed.
419             return true;
420         }
421 
422         // Loop early exit means wrong ordering.
423         if (abort) {
424             throw new NonMonotonicSequenceException(val[index], previous, index, dir, strict);
425         } else {
426             return false;
427         }
428     }
429 
430     /**
431      * Check that the given array is sorted.
432      *
433      * @param val Values.
434      * @param dir Ordering direction.
435      * @param strict Whether the order should be strict.
436      * @throws NonMonotonicSequenceException if the array is not sorted.
437      * @since 2.2
438      */
439     public static void checkOrder(double[] val, OrderDirection dir,
440                                   boolean strict) throws NonMonotonicSequenceException {
441         checkOrder(val, dir, strict, true);
442     }
443 
444     /**
445      * Check that the given array is sorted in strictly increasing order.
446      *
447      * @param val Values.
448      * @throws NonMonotonicSequenceException if the array is not sorted.
449      * @since 2.2
450      */
451     public static void checkOrder(double[] val) throws NonMonotonicSequenceException {
452         checkOrder(val, OrderDirection.INCREASING, true);
453     }
454 
455     /**
456      * Throws DimensionMismatchException if the input array is not rectangular.
457      *
458      * @param in array to be tested
459      * @throws NullArgumentException if input array is null
460      * @throws DimensionMismatchException if input array is not rectangular
461      * @since 3.1
462      */
463     public static void checkRectangular(final long[][] in)
464         throws NullArgumentException, DimensionMismatchException {
465         MathUtils.checkNotNull(in);
466         for (int i = 1; i < in.length; i++) {
467             if (in[i].length != in[0].length) {
468                 throw new DimensionMismatchException(
469                         LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
470                         in[i].length, in[0].length);
471             }
472         }
473     }
474 
475     /**
476      * Check that all entries of the input array are strictly positive.
477      *
478      * @param in Array to be tested
479      * @throws NotStrictlyPositiveException if any entries of the array are not
480      * strictly positive.
481      * @since 3.1
482      */
483     public static void checkPositive(final double[] in)
484         throws NotStrictlyPositiveException {
485         for (int i = 0; i < in.length; i++) {
486             if (in[i] <= 0) {
487                 throw new NotStrictlyPositiveException(in[i]);
488             }
489         }
490     }
491 
492     /**
493      * Check that no entry of the input array is {@code NaN}.
494      *
495      * @param in Array to be tested.
496      * @throws NotANumberException if an entry is {@code NaN}.
497      * @since 3.4
498      */
499     public static void checkNotNaN(final double[] in)
500         throws NotANumberException {
501         for(int i = 0; i < in.length; i++) {
502             if (Double.isNaN(in[i])) {
503                 throw new NotANumberException();
504             }
505         }
506     }
507 
508     /**
509      * Check that all entries of the input array are >= 0.
510      *
511      * @param in Array to be tested
512      * @throws NotPositiveException if any array entries are less than 0.
513      * @since 3.1
514      */
515     public static void checkNonNegative(final long[] in)
516         throws NotPositiveException {
517         for (int i = 0; i < in.length; i++) {
518             if (in[i] < 0) {
519                 throw new NotPositiveException(in[i]);
520             }
521         }
522     }
523 
524     /**
525      * Check all entries of the input array are >= 0.
526      *
527      * @param in Array to be tested
528      * @throws NotPositiveException if any array entries are less than 0.
529      * @since 3.1
530      */
531     public static void checkNonNegative(final long[][] in)
532         throws NotPositiveException {
533         for (int i = 0; i < in.length; i ++) {
534             for (int j = 0; j < in[i].length; j++) {
535                 if (in[i][j] < 0) {
536                     throw new NotPositiveException(in[i][j]);
537                 }
538             }
539         }
540     }
541 
542     /**
543      * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
544      * Translation of the minpack enorm subroutine.
545      *
546      * The redistribution policy for MINPACK is available
547      * <a href="http://www.netlib.org/minpack/disclaimer">here</a>, for
548      * convenience, it is reproduced below.</p>
549      *
550      * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
551      * <tr><td>
552      *    Minpack Copyright Notice (1999) University of Chicago.
553      *    All rights reserved
554      * </td></tr>
555      * <tr><td>
556      * Redistribution and use in source and binary forms, with or without
557      * modification, are permitted provided that the following conditions
558      * are met:
559      * <ol>
560      *  <li>Redistributions of source code must retain the above copyright
561      *      notice, this list of conditions and the following disclaimer.</li>
562      * <li>Redistributions in binary form must reproduce the above
563      *     copyright notice, this list of conditions and the following
564      *     disclaimer in the documentation and/or other materials provided
565      *     with the distribution.</li>
566      * <li>The end-user documentation included with the redistribution, if any,
567      *     must include the following acknowledgment:
568      *     {@code This product includes software developed by the University of
569      *           Chicago, as Operator of Argonne National Laboratory.}
570      *     Alternately, this acknowledgment may appear in the software itself,
571      *     if and wherever such third-party acknowledgments normally appear.</li>
572      * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
573      *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
574      *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
575      *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
576      *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
577      *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
578      *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
579      *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
580      *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
581      *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
582      *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
583      *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
584      *     BE CORRECTED.</strong></li>
585      * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
586      *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
587      *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
588      *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
589      *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
590      *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
591      *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
592      *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
593      *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
594      *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
595      * <ol></td></tr>
596      * </table>
597      *
598      * @param v Vector of doubles.
599      * @return the 2-norm of the vector.
600      * @since 2.2
601      */
602     public static double safeNorm(double[] v) {
603         double rdwarf = 3.834e-20;
604         double rgiant = 1.304e+19;
605         double s1 = 0;
606         double s2 = 0;
607         double s3 = 0;
608         double x1max = 0;
609         double x3max = 0;
610         double floatn = v.length;
611         double agiant = rgiant / floatn;
612         for (int i = 0; i < v.length; i++) {
613             double xabs = FastMath.abs(v[i]);
614             if (xabs < rdwarf || xabs > agiant) {
615                 if (xabs > rdwarf) {
616                     if (xabs > x1max) {
617                         double r = x1max / xabs;
618                         s1= 1 + s1 * r * r;
619                         x1max = xabs;
620                     } else {
621                         double r = xabs / x1max;
622                         s1 += r * r;
623                     }
624                 } else {
625                     if (xabs > x3max) {
626                         double r = x3max / xabs;
627                         s3= 1 + s3 * r * r;
628                         x3max = xabs;
629                     } else {
630                         if (xabs != 0) {
631                             double r = xabs / x3max;
632                             s3 += r * r;
633                         }
634                     }
635                 }
636             } else {
637                 s2 += xabs * xabs;
638             }
639         }
640         double norm;
641         if (s1 != 0) {
642             norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
643         } else {
644             if (s2 == 0) {
645                 norm = x3max * Math.sqrt(s3);
646             } else {
647                 if (s2 >= x3max) {
648                     norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
649                 } else {
650                     norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
651                 }
652             }
653         }
654         return norm;
655     }
656 
657     /**
658      * Sort an array in ascending order in place and perform the same reordering
659      * of entries on other arrays. For example, if
660      * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then
661      * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]},
662      * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}.
663      *
664      * @param x Array to be sorted and used as a pattern for permutation
665      * of the other arrays.
666      * @param yList Set of arrays whose permutations of entries will follow
667      * those performed on {@code x}.
668      * @throws DimensionMismatchException if any {@code y} is not the same
669      * size as {@code x}.
670      * @throws NullArgumentException if {@code x} or any {@code y} is null.
671      * @since 3.0
672      */
673     public static void sortInPlace(double[] x, double[] ... yList)
674         throws DimensionMismatchException, NullArgumentException {
675         sortInPlace(x, OrderDirection.INCREASING, yList);
676     }
677 
678     /**
679      * Sort an array in place and perform the same reordering of entries on
680      * other arrays.  This method works the same as the other
681      * {@link #sortInPlace(double[], double[][]) sortInPlace} method, but
682      * allows the order of the sort to be provided in the {@code dir}
683      * parameter.
684      *
685      * @param x Array to be sorted and used as a pattern for permutation
686      * of the other arrays.
687      * @param dir Order direction.
688      * @param yList Set of arrays whose permutations of entries will follow
689      * those performed on {@code x}.
690      * @throws DimensionMismatchException if any {@code y} is not the same
691      * size as {@code x}.
692      * @throws NullArgumentException if {@code x} or any {@code y} is null
693      * @since 3.0
694      */
695     public static void sortInPlace(double[] x,
696                                    final OrderDirection dir,
697                                    double[] ... yList)
698         throws NullArgumentException,
699                DimensionMismatchException {
700 
701         // Consistency checks.
702         if (x == null) {
703             throw new NullArgumentException();
704         }
705 
706         final int yListLen = yList.length;
707         final int len = x.length;
708 
709         for (int j = 0; j < yListLen; j++) {
710             final double[] y = yList[j];
711             if (y == null) {
712                 throw new NullArgumentException();
713             }
714             if (y.length != len) {
715                 throw new DimensionMismatchException(y.length, len);
716             }
717         }
718 
719         // Associate each abscissa "x[i]" with its index "i".
720         final List<Pair<Double, Integer>> list
721             = new ArrayList<Pair<Double, Integer>>(len);
722         for (int i = 0; i < len; i++) {
723             list.add(new Pair<Double, Integer>(x[i], i));
724         }
725 
726         // Create comparators for increasing and decreasing orders.
727         final Comparator<Pair<Double, Integer>> comp
728             = dir == MathArrays.OrderDirection.INCREASING ?
729             new Comparator<Pair<Double, Integer>>() {
730             public int compare(Pair<Double, Integer> o1,
731                                Pair<Double, Integer> o2) {
732                 return o1.getKey().compareTo(o2.getKey());
733             }
734         } : new Comparator<Pair<Double,Integer>>() {
735             public int compare(Pair<Double, Integer> o1,
736                                Pair<Double, Integer> o2) {
737                 return o2.getKey().compareTo(o1.getKey());
738             }
739         };
740 
741         // Sort.
742         Collections.sort(list, comp);
743 
744         // Modify the original array so that its elements are in
745         // the prescribed order.
746         // Retrieve indices of original locations.
747         final int[] indices = new int[len];
748         for (int i = 0; i < len; i++) {
749             final Pair<Double, Integer> e = list.get(i);
750             x[i] = e.getKey();
751             indices[i] = e.getValue();
752         }
753 
754         // In each of the associated arrays, move the
755         // elements to their new location.
756         for (int j = 0; j < yListLen; j++) {
757             // Input array will be modified in place.
758             final double[] yInPlace = yList[j];
759             final double[] yOrig = yInPlace.clone();
760 
761             for (int i = 0; i < len; i++) {
762                 yInPlace[i] = yOrig[indices[i]];
763             }
764         }
765     }
766 
767     /**
768      * Creates a copy of the {@code source} array.
769      *
770      * @param source Array to be copied.
771      * @return the copied array.
772      */
773      public static int[] copyOf(int[] source) {
774          return copyOf(source, source.length);
775      }
776 
777     /**
778      * Creates a copy of the {@code source} array.
779      *
780      * @param source Array to be copied.
781      * @return the copied array.
782      */
783      public static double[] copyOf(double[] source) {
784          return copyOf(source, source.length);
785      }
786 
787     /**
788      * Creates a copy of the {@code source} array.
789      *
790      * @param source Array to be copied.
791      * @param len Number of entries to copy. If smaller then the source
792      * length, the copy will be truncated, if larger it will padded with
793      * zeroes.
794      * @return the copied array.
795      */
796     public static int[] copyOf(int[] source, int len) {
797          final int[] output = new int[len];
798          System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
799          return output;
800      }
801 
802     /**
803      * Creates a copy of the {@code source} array.
804      *
805      * @param source Array to be copied.
806      * @param len Number of entries to copy. If smaller then the source
807      * length, the copy will be truncated, if larger it will padded with
808      * zeroes.
809      * @return the copied array.
810      */
811     public static double[] copyOf(double[] source, int len) {
812          final double[] output = new double[len];
813          System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
814          return output;
815      }
816 
817     /**
818      * Creates a copy of the {@code source} array.
819      *
820      * @param source Array to be copied.
821      * @param from Initial index of the range to be copied, inclusive.
822      * @param to Final index of the range to be copied, exclusive. (This index may lie outside the array.)
823      * @return the copied array.
824      */
825     public static double[] copyOfRange(double[] source, int from, int to) {
826         final int len = to - from;
827         final double[] output = new double[len];
828         System.arraycopy(source, from, output, 0, FastMath.min(len, source.length - from));
829         return output;
830      }
831 
832     /**
833      * Compute a linear combination accurately.
834      * This method computes the sum of the products
835      * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
836      * It does so by using specific multiplication and addition algorithms to
837      * preserve accuracy and reduce cancellation effects.
838      * <br/>
839      * It is based on the 2005 paper
840      * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
841      * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
842      * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
843      *
844      * @param a Factors.
845      * @param b Factors.
846      * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
847      * @throws DimensionMismatchException if arrays dimensions don't match
848      */
849     public static double linearCombination(final double[] a, final double[] b)
850         throws DimensionMismatchException {
851         final int len = a.length;
852         if (len != b.length) {
853             throw new DimensionMismatchException(len, b.length);
854         }
855 
856         if (len == 1) {
857             // Revert to scalar multiplication.
858             return a[0] * b[0];
859         }
860 
861         final double[] prodHigh = new double[len];
862         double prodLowSum = 0;
863 
864         for (int i = 0; i < len; i++) {
865             final double ai = a[i];
866             final double ca = SPLIT_FACTOR * ai;
867             final double aHigh = ca - (ca - ai);
868             final double aLow = ai - aHigh;
869 
870             final double bi = b[i];
871             final double cb = SPLIT_FACTOR * bi;
872             final double bHigh = cb - (cb - bi);
873             final double bLow = bi - bHigh;
874             prodHigh[i] = ai * bi;
875             final double prodLow = aLow * bLow - (((prodHigh[i] -
876                                                     aHigh * bHigh) -
877                                                    aLow * bHigh) -
878                                                   aHigh * bLow);
879             prodLowSum += prodLow;
880         }
881 
882 
883         final double prodHighCur = prodHigh[0];
884         double prodHighNext = prodHigh[1];
885         double sHighPrev = prodHighCur + prodHighNext;
886         double sPrime = sHighPrev - prodHighNext;
887         double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime);
888 
889         final int lenMinusOne = len - 1;
890         for (int i = 1; i < lenMinusOne; i++) {
891             prodHighNext = prodHigh[i + 1];
892             final double sHighCur = sHighPrev + prodHighNext;
893             sPrime = sHighCur - prodHighNext;
894             sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime);
895             sHighPrev = sHighCur;
896         }
897 
898         double result = sHighPrev + (prodLowSum + sLowSum);
899 
900         if (Double.isNaN(result)) {
901             // either we have split infinite numbers or some coefficients were NaNs,
902             // just rely on the naive implementation and let IEEE754 handle this
903             result = 0;
904             for (int i = 0; i < len; ++i) {
905                 result += a[i] * b[i];
906             }
907         }
908 
909         return result;
910     }
911 
912     /**
913      * Compute a linear combination accurately.
914      * <p>
915      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
916      * a<sub>2</sub>&times;b<sub>2</sub> to high accuracy. It does
917      * so by using specific multiplication and addition algorithms to
918      * preserve accuracy and reduce cancellation effects. It is based
919      * on the 2005 paper <a
920      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
921      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
922      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
923      * </p>
924      * @param a1 first factor of the first term
925      * @param b1 second factor of the first term
926      * @param a2 first factor of the second term
927      * @param b2 second factor of the second term
928      * @return a<sub>1</sub>&times;b<sub>1</sub> +
929      * a<sub>2</sub>&times;b<sub>2</sub>
930      * @see #linearCombination(double, double, double, double, double, double)
931      * @see #linearCombination(double, double, double, double, double, double, double, double)
932      */
933     public static double linearCombination(final double a1, final double b1,
934                                            final double a2, final double b2) {
935 
936         // the code below is split in many additions/subtractions that may
937         // appear redundant. However, they should NOT be simplified, as they
938         // use IEEE754 floating point arithmetic rounding properties.
939         // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
940         // The variable naming conventions are that xyzHigh contains the most significant
941         // bits of xyz and xyzLow contains its least significant bits. So theoretically
942         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
943         // be represented in only one double precision number so we preserve two numbers
944         // to hold it as long as we can, combining the high and low order bits together
945         // only at the end, after cancellation may have occurred on high order bits
946 
947         // split a1 and b1 as two 26 bits numbers
948         final double ca1        = SPLIT_FACTOR * a1;
949         final double a1High     = ca1 - (ca1 - a1);
950         final double a1Low      = a1 - a1High;
951         final double cb1        = SPLIT_FACTOR * b1;
952         final double b1High     = cb1 - (cb1 - b1);
953         final double b1Low      = b1 - b1High;
954 
955         // accurate multiplication a1 * b1
956         final double prod1High  = a1 * b1;
957         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
958 
959         // split a2 and b2 as two 26 bits numbers
960         final double ca2        = SPLIT_FACTOR * a2;
961         final double a2High     = ca2 - (ca2 - a2);
962         final double a2Low      = a2 - a2High;
963         final double cb2        = SPLIT_FACTOR * b2;
964         final double b2High     = cb2 - (cb2 - b2);
965         final double b2Low      = b2 - b2High;
966 
967         // accurate multiplication a2 * b2
968         final double prod2High  = a2 * b2;
969         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
970 
971         // accurate addition a1 * b1 + a2 * b2
972         final double s12High    = prod1High + prod2High;
973         final double s12Prime   = s12High - prod2High;
974         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
975 
976         // final rounding, s12 may have suffered many cancellations, we try
977         // to recover some bits from the extra words we have saved up to now
978         double result = s12High + (prod1Low + prod2Low + s12Low);
979 
980         if (Double.isNaN(result)) {
981             // either we have split infinite numbers or some coefficients were NaNs,
982             // just rely on the naive implementation and let IEEE754 handle this
983             result = a1 * b1 + a2 * b2;
984         }
985 
986         return result;
987     }
988 
989     /**
990      * Compute a linear combination accurately.
991      * <p>
992      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
993      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
994      * to high accuracy. It does so by using specific multiplication and
995      * addition algorithms to preserve accuracy and reduce cancellation effects.
996      * It is based on the 2005 paper <a
997      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
998      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
999      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
1000      * </p>
1001      * @param a1 first factor of the first term
1002      * @param b1 second factor of the first term
1003      * @param a2 first factor of the second term
1004      * @param b2 second factor of the second term
1005      * @param a3 first factor of the third term
1006      * @param b3 second factor of the third term
1007      * @return a<sub>1</sub>&times;b<sub>1</sub> +
1008      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
1009      * @see #linearCombination(double, double, double, double)
1010      * @see #linearCombination(double, double, double, double, double, double, double, double)
1011      */
1012     public static double linearCombination(final double a1, final double b1,
1013                                            final double a2, final double b2,
1014                                            final double a3, final double b3) {
1015 
1016         // the code below is split in many additions/subtractions that may
1017         // appear redundant. However, they should NOT be simplified, as they
1018         // do use IEEE754 floating point arithmetic rounding properties.
1019         // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
1020         // The variables naming conventions are that xyzHigh contains the most significant
1021         // bits of xyz and xyzLow contains its least significant bits. So theoretically
1022         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
1023         // be represented in only one double precision number so we preserve two numbers
1024         // to hold it as long as we can, combining the high and low order bits together
1025         // only at the end, after cancellation may have occurred on high order bits
1026 
1027         // split a1 and b1 as two 26 bits numbers
1028         final double ca1        = SPLIT_FACTOR * a1;
1029         final double a1High     = ca1 - (ca1 - a1);
1030         final double a1Low      = a1 - a1High;
1031         final double cb1        = SPLIT_FACTOR * b1;
1032         final double b1High     = cb1 - (cb1 - b1);
1033         final double b1Low      = b1 - b1High;
1034 
1035         // accurate multiplication a1 * b1
1036         final double prod1High  = a1 * b1;
1037         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
1038 
1039         // split a2 and b2 as two 26 bits numbers
1040         final double ca2        = SPLIT_FACTOR * a2;
1041         final double a2High     = ca2 - (ca2 - a2);
1042         final double a2Low      = a2 - a2High;
1043         final double cb2        = SPLIT_FACTOR * b2;
1044         final double b2High     = cb2 - (cb2 - b2);
1045         final double b2Low      = b2 - b2High;
1046 
1047         // accurate multiplication a2 * b2
1048         final double prod2High  = a2 * b2;
1049         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
1050 
1051         // split a3 and b3 as two 26 bits numbers
1052         final double ca3        = SPLIT_FACTOR * a3;
1053         final double a3High     = ca3 - (ca3 - a3);
1054         final double a3Low      = a3 - a3High;
1055         final double cb3        = SPLIT_FACTOR * b3;
1056         final double b3High     = cb3 - (cb3 - b3);
1057         final double b3Low      = b3 - b3High;
1058 
1059         // accurate multiplication a3 * b3
1060         final double prod3High  = a3 * b3;
1061         final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
1062 
1063         // accurate addition a1 * b1 + a2 * b2
1064         final double s12High    = prod1High + prod2High;
1065         final double s12Prime   = s12High - prod2High;
1066         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
1067 
1068         // accurate addition a1 * b1 + a2 * b2 + a3 * b3
1069         final double s123High   = s12High + prod3High;
1070         final double s123Prime  = s123High - prod3High;
1071         final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
1072 
1073         // final rounding, s123 may have suffered many cancellations, we try
1074         // to recover some bits from the extra words we have saved up to now
1075         double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low);
1076 
1077         if (Double.isNaN(result)) {
1078             // either we have split infinite numbers or some coefficients were NaNs,
1079             // just rely on the naive implementation and let IEEE754 handle this
1080             result = a1 * b1 + a2 * b2 + a3 * b3;
1081         }
1082 
1083         return result;
1084     }
1085 
1086     /**
1087      * Compute a linear combination accurately.
1088      * <p>
1089      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
1090      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
1091      * a<sub>4</sub>&times;b<sub>4</sub>
1092      * to high accuracy. It does so by using specific multiplication and
1093      * addition algorithms to preserve accuracy and reduce cancellation effects.
1094      * It is based on the 2005 paper <a
1095      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
1096      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
1097      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
1098      * </p>
1099      * @param a1 first factor of the first term
1100      * @param b1 second factor of the first term
1101      * @param a2 first factor of the second term
1102      * @param b2 second factor of the second term
1103      * @param a3 first factor of the third term
1104      * @param b3 second factor of the third term
1105      * @param a4 first factor of the third term
1106      * @param b4 second factor of the third term
1107      * @return a<sub>1</sub>&times;b<sub>1</sub> +
1108      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
1109      * a<sub>4</sub>&times;b<sub>4</sub>
1110      * @see #linearCombination(double, double, double, double)
1111      * @see #linearCombination(double, double, double, double, double, double)
1112      */
1113     public static double linearCombination(final double a1, final double b1,
1114                                            final double a2, final double b2,
1115                                            final double a3, final double b3,
1116                                            final double a4, final double b4) {
1117 
1118         // the code below is split in many additions/subtractions that may
1119         // appear redundant. However, they should NOT be simplified, as they
1120         // do use IEEE754 floating point arithmetic rounding properties.
1121         // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
1122         // The variables naming conventions are that xyzHigh contains the most significant
1123         // bits of xyz and xyzLow contains its least significant bits. So theoretically
1124         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
1125         // be represented in only one double precision number so we preserve two numbers
1126         // to hold it as long as we can, combining the high and low order bits together
1127         // only at the end, after cancellation may have occurred on high order bits
1128 
1129         // split a1 and b1 as two 26 bits numbers
1130         final double ca1        = SPLIT_FACTOR * a1;
1131         final double a1High     = ca1 - (ca1 - a1);
1132         final double a1Low      = a1 - a1High;
1133         final double cb1        = SPLIT_FACTOR * b1;
1134         final double b1High     = cb1 - (cb1 - b1);
1135         final double b1Low      = b1 - b1High;
1136 
1137         // accurate multiplication a1 * b1
1138         final double prod1High  = a1 * b1;
1139         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
1140 
1141         // split a2 and b2 as two 26 bits numbers
1142         final double ca2        = SPLIT_FACTOR * a2;
1143         final double a2High     = ca2 - (ca2 - a2);
1144         final double a2Low      = a2 - a2High;
1145         final double cb2        = SPLIT_FACTOR * b2;
1146         final double b2High     = cb2 - (cb2 - b2);
1147         final double b2Low      = b2 - b2High;
1148 
1149         // accurate multiplication a2 * b2
1150         final double prod2High  = a2 * b2;
1151         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
1152 
1153         // split a3 and b3 as two 26 bits numbers
1154         final double ca3        = SPLIT_FACTOR * a3;
1155         final double a3High     = ca3 - (ca3 - a3);
1156         final double a3Low      = a3 - a3High;
1157         final double cb3        = SPLIT_FACTOR * b3;
1158         final double b3High     = cb3 - (cb3 - b3);
1159         final double b3Low      = b3 - b3High;
1160 
1161         // accurate multiplication a3 * b3
1162         final double prod3High  = a3 * b3;
1163         final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
1164 
1165         // split a4 and b4 as two 26 bits numbers
1166         final double ca4        = SPLIT_FACTOR * a4;
1167         final double a4High     = ca4 - (ca4 - a4);
1168         final double a4Low      = a4 - a4High;
1169         final double cb4        = SPLIT_FACTOR * b4;
1170         final double b4High     = cb4 - (cb4 - b4);
1171         final double b4Low      = b4 - b4High;
1172 
1173         // accurate multiplication a4 * b4
1174         final double prod4High  = a4 * b4;
1175         final double prod4Low   = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low);
1176 
1177         // accurate addition a1 * b1 + a2 * b2
1178         final double s12High    = prod1High + prod2High;
1179         final double s12Prime   = s12High - prod2High;
1180         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
1181 
1182         // accurate addition a1 * b1 + a2 * b2 + a3 * b3
1183         final double s123High   = s12High + prod3High;
1184         final double s123Prime  = s123High - prod3High;
1185         final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
1186 
1187         // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4
1188         final double s1234High  = s123High + prod4High;
1189         final double s1234Prime = s1234High - prod4High;
1190         final double s1234Low   = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime);
1191 
1192         // final rounding, s1234 may have suffered many cancellations, we try
1193         // to recover some bits from the extra words we have saved up to now
1194         double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low);
1195 
1196         if (Double.isNaN(result)) {
1197             // either we have split infinite numbers or some coefficients were NaNs,
1198             // just rely on the naive implementation and let IEEE754 handle this
1199             result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
1200         }
1201 
1202         return result;
1203     }
1204 
1205     /**
1206      * Returns true iff both arguments are null or have same dimensions and all
1207      * their elements are equal as defined by
1208      * {@link Precision#equals(float,float)}.
1209      *
1210      * @param x first array
1211      * @param y second array
1212      * @return true if the values are both null or have same dimension
1213      * and equal elements.
1214      */
1215     public static boolean equals(float[] x, float[] y) {
1216         if ((x == null) || (y == null)) {
1217             return !((x == null) ^ (y == null));
1218         }
1219         if (x.length != y.length) {
1220             return false;
1221         }
1222         for (int i = 0; i < x.length; ++i) {
1223             if (!Precision.equals(x[i], y[i])) {
1224                 return false;
1225             }
1226         }
1227         return true;
1228     }
1229 
1230     /**
1231      * Returns true iff both arguments are null or have same dimensions and all
1232      * their elements are equal as defined by
1233      * {@link Precision#equalsIncludingNaN(double,double) this method}.
1234      *
1235      * @param x first array
1236      * @param y second array
1237      * @return true if the values are both null or have same dimension and
1238      * equal elements
1239      * @since 2.2
1240      */
1241     public static boolean equalsIncludingNaN(float[] x, float[] y) {
1242         if ((x == null) || (y == null)) {
1243             return !((x == null) ^ (y == null));
1244         }
1245         if (x.length != y.length) {
1246             return false;
1247         }
1248         for (int i = 0; i < x.length; ++i) {
1249             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
1250                 return false;
1251             }
1252         }
1253         return true;
1254     }
1255 
1256     /**
1257      * Returns {@code true} iff both arguments are {@code null} or have same
1258      * dimensions and all their elements are equal as defined by
1259      * {@link Precision#equals(double,double)}.
1260      *
1261      * @param x First array.
1262      * @param y Second array.
1263      * @return {@code true} if the values are both {@code null} or have same
1264      * dimension and equal elements.
1265      */
1266     public static boolean equals(double[] x, double[] y) {
1267         if ((x == null) || (y == null)) {
1268             return !((x == null) ^ (y == null));
1269         }
1270         if (x.length != y.length) {
1271             return false;
1272         }
1273         for (int i = 0; i < x.length; ++i) {
1274             if (!Precision.equals(x[i], y[i])) {
1275                 return false;
1276             }
1277         }
1278         return true;
1279     }
1280 
1281     /**
1282      * Returns {@code true} iff both arguments are {@code null} or have same
1283      * dimensions and all their elements are equal as defined by
1284      * {@link Precision#equalsIncludingNaN(double,double) this method}.
1285      *
1286      * @param x First array.
1287      * @param y Second array.
1288      * @return {@code true} if the values are both {@code null} or have same
1289      * dimension and equal elements.
1290      * @since 2.2
1291      */
1292     public static boolean equalsIncludingNaN(double[] x, double[] y) {
1293         if ((x == null) || (y == null)) {
1294             return !((x == null) ^ (y == null));
1295         }
1296         if (x.length != y.length) {
1297             return false;
1298         }
1299         for (int i = 0; i < x.length; ++i) {
1300             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
1301                 return false;
1302             }
1303         }
1304         return true;
1305     }
1306 
1307     /**
1308      * Normalizes an array to make it sum to a specified value.
1309      * Returns the result of the transformation
1310      * <pre>
1311      *    x |-> x * normalizedSum / sum
1312      * </pre>
1313      * applied to each non-NaN element x of the input array, where sum is the
1314      * sum of the non-NaN entries in the input array.
1315      * <p>
1316      * Throws IllegalArgumentException if {@code normalizedSum} is infinite
1317      * or NaN and ArithmeticException if the input array contains any infinite elements
1318      * or sums to 0.
1319      * <p>
1320      * Ignores (i.e., copies unchanged to the output array) NaNs in the input array.
1321      *
1322      * @param values Input array to be normalized
1323      * @param normalizedSum Target sum for the normalized array
1324      * @return the normalized array.
1325      * @throws MathArithmeticException if the input array contains infinite
1326      * elements or sums to zero.
1327      * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}.
1328      * @since 2.1
1329      */
1330     public static double[] normalizeArray(double[] values, double normalizedSum)
1331         throws MathIllegalArgumentException, MathArithmeticException {
1332         if (Double.isInfinite(normalizedSum)) {
1333             throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_INFINITE);
1334         }
1335         if (Double.isNaN(normalizedSum)) {
1336             throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_NAN);
1337         }
1338         double sum = 0d;
1339         final int len = values.length;
1340         double[] out = new double[len];
1341         for (int i = 0; i < len; i++) {
1342             if (Double.isInfinite(values[i])) {
1343                 throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
1344             }
1345             if (!Double.isNaN(values[i])) {
1346                 sum += values[i];
1347             }
1348         }
1349         if (sum == 0) {
1350             throw new MathArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
1351         }
1352         for (int i = 0; i < len; i++) {
1353             if (Double.isNaN(values[i])) {
1354                 out[i] = Double.NaN;
1355             } else {
1356                 out[i] = values[i] * normalizedSum / sum;
1357             }
1358         }
1359         return out;
1360     }
1361 
1362     /** Build an array of elements.
1363      * <p>
1364      * Arrays are filled with field.getZero()
1365      *
1366      * @param <T> the type of the field elements
1367      * @param field field to which array elements belong
1368      * @param length of the array
1369      * @return a new array
1370      * @since 3.2
1371      */
1372     public static <T> T[] buildArray(final Field<T> field, final int length) {
1373         @SuppressWarnings("unchecked") // OK because field must be correct class
1374         T[] array = (T[]) Array.newInstance(field.getRuntimeClass(), length);
1375         Arrays.fill(array, field.getZero());
1376         return array;
1377     }
1378 
1379     /** Build a double dimension  array of elements.
1380      * <p>
1381      * Arrays are filled with field.getZero()
1382      *
1383      * @param <T> the type of the field elements
1384      * @param field field to which array elements belong
1385      * @param rows number of rows in the array
1386      * @param columns number of columns (may be negative to build partial
1387      * arrays in the same way <code>new Field[rows][]</code> works)
1388      * @return a new array
1389      * @since 3.2
1390      */
1391     @SuppressWarnings("unchecked")
1392     public static <T> T[][] buildArray(final Field<T> field, final int rows, final int columns) {
1393         final T[][] array;
1394         if (columns < 0) {
1395             T[] dummyRow = buildArray(field, 0);
1396             array = (T[][]) Array.newInstance(dummyRow.getClass(), rows);
1397         } else {
1398             array = (T[][]) Array.newInstance(field.getRuntimeClass(),
1399                                               new int[] {
1400                                                   rows, columns
1401                                               });
1402             for (int i = 0; i < rows; ++i) {
1403                 Arrays.fill(array[i], field.getZero());
1404             }
1405         }
1406         return array;
1407     }
1408 
1409     /**
1410      * Calculates the <a href="http://en.wikipedia.org/wiki/Convolution">
1411      * convolution</a> between two sequences.
1412      * <p>
1413      * The solution is obtained via straightforward computation of the
1414      * convolution sum (and not via FFT). Whenever the computation needs
1415      * an element that would be located at an index outside the input arrays,
1416      * the value is assumed to be zero.
1417      *
1418      * @param x First sequence.
1419      * Typically, this sequence will represent an input signal to a system.
1420      * @param h Second sequence.
1421      * Typically, this sequence will represent the impulse response of the system.
1422      * @return the convolution of {@code x} and {@code h}.
1423      * This array's length will be {@code x.length + h.length - 1}.
1424      * @throws NullArgumentException if either {@code x} or {@code h} is {@code null}.
1425      * @throws NoDataException if either {@code x} or {@code h} is empty.
1426      *
1427      * @since 3.3
1428      */
1429     public static double[] convolve(double[] x, double[] h)
1430         throws NullArgumentException,
1431                NoDataException {
1432         MathUtils.checkNotNull(x);
1433         MathUtils.checkNotNull(h);
1434 
1435         final int xLen = x.length;
1436         final int hLen = h.length;
1437 
1438         if (xLen == 0 || hLen == 0) {
1439             throw new NoDataException();
1440         }
1441 
1442         // initialize the output array
1443         final int totalLength = xLen + hLen - 1;
1444         final double[] y = new double[totalLength];
1445 
1446         // straightforward implementation of the convolution sum
1447         for (int n = 0; n < totalLength; n++) {
1448             double yn = 0;
1449             int k = FastMath.max(0, n + 1 - xLen);
1450             int j = n - k;
1451             while (k < hLen && j >= 0) {
1452                 yn += x[j--] * h[k++];
1453             }
1454             y[n] = yn;
1455         }
1456 
1457         return y;
1458     }
1459 
1460     /**
1461      * Specification for indicating that some operation applies
1462      * before or after a given index.
1463      */
1464     public static enum Position {
1465         /** Designates the beginning of the array (near index 0). */
1466         HEAD,
1467         /** Designates the end of the array. */
1468         TAIL
1469     }
1470 
1471     /**
1472      * Shuffle the entries of the given array.
1473      * The {@code start} and {@code pos} parameters select which portion
1474      * of the array is randomized and which is left untouched.
1475      *
1476      * @see #shuffle(int[],int,Position,RandomGenerator)
1477      *
1478      * @param list Array whose entries will be shuffled (in-place).
1479      * @param start Index at which shuffling begins.
1480      * @param pos Shuffling is performed for index positions between
1481      * {@code start} and either the end (if {@link Position#TAIL})
1482      * or the beginning (if {@link Position#HEAD}) of the array.
1483      */
1484     public static void shuffle(int[] list,
1485                                int start,
1486                                Position pos) {
1487         shuffle(list, start, pos, new Well19937c());
1488     }
1489 
1490     /**
1491      * Shuffle the entries of the given array, using the
1492      * <a href="http://en.wikipedia.org/wiki/Fisher–Yates_shuffle#The_modern_algorithm">
1493      * Fisher–Yates</a> algorithm.
1494      * The {@code start} and {@code pos} parameters select which portion
1495      * of the array is randomized and which is left untouched.
1496      *
1497      * @param list Array whose entries will be shuffled (in-place).
1498      * @param start Index at which shuffling begins.
1499      * @param pos Shuffling is performed for index positions between
1500      * {@code start} and either the end (if {@link Position#TAIL})
1501      * or the beginning (if {@link Position#HEAD}) of the array.
1502      * @param rng Random number generator.
1503      */
1504     public static void shuffle(int[] list,
1505                                int start,
1506                                Position pos,
1507                                RandomGenerator rng) {
1508         switch (pos) {
1509         case TAIL: {
1510             for (int i = list.length - 1; i >= start; i--) {
1511                 final int target;
1512                 if (i == start) {
1513                     target = start;
1514                 } else {
1515                     // NumberIsTooLargeException cannot occur.
1516                     target = new UniformIntegerDistribution(rng, start, i).sample();
1517                 }
1518                 final int temp = list[target];
1519                 list[target] = list[i];
1520                 list[i] = temp;
1521             }
1522         }
1523             break;
1524         case HEAD: {
1525             for (int i = 0; i <= start; i++) {
1526                 final int target;
1527                 if (i == start) {
1528                     target = start;
1529                 } else {
1530                     // NumberIsTooLargeException cannot occur.
1531                     target = new UniformIntegerDistribution(rng, i, start).sample();
1532                 }
1533                 final int temp = list[target];
1534                 list[target] = list[i];
1535                 list[i] = temp;
1536             }
1537         }
1538             break;
1539         default:
1540             throw new MathInternalError(); // Should never happen.
1541         }
1542     }
1543 
1544     /**
1545      * Shuffle the entries of the given array.
1546      *
1547      * @see #shuffle(int[],int,Position,RandomGenerator)
1548      *
1549      * @param list Array whose entries will be shuffled (in-place).
1550      * @param rng Random number generator.
1551      */
1552     public static void shuffle(int[] list,
1553                                RandomGenerator rng) {
1554         shuffle(list, 0, Position.TAIL, rng);
1555     }
1556 
1557     /**
1558      * Shuffle the entries of the given array.
1559      *
1560      * @see #shuffle(int[],int,Position,RandomGenerator)
1561      *
1562      * @param list Array whose entries will be shuffled (in-place).
1563      */
1564     public static void shuffle(int[] list) {
1565         shuffle(list, new Well19937c());
1566     }
1567 
1568     /**
1569      * Returns an array representing the natural number {@code n}.
1570      *
1571      * @param n Natural number.
1572      * @return an array whose entries are the numbers 0, 1, ..., {@code n}-1.
1573      * If {@code n == 0}, the returned array is empty.
1574      */
1575     public static int[] natural(int n) {
1576         final int[] a = new int[n];
1577         for (int i = 0; i < n; i++) {
1578             a[i] = i;
1579         }
1580         return a;
1581     }
1582     /**
1583      * This method is used
1584      * to verify that the input parameters designate a subarray of positive length.
1585      * <p>
1586      * <ul>
1587      * <li>returns <code>true</code> iff the parameters designate a subarray of
1588      * positive length</li>
1589      * <li>throws <code>MathIllegalArgumentException</code> if the array is null or
1590      * or the indices are invalid</li>
1591      * <li>returns <code>false</li> if the array is non-null, but
1592      * <code>length</code> is 0.
1593      * </ul></p>
1594      *
1595      * @param values the input array
1596      * @param begin index of the first array element to include
1597      * @param length the number of elements to include
1598      * @return true if the parameters are valid and designate a subarray of positive length
1599      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
1600      * @since 3.3
1601      */
1602     public static boolean verifyValues(final double[] values, final int begin, final int length)
1603             throws MathIllegalArgumentException {
1604         return verifyValues(values, begin, length, false);
1605     }
1606 
1607     /**
1608      * This method is used
1609      * to verify that the input parameters designate a subarray of positive length.
1610      * <p>
1611      * <ul>
1612      * <li>returns <code>true</code> iff the parameters designate a subarray of
1613      * non-negative length</li>
1614      * <li>throws <code>IllegalArgumentException</code> if the array is null or
1615      * or the indices are invalid</li>
1616      * <li>returns <code>false</li> if the array is non-null, but
1617      * <code>length</code> is 0 unless <code>allowEmpty</code> is <code>true</code>
1618      * </ul></p>
1619      *
1620      * @param values the input array
1621      * @param begin index of the first array element to include
1622      * @param length the number of elements to include
1623      * @param allowEmpty if <code>true</code> then zero length arrays are allowed
1624      * @return true if the parameters are valid
1625      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
1626      * @since 3.3
1627      */
1628     public static boolean verifyValues(final double[] values, final int begin,
1629             final int length, final boolean allowEmpty) throws MathIllegalArgumentException {
1630 
1631         if (values == null) {
1632             throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
1633         }
1634 
1635         if (begin < 0) {
1636             throw new NotPositiveException(LocalizedFormats.START_POSITION, Integer.valueOf(begin));
1637         }
1638 
1639         if (length < 0) {
1640             throw new NotPositiveException(LocalizedFormats.LENGTH, Integer.valueOf(length));
1641         }
1642 
1643         if (begin + length > values.length) {
1644             throw new NumberIsTooLargeException(LocalizedFormats.SUBARRAY_ENDS_AFTER_ARRAY_END,
1645                     Integer.valueOf(begin + length), Integer.valueOf(values.length), true);
1646         }
1647 
1648         if (length == 0 && !allowEmpty) {
1649             return false;
1650         }
1651 
1652         return true;
1653 
1654     }
1655 
1656     /**
1657      * This method is used
1658      * to verify that the begin and length parameters designate a subarray of positive length
1659      * and the weights are all non-negative, non-NaN, finite, and not all zero.
1660      * <p>
1661      * <ul>
1662      * <li>returns <code>true</code> iff the parameters designate a subarray of
1663      * positive length and the weights array contains legitimate values.</li>
1664      * <li>throws <code>IllegalArgumentException</code> if any of the following are true:
1665      * <ul><li>the values array is null</li>
1666      *     <li>the weights array is null</li>
1667      *     <li>the weights array does not have the same length as the values array</li>
1668      *     <li>the weights array contains one or more infinite values</li>
1669      *     <li>the weights array contains one or more NaN values</li>
1670      *     <li>the weights array contains negative values</li>
1671      *     <li>the start and length arguments do not determine a valid array</li></ul>
1672      * </li>
1673      * <li>returns <code>false</li> if the array is non-null, but
1674      * <code>length</code> is 0.
1675      * </ul></p>
1676      *
1677      * @param values the input array
1678      * @param weights the weights array
1679      * @param begin index of the first array element to include
1680      * @param length the number of elements to include
1681      * @return true if the parameters are valid and designate a subarray of positive length
1682      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
1683      * @since 3.3
1684      */
1685     public static boolean verifyValues(
1686         final double[] values,
1687         final double[] weights,
1688         final int begin,
1689         final int length) throws MathIllegalArgumentException {
1690         return verifyValues(values, weights, begin, length, false);
1691     }
1692 
1693     /**
1694      * This method is used
1695      * to verify that the begin and length parameters designate a subarray of positive length
1696      * and the weights are all non-negative, non-NaN, finite, and not all zero.
1697      * <p>
1698      * <ul>
1699      * <li>returns <code>true</code> iff the parameters designate a subarray of
1700      * non-negative length and the weights array contains legitimate values.</li>
1701      * <li>throws <code>MathIllegalArgumentException</code> if any of the following are true:
1702      * <ul><li>the values array is null</li>
1703      *     <li>the weights array is null</li>
1704      *     <li>the weights array does not have the same length as the values array</li>
1705      *     <li>the weights array contains one or more infinite values</li>
1706      *     <li>the weights array contains one or more NaN values</li>
1707      *     <li>the weights array contains negative values</li>
1708      *     <li>the start and length arguments do not determine a valid array</li></ul>
1709      * </li>
1710      * <li>returns <code>false</li> if the array is non-null, but
1711      * <code>length</code> is 0 unless <code>allowEmpty</code> is <code>true</code>.
1712      * </ul></p>
1713      *
1714      * @param values the input array.
1715      * @param weights the weights array.
1716      * @param begin index of the first array element to include.
1717      * @param length the number of elements to include.
1718      * @param allowEmpty if {@code true} than allow zero length arrays to pass.
1719      * @return {@code true} if the parameters are valid.
1720      * @throws NullArgumentException if either of the arrays are null
1721      * @throws MathIllegalArgumentException if the array indices are not valid,
1722      * the weights array contains NaN, infinite or negative elements, or there
1723      * are no positive weights.
1724      * @since 3.3
1725      */
1726     public static boolean verifyValues(final double[] values, final double[] weights,
1727             final int begin, final int length, final boolean allowEmpty) throws MathIllegalArgumentException {
1728 
1729         if (weights == null || values == null) {
1730             throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
1731         }
1732 
1733         if (weights.length != values.length) {
1734             throw new DimensionMismatchException(weights.length, values.length);
1735         }
1736 
1737         boolean containsPositiveWeight = false;
1738         for (int i = begin; i < begin + length; i++) {
1739             final double weight = weights[i];
1740             if (Double.isNaN(weight)) {
1741                 throw new MathIllegalArgumentException(LocalizedFormats.NAN_ELEMENT_AT_INDEX, Integer.valueOf(i));
1742             }
1743             if (Double.isInfinite(weight)) {
1744                 throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, Double.valueOf(weight), Integer.valueOf(i));
1745             }
1746             if (weight < 0) {
1747                 throw new MathIllegalArgumentException(LocalizedFormats.NEGATIVE_ELEMENT_AT_INDEX, Integer.valueOf(i), Double.valueOf(weight));
1748             }
1749             if (!containsPositiveWeight && weight > 0.0) {
1750                 containsPositiveWeight = true;
1751             }
1752         }
1753 
1754         if (!containsPositiveWeight) {
1755             throw new MathIllegalArgumentException(LocalizedFormats.WEIGHT_AT_LEAST_ONE_NON_ZERO);
1756         }
1757 
1758         return verifyValues(values, begin, length, allowEmpty);
1759     }
1760 }