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