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