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