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.Iterator;
026import java.util.List;
027import java.util.TreeSet;
028
029import org.apache.commons.math3.Field;
030import org.apache.commons.math3.random.RandomGenerator;
031import org.apache.commons.math3.random.Well19937c;
032import org.apache.commons.math3.distribution.UniformIntegerDistribution;
033import org.apache.commons.math3.exception.DimensionMismatchException;
034import org.apache.commons.math3.exception.MathArithmeticException;
035import org.apache.commons.math3.exception.MathIllegalArgumentException;
036import org.apache.commons.math3.exception.MathInternalError;
037import org.apache.commons.math3.exception.NoDataException;
038import org.apache.commons.math3.exception.NonMonotonicSequenceException;
039import org.apache.commons.math3.exception.NotPositiveException;
040import org.apache.commons.math3.exception.NotStrictlyPositiveException;
041import org.apache.commons.math3.exception.NullArgumentException;
042import org.apache.commons.math3.exception.NumberIsTooLargeException;
043import org.apache.commons.math3.exception.NotANumberException;
044import org.apache.commons.math3.exception.util.LocalizedFormats;
045
046/**
047 * Arrays utilities.
048 *
049 * @since 3.0
050 */
051public class MathArrays {
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        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}