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