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