    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
    018    package org.apache.commons.math3.util;
    019
    020    import java.util.List;
    021    import java.util.ArrayList;
    022    import java.util.Comparator;
    023    import java.util.Collections;
    024
    025    import org.apache.commons.math3.exception.DimensionMismatchException;
    026    import org.apache.commons.math3.exception.MathInternalError;
    027    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
    028    import org.apache.commons.math3.exception.NotPositiveException;
    029    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    030    import org.apache.commons.math3.exception.NullArgumentException;
    031    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    032    import org.apache.commons.math3.exception.util.LocalizedFormats;
    033    import org.apache.commons.math3.exception.MathArithmeticException;
    034
    035    /**
    036     * Arrays utilities.
    037     *
    038     * @since 3.0
    039     * @version $Id: MathArrays.java 1422313 2012-12-15 18:53:41Z psteitz $
    040     */
    041    public class MathArrays {
    042        /** Factor used for splitting double numbers: n = 2^27 + 1 (i.e. {@value}). */
    043        private static final int SPLIT_FACTOR = 0x8000001;
    044
    045        /**
    046         * Private constructor.
    047         */
    048        private MathArrays() {}
    049
    050        /**
    051         * Real-valued function that operate on an array or a part of it.
    052         * @since 3.1
    053         */
    054        public interface Function {
    055            /**
    056             * Operates on an entire array.
    057             *
    058             * @param array Array to operate on.
    059             * @return the result of the operation.
    060             */
    061            double evaluate(double[] array);
    062            /**
    063             * @param array Array to operate on.
    064             * @param startIndex Index of the first element to take into account.
    065             * @param numElements Number of elements to take into account.
    066             * @return the result of the operation.
    067             */
    068            double evaluate(double[] array,
    069                            int startIndex,
    070                            int numElements);
    071        }
    072
    073        /**
    074         * Creates an array whose contents will be the element-by-element
    075         * addition of the arguments.
    076         *
    077         * @param a First term of the addition.
    078         * @param b Second term of the addition.
    079         * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}.
    080         * @throws DimensionMismatchException if the array lengths differ.
    081         * @since 3.1
    082         */
    083        public static double[] ebeAdd(double[] a,
    084                                      double[] b) {
    085            if (a.length != b.length) {
    086                throw new DimensionMismatchException(a.length, b.length);
    087            }
    088
    089            final double[] result = a.clone();
    090            for (int i = 0; i < a.length; i++) {
    091                result[i] += b[i];
    092            }
    093            return result;
    094        }
    095        /**
    096         * Creates an array whose contents will be the element-by-element
    097         * subtraction of the second argument from the first.
    098         *
    099         * @param a First term.
    100         * @param b Element to be subtracted.
    101         * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}.
    102         * @throws DimensionMismatchException if the array lengths differ.
    103         * @since 3.1
    104         */
    105        public static double[] ebeSubtract(double[] a,
    106                                           double[] b) {
    107            if (a.length != b.length) {
    108                throw new DimensionMismatchException(a.length, b.length);
    109            }
    110
    111            final double[] result = a.clone();
    112            for (int i = 0; i < a.length; i++) {
    113                result[i] -= b[i];
    114            }
    115            return result;
    116        }
    117        /**
    118         * Creates an array whose contents will be the element-by-element
    119         * multiplication of the arguments.
    120         *
    121         * @param a First factor of the multiplication.
    122         * @param b Second factor of the multiplication.
    123         * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}.
    124         * @throws DimensionMismatchException if the array lengths differ.
    125         * @since 3.1
    126         */
    127        public static double[] ebeMultiply(double[] a,
    128                                           double[] b) {
    129            if (a.length != b.length) {
    130                throw new DimensionMismatchException(a.length, b.length);
    131            }
    132
    133            final double[] result = a.clone();
    134            for (int i = 0; i < a.length; i++) {
    135                result[i] *= b[i];
    136            }
    137            return result;
    138        }
    139        /**
    140         * Creates an array whose contents will be the element-by-element
    141         * division of the first argument by the second.
    142         *
    143         * @param a Numerator of the division.
    144         * @param b Denominator of the division.
    145         * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}.
    146         * @throws DimensionMismatchException if the array lengths differ.
    147         * @since 3.1
    148         */
    149        public static double[] ebeDivide(double[] a,
    150                                         double[] b) {
    151            if (a.length != b.length) {
    152                throw new DimensionMismatchException(a.length, b.length);
    153            }
    154
    155            final double[] result = a.clone();
    156            for (int i = 0; i < a.length; i++) {
    157                result[i] /= b[i];
    158            }
    159            return result;
    160        }
    161
    162        /**
    163         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
    164         *
    165         * @param p1 the first point
    166         * @param p2 the second point
    167         * @return the L<sub>1</sub> distance between the two points
    168         */
    169        public static double distance1(double[] p1, double[] p2) {
    170            double sum = 0;
    171            for (int i = 0; i < p1.length; i++) {
    172                sum += FastMath.abs(p1[i] - p2[i]);
    173            }
    174            return sum;
    175        }
    176
    177        /**
    178         * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
    179         *
    180         * @param p1 the first point
    181         * @param p2 the second point
    182         * @return the L<sub>1</sub> distance between the two points
    183         */
    184        public static int distance1(int[] p1, int[] p2) {
    185          int sum = 0;
    186          for (int i = 0; i < p1.length; i++) {
    187              sum += FastMath.abs(p1[i] - p2[i]);
    188          }
    189          return sum;
    190        }
    191
    192        /**
    193         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
    194         *
    195         * @param p1 the first point
    196         * @param p2 the second point
    197         * @return the L<sub>2</sub> distance between the two points
    198         */
    199        public static double distance(double[] p1, double[] p2) {
    200            double sum = 0;
    201            for (int i = 0; i < p1.length; i++) {
    202                final double dp = p1[i] - p2[i];
    203                sum += dp * dp;
    204            }
    205            return FastMath.sqrt(sum);
    206        }
    207
    208        /**
    209         * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
    210         *
    211         * @param p1 the first point
    212         * @param p2 the second point
    213         * @return the L<sub>2</sub> distance between the two points
    214         */
    215        public static double distance(int[] p1, int[] p2) {
    216          double sum = 0;
    217          for (int i = 0; i < p1.length; i++) {
    218              final double dp = p1[i] - p2[i];
    219              sum += dp * dp;
    220          }
    221          return FastMath.sqrt(sum);
    222        }
    223
    224        /**
    225         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
    226         *
    227         * @param p1 the first point
    228         * @param p2 the second point
    229         * @return the L<sub>&infin;</sub> distance between the two points
    230         */
    231        public static double distanceInf(double[] p1, double[] p2) {
    232            double max = 0;
    233            for (int i = 0; i < p1.length; i++) {
    234                max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
    235            }
    236            return max;
    237        }
    238
    239        /**
    240         * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
    241         *
    242         * @param p1 the first point
    243         * @param p2 the second point
    244         * @return the L<sub>&infin;</sub> distance between the two points
    245         */
    246        public static int distanceInf(int[] p1, int[] p2) {
    247            int max = 0;
    248            for (int i = 0; i < p1.length; i++) {
    249                max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
    250            }
    251            return max;
    252        }
    253
    254        /**
    255         * Specification of ordering direction.
    256         */
    257        public static enum OrderDirection {
    258            /** Constant for increasing direction. */
    259            INCREASING,
    260            /** Constant for decreasing direction. */
    261            DECREASING
    262        }
    263
    264        /**
    265         * Check that an array is monotonically increasing or decreasing.
    266         *
    267         * @param <T> the type of the elements in the specified array
    268         * @param val Values.
    269         * @param dir Ordering direction.
    270         * @param strict Whether the order should be strict.
    271         * @return {@code true} if sorted, {@code false} otherwise.
    272         */
    273        public static  <T extends Comparable<? super T>> boolean isMonotonic(T[] val,
    274                                          OrderDirection dir,
    275                                          boolean strict) {
    276            T previous = val[0];
    277            final int max = val.length;
    278            for (int i = 1; i < max; i++) {
    279                final int comp;
    280                switch (dir) {
    281                case INCREASING:
    282                    comp = previous.compareTo(val[i]);
    283                    if (strict) {
    284                        if (comp >= 0) {
    285                            return false;
    286                        }
    287                    } else {
    288                        if (comp > 0) {
    289                            return false;
    290                        }
    291                    }
    292                    break;
    293                case DECREASING:
    294                    comp = val[i].compareTo(previous);
    295                    if (strict) {
    296                        if (comp >= 0) {
    297                            return false;
    298                        }
    299                    } else {
    300                        if (comp > 0) {
    301                           return false;
    302                        }
    303                    }
    304                    break;
    305                default:
    306                    // Should never happen.
    307                    throw new MathInternalError();
    308                }
    309
    310                previous = val[i];
    311            }
    312            return true;
    313        }
    314
    315        /**
    316         * Check that an array is monotonically increasing or decreasing.
    317         *
    318         * @param val Values.
    319         * @param dir Ordering direction.
    320         * @param strict Whether the order should be strict.
    321         * @return {@code true} if sorted, {@code false} otherwise.
    322         */
    323        public static boolean isMonotonic(double[] val,
    324                                          OrderDirection dir,
    325                                          boolean strict) {
    326            return checkOrder(val, dir, strict, false);
    327        }
    328
    329        /**
    330         * Check that the given array is sorted.
    331         *
    332         * @param val Values.
    333         * @param dir Ordering direction.
    334         * @param strict Whether the order should be strict.
    335         * @param abort Whether to throw an exception if the check fails.
    336         * @return {@code true} if the array is sorted.
    337         * @throws NonMonotonicSequenceException if the array is not sorted
    338         * and {@code abort} is {@code true}.
    339         */
    340        public static boolean checkOrder(double[] val, OrderDirection dir,
    341                                         boolean strict, boolean abort)
    342            throws NonMonotonicSequenceException {
    343            double previous = val[0];
    344            final int max = val.length;
    345
    346            int index;
    347            ITEM:
    348            for (index = 1; index < max; index++) {
    349                switch (dir) {
    350                case INCREASING:
    351                    if (strict) {
    352                        if (val[index] <= previous) {
    353                            break ITEM;
    354                        }
    355                    } else {
    356                        if (val[index] < previous) {
    357                            break ITEM;
    358                        }
    359                    }
    360                    break;
    361                case DECREASING:
    362                    if (strict) {
    363                        if (val[index] >= previous) {
    364                            break ITEM;
    365                        }
    366                    } else {
    367                        if (val[index] > previous) {
    368                            break ITEM;
    369                        }
    370                    }
    371                    break;
    372                default:
    373                    // Should never happen.
    374                    throw new MathInternalError();
    375                }
    376
    377                previous = val[index];
    378            }
    379
    380            if (index == max) {
    381                // Loop completed.
    382                return true;
    383            }
    384
    385            // Loop early exit means wrong ordering.
    386            if (abort) {
    387                throw new NonMonotonicSequenceException(val[index], previous, index, dir, strict);
    388            } else {
    389                return false;
    390            }
    391        }
    392
    393        /**
    394         * Check that the given array is sorted.
    395         *
    396         * @param val Values.
    397         * @param dir Ordering direction.
    398         * @param strict Whether the order should be strict.
    399         * @throws NonMonotonicSequenceException if the array is not sorted.
    400         * @since 2.2
    401         */
    402        public static void checkOrder(double[] val, OrderDirection dir,
    403                                      boolean strict) throws NonMonotonicSequenceException {
    404            checkOrder(val, dir, strict, true);
    405        }
    406
    407        /**
    408         * Check that the given array is sorted in strictly increasing order.
    409         *
    410         * @param val Values.
    411         * @throws NonMonotonicSequenceException if the array is not sorted.
    412         * @since 2.2
    413         */
    414        public static void checkOrder(double[] val) throws NonMonotonicSequenceException {
    415            checkOrder(val, OrderDirection.INCREASING, true);
    416        }
    417
    418        /**
    419         * Throws DimensionMismatchException if the input array is not rectangular.
    420         *
    421         * @param in array to be tested
    422         * @throws NullArgumentException if input array is null
    423         * @throws DimensionMismatchException if input array is not rectangular
    424         * @since 3.1
    425         */
    426        public static void checkRectangular(final long[][] in)
    427            throws NullArgumentException, DimensionMismatchException {
    428            MathUtils.checkNotNull(in);
    429            for (int i = 1; i < in.length; i++) {
    430                if (in[i].length != in[0].length) {
    431                    throw new DimensionMismatchException(
    432                            LocalizedFormats.DIFFERENT_ROWS_LENGTHS,
    433                            in[i].length, in[0].length);
    434                }
    435            }
    436        }
    437
    438        /**
    439         * Check that all entries of the input array are strictly positive.
    440         *
    441         * @param in Array to be tested
    442         * @throws NotStrictlyPositiveException if any entries of the array are not
    443         * strictly positive.
    444         * @since 3.1
    445         */
    446        public static void checkPositive(final double[] in)
    447            throws NotStrictlyPositiveException {
    448            for (int i = 0; i < in.length; i++) {
    449                if (in[i] <= 0) {
    450                    throw new NotStrictlyPositiveException(in[i]);
    451                }
    452            }
    453        }
    454
    455        /**
    456         * Check that all entries of the input array are >= 0.
    457         *
    458         * @param in Array to be tested
    459         * @throws NotPositiveException if any array entries are less than 0.
    460         * @since 3.1
    461         */
    462        public static void checkNonNegative(final long[] in)
    463            throws NotPositiveException {
    464            for (int i = 0; i < in.length; i++) {
    465                if (in[i] < 0) {
    466                    throw new NotPositiveException(in[i]);
    467                }
    468            }
    469        }
    470
    471        /**
    472         * Check all entries of the input array are >= 0.
    473         *
    474         * @param in Array to be tested
    475         * @throws NotPositiveException if any array entries are less than 0.
    476         * @since 3.1
    477         */
    478        public static void checkNonNegative(final long[][] in)
    479            throws NotPositiveException {
    480            for (int i = 0; i < in.length; i ++) {
    481                for (int j = 0; j < in[i].length; j++) {
    482                    if (in[i][j] < 0) {
    483                        throw new NotPositiveException(in[i][j]);
    484                    }
    485                }
    486            }
    487        }
    488
    489        /**
    490         * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
    491         * Translation of the minpack enorm subroutine.
    492         *
    493         * The redistribution policy for MINPACK is available
    494         * <a href="http://www.netlib.org/minpack/disclaimer">here</a>, for
    495         * convenience, it is reproduced below.</p>
    496         *
    497         * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
    498         * <tr><td>
    499         *    Minpack Copyright Notice (1999) University of Chicago.
    500         *    All rights reserved
    501         * </td></tr>
    502         * <tr><td>
    503         * Redistribution and use in source and binary forms, with or without
    504         * modification, are permitted provided that the following conditions
    505         * are met:
    506         * <ol>
    507         *  <li>Redistributions of source code must retain the above copyright
    508         *      notice, this list of conditions and the following disclaimer.</li>
    509         * <li>Redistributions in binary form must reproduce the above
    510         *     copyright notice, this list of conditions and the following
    511         *     disclaimer in the documentation and/or other materials provided
    512         *     with the distribution.</li>
    513         * <li>The end-user documentation included with the redistribution, if any,
    514         *     must include the following acknowledgment:
    515         *     {@code This product includes software developed by the University of
    516         *           Chicago, as Operator of Argonne National Laboratory.}
    517         *     Alternately, this acknowledgment may appear in the software itself,
    518         *     if and wherever such third-party acknowledgments normally appear.</li>
    519         * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
    520         *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
    521         *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
    522         *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
    523         *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
    524         *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
    525         *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
    526         *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
    527         *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
    528         *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
    529         *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
    530         *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
    531         *     BE CORRECTED.</strong></li>
    532         * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
    533         *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
    534         *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
    535         *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
    536         *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
    537         *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
    538         *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
    539         *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
    540         *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
    541         *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
    542         * <ol></td></tr>
    543         * </table>
    544         *
    545         * @param v Vector of doubles.
    546         * @return the 2-norm of the vector.
    547         * @since 2.2
    548         */
    549        public static double safeNorm(double[] v) {
    550            double rdwarf = 3.834e-20;
    551            double rgiant = 1.304e+19;
    552            double s1 = 0;
    553            double s2 = 0;
    554            double s3 = 0;
    555            double x1max = 0;
    556            double x3max = 0;
    557            double floatn = v.length;
    558            double agiant = rgiant / floatn;
    559            for (int i = 0; i < v.length; i++) {
    560                double xabs = Math.abs(v[i]);
    561                if (xabs < rdwarf || xabs > agiant) {
    562                    if (xabs > rdwarf) {
    563                        if (xabs > x1max) {
    564                            double r = x1max / xabs;
    565                            s1= 1 + s1 * r * r;
    566                            x1max = xabs;
    567                        } else {
    568                            double r = xabs / x1max;
    569                            s1 += r * r;
    570                        }
    571                    } else {
    572                        if (xabs > x3max) {
    573                            double r = x3max / xabs;
    574                            s3= 1 + s3 * r * r;
    575                            x3max = xabs;
    576                        } else {
    577                            if (xabs != 0) {
    578                                double r = xabs / x3max;
    579                                s3 += r * r;
    580                            }
    581                        }
    582                    }
    583                } else {
    584                    s2 += xabs * xabs;
    585                }
    586            }
    587            double norm;
    588            if (s1 != 0) {
    589                norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
    590            } else {
    591                if (s2 == 0) {
    592                    norm = x3max * Math.sqrt(s3);
    593                } else {
    594                    if (s2 >= x3max) {
    595                        norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
    596                    } else {
    597                        norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
    598                    }
    599                }
    600            }
    601            return norm;
    602        }
    603
    604        /**
    605         * Sort an array in ascending order in place and perform the same reordering
    606         * of entries on other arrays. For example, if
    607         * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then
    608         * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]},
    609         * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}.
    610         *
    611         * @param x Array to be sorted and used as a pattern for permutation
    612         * of the other arrays.
    613         * @param yList Set of arrays whose permutations of entries will follow
    614         * those performed on {@code x}.
    615         * @throws DimensionMismatchException if any {@code y} is not the same
    616         * size as {@code x}.
    617         * @throws NullArgumentException if {@code x} or any {@code y} is null.
    618         * @since 3.0
    619         */
    620        public static void sortInPlace(double[] x, double[] ... yList)
    621            throws DimensionMismatchException, NullArgumentException {
    622            sortInPlace(x, OrderDirection.INCREASING, yList);
    623        }
    624
    625        /**
    626         * Sort an array in place and perform the same reordering of entries on
    627         * other arrays.  This method works the same as the other
    628         * {@link #sortInPlace(double[], double[][]) sortInPlace} method, but
    629         * allows the order of the sort to be provided in the {@code dir}
    630         * parameter.
    631         *
    632         * @param x Array to be sorted and used as a pattern for permutation
    633         * of the other arrays.
    634         * @param dir Order direction.
    635         * @param yList Set of arrays whose permutations of entries will follow
    636         * those performed on {@code x}.
    637         * @throws DimensionMismatchException if any {@code y} is not the same
    638         * size as {@code x}.
    639         * @throws NullArgumentException if {@code x} or any {@code y} is null
    640         * @since 3.0
    641         */
    642        public static void sortInPlace(double[] x,
    643                                       final OrderDirection dir,
    644                                       double[] ... yList)
    645            throws NullArgumentException, DimensionMismatchException {
    646            if (x == null) {
    647                throw new NullArgumentException();
    648            }
    649
    650            final int len = x.length;
    651            final List<Pair<Double, double[]>> list
    652                = new ArrayList<Pair<Double, double[]>>(len);
    653
    654            final int yListLen = yList.length;
    655            for (int i = 0; i < len; i++) {
    656                final double[] yValues = new double[yListLen];
    657                for (int j = 0; j < yListLen; j++) {
    658                    double[] y = yList[j];
    659                    if (y == null) {
    660                        throw new NullArgumentException();
    661                    }
    662                    if (y.length != len) {
    663                        throw new DimensionMismatchException(y.length, len);
    664                    }
    665                    yValues[j] = y[i];
    666                }
    667                list.add(new Pair<Double, double[]>(x[i], yValues));
    668            }
    669
    670            final Comparator<Pair<Double, double[]>> comp
    671                = new Comparator<Pair<Double, double[]>>() {
    672                public int compare(Pair<Double, double[]> o1,
    673                                   Pair<Double, double[]> o2) {
    674                    int val;
    675                    switch (dir) {
    676                    case INCREASING:
    677                        val = o1.getKey().compareTo(o2.getKey());
    678                    break;
    679                    case DECREASING:
    680                        val = o2.getKey().compareTo(o1.getKey());
    681                    break;
    682                    default:
    683                        // Should never happen.
    684                        throw new MathInternalError();
    685                    }
    686                    return val;
    687                }
    688            };
    689
    690            Collections.sort(list, comp);
    691
    692            for (int i = 0; i < len; i++) {
    693                final Pair<Double, double[]> e = list.get(i);
    694                x[i] = e.getKey();
    695                final double[] yValues = e.getValue();
    696                for (int j = 0; j < yListLen; j++) {
    697                    yList[j][i] = yValues[j];
    698                }
    699            }
    700        }
    701
    702        /**
    703         * Creates a copy of the {@code source} array.
    704         *
    705         * @param source Array to be copied.
    706         * @return the copied array.
    707         */
    708         public static int[] copyOf(int[] source) {
    709             return copyOf(source, source.length);
    710         }
    711
    712        /**
    713         * Creates a copy of the {@code source} array.
    714         *
    715         * @param source Array to be copied.
    716         * @return the copied array.
    717         */
    718         public static double[] copyOf(double[] source) {
    719             return copyOf(source, source.length);
    720         }
    721
    722        /**
    723         * Creates a copy of the {@code source} array.
    724         *
    725         * @param source Array to be copied.
    726         * @param len Number of entries to copy. If smaller then the source
    727         * length, the copy will be truncated, if larger it will padded with
    728         * zeroes.
    729         * @return the copied array.
    730         */
    731        public static int[] copyOf(int[] source, int len) {
    732             final int[] output = new int[len];
    733             System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
    734             return output;
    735         }
    736
    737        /**
    738         * Creates a copy of the {@code source} array.
    739         *
    740         * @param source Array to be copied.
    741         * @param len Number of entries to copy. If smaller then the source
    742         * length, the copy will be truncated, if larger it will padded with
    743         * zeroes.
    744         * @return the copied array.
    745         */
    746        public static double[] copyOf(double[] source, int len) {
    747             final double[] output = new double[len];
    748             System.arraycopy(source, 0, output, 0, FastMath.min(len, source.length));
    749             return output;
    750         }
    751
    752        /**
    753         * Compute a linear combination accurately.
    754         * This method computes the sum of the products
    755         * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
    756         * It does so by using specific multiplication and addition algorithms to
    757         * preserve accuracy and reduce cancellation effects.
    758         * <br/>
    759         * It is based on the 2005 paper
    760         * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
    761         * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
    762         * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
    763         *
    764         * @param a Factors.
    765         * @param b Factors.
    766         * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
    767         * @throws DimensionMismatchException if arrays dimensions don't match
    768         */
    769        public static double linearCombination(final double[] a, final double[] b)
    770            throws DimensionMismatchException {
    771            final int len = a.length;
    772            if (len != b.length) {
    773                throw new DimensionMismatchException(len, b.length);
    774            }
    775
    776            final double[] prodHigh = new double[len];
    777            double prodLowSum = 0;
    778
    779            for (int i = 0; i < len; i++) {
    780                final double ai = a[i];
    781                final double ca = SPLIT_FACTOR * ai;
    782                final double aHigh = ca - (ca - ai);
    783                final double aLow = ai - aHigh;
    784
    785                final double bi = b[i];
    786                final double cb = SPLIT_FACTOR * bi;
    787                final double bHigh = cb - (cb - bi);
    788                final double bLow = bi - bHigh;
    789                prodHigh[i] = ai * bi;
    790                final double prodLow = aLow * bLow - (((prodHigh[i] -
    791                                                        aHigh * bHigh) -
    792                                                       aLow * bHigh) -
    793                                                      aHigh * bLow);
    794                prodLowSum += prodLow;
    795            }
    796
    797
    798            final double prodHighCur = prodHigh[0];
    799            double prodHighNext = prodHigh[1];
    800            double sHighPrev = prodHighCur + prodHighNext;
    801            double sPrime = sHighPrev - prodHighNext;
    802            double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime);
    803
    804            final int lenMinusOne = len - 1;
    805            for (int i = 1; i < lenMinusOne; i++) {
    806                prodHighNext = prodHigh[i + 1];
    807                final double sHighCur = sHighPrev + prodHighNext;
    808                sPrime = sHighCur - prodHighNext;
    809                sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime);
    810                sHighPrev = sHighCur;
    811            }
    812
    813            double result = sHighPrev + (prodLowSum + sLowSum);
    814
    815            if (Double.isNaN(result)) {
    816                // either we have split infinite numbers or some coefficients were NaNs,
    817                // just rely on the naive implementation and let IEEE754 handle this
    818                result = 0;
    819                for (int i = 0; i < len; ++i) {
    820                    result += a[i] * b[i];
    821                }
    822            }
    823
    824            return result;
    825        }
    826
    827        /**
    828         * Compute a linear combination accurately.
    829         * <p>
    830         * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
    831         * a<sub>2</sub>&times;b<sub>2</sub> to high accuracy. It does
    832         * so by using specific multiplication and addition algorithms to
    833         * preserve accuracy and reduce cancellation effects. It is based
    834         * on the 2005 paper <a
    835         * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
    836         * Accurate Sum and Dot Product</a> by Takeshi Ogita,
    837         * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
    838         * </p>
    839         * @param a1 first factor of the first term
    840         * @param b1 second factor of the first term
    841         * @param a2 first factor of the second term
    842         * @param b2 second factor of the second term
    843         * @return a<sub>1</sub>&times;b<sub>1</sub> +
    844         * a<sub>2</sub>&times;b<sub>2</sub>
    845         * @see #linearCombination(double, double, double, double, double, double)
    846         * @see #linearCombination(double, double, double, double, double, double, double, double)
    847         */
    848        public static double linearCombination(final double a1, final double b1,
    849                                               final double a2, final double b2) {
    850
    851            // the code below is split in many additions/subtractions that may
    852            // appear redundant. However, they should NOT be simplified, as they
    853            // use IEEE754 floating point arithmetic rounding properties.
    854            // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
    855            // The variable naming conventions are that xyzHigh contains the most significant
    856            // bits of xyz and xyzLow contains its least significant bits. So theoretically
    857            // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
    858            // be represented in only one double precision number so we preserve two numbers
    859            // to hold it as long as we can, combining the high and low order bits together
    860            // only at the end, after cancellation may have occurred on high order bits
    861
    862            // split a1 and b1 as two 26 bits numbers
    863            final double ca1        = SPLIT_FACTOR * a1;
    864            final double a1High     = ca1 - (ca1 - a1);
    865            final double a1Low      = a1 - a1High;
    866            final double cb1        = SPLIT_FACTOR * b1;
    867            final double b1High     = cb1 - (cb1 - b1);
    868            final double b1Low      = b1 - b1High;
    869
    870            // accurate multiplication a1 * b1
    871            final double prod1High  = a1 * b1;
    872            final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
    873
    874            // split a2 and b2 as two 26 bits numbers
    875            final double ca2        = SPLIT_FACTOR * a2;
    876            final double a2High     = ca2 - (ca2 - a2);
    877            final double a2Low      = a2 - a2High;
    878            final double cb2        = SPLIT_FACTOR * b2;
    879            final double b2High     = cb2 - (cb2 - b2);
    880            final double b2Low      = b2 - b2High;
    881
    882            // accurate multiplication a2 * b2
    883            final double prod2High  = a2 * b2;
    884            final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
    885
    886            // accurate addition a1 * b1 + a2 * b2
    887            final double s12High    = prod1High + prod2High;
    888            final double s12Prime   = s12High - prod2High;
    889            final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
    890
    891            // final rounding, s12 may have suffered many cancellations, we try
    892            // to recover some bits from the extra words we have saved up to now
    893            double result = s12High + (prod1Low + prod2Low + s12Low);
    894
    895            if (Double.isNaN(result)) {
    896                // either we have split infinite numbers or some coefficients were NaNs,
    897                // just rely on the naive implementation and let IEEE754 handle this
    898                result = a1 * b1 + a2 * b2;
    899            }
    900
    901            return result;
    902        }
    903
    904        /**
    905         * Compute a linear combination accurately.
    906         * <p>
    907         * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
    908         * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
    909         * to high accuracy. It does so by using specific multiplication and
    910         * addition algorithms to preserve accuracy and reduce cancellation effects.
    911         * It is based on the 2005 paper <a
    912         * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
    913         * Accurate Sum and Dot Product</a> by Takeshi Ogita,
    914         * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
    915         * </p>
    916         * @param a1 first factor of the first term
    917         * @param b1 second factor of the first term
    918         * @param a2 first factor of the second term
    919         * @param b2 second factor of the second term
    920         * @param a3 first factor of the third term
    921         * @param b3 second factor of the third term
    922         * @return a<sub>1</sub>&times;b<sub>1</sub> +
    923         * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
    924         * @see #linearCombination(double, double, double, double)
    925         * @see #linearCombination(double, double, double, double, double, double, double, double)
    926         */
    927        public static double linearCombination(final double a1, final double b1,
    928                                               final double a2, final double b2,
    929                                               final double a3, final double b3) {
    930
    931            // the code below is split in many additions/subtractions that may
    932            // appear redundant. However, they should NOT be simplified, as they
    933            // do use IEEE754 floating point arithmetic rounding properties.
    934            // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
    935            // The variables naming conventions are that xyzHigh contains the most significant
    936            // bits of xyz and xyzLow contains its least significant bits. So theoretically
    937            // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
    938            // be represented in only one double precision number so we preserve two numbers
    939            // to hold it as long as we can, combining the high and low order bits together
    940            // only at the end, after cancellation may have occurred on high order bits
    941
    942            // split a1 and b1 as two 26 bits numbers
    943            final double ca1        = SPLIT_FACTOR * a1;
    944            final double a1High     = ca1 - (ca1 - a1);
    945            final double a1Low      = a1 - a1High;
    946            final double cb1        = SPLIT_FACTOR * b1;
    947            final double b1High     = cb1 - (cb1 - b1);
    948            final double b1Low      = b1 - b1High;
    949
    950            // accurate multiplication a1 * b1
    951            final double prod1High  = a1 * b1;
    952            final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
    953
    954            // split a2 and b2 as two 26 bits numbers
    955            final double ca2        = SPLIT_FACTOR * a2;
    956            final double a2High     = ca2 - (ca2 - a2);
    957            final double a2Low      = a2 - a2High;
    958            final double cb2        = SPLIT_FACTOR * b2;
    959            final double b2High     = cb2 - (cb2 - b2);
    960            final double b2Low      = b2 - b2High;
    961
    962            // accurate multiplication a2 * b2
    963            final double prod2High  = a2 * b2;
    964            final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
    965
    966            // split a3 and b3 as two 26 bits numbers
    967            final double ca3        = SPLIT_FACTOR * a3;
    968            final double a3High     = ca3 - (ca3 - a3);
    969            final double a3Low      = a3 - a3High;
    970            final double cb3        = SPLIT_FACTOR * b3;
    971            final double b3High     = cb3 - (cb3 - b3);
    972            final double b3Low      = b3 - b3High;
    973
    974            // accurate multiplication a3 * b3
    975            final double prod3High  = a3 * b3;
    976            final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
    977
    978            // accurate addition a1 * b1 + a2 * b2
    979            final double s12High    = prod1High + prod2High;
    980            final double s12Prime   = s12High - prod2High;
    981            final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
    982
    983            // accurate addition a1 * b1 + a2 * b2 + a3 * b3
    984            final double s123High   = s12High + prod3High;
    985            final double s123Prime  = s123High - prod3High;
    986            final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
    987
    988            // final rounding, s123 may have suffered many cancellations, we try
    989            // to recover some bits from the extra words we have saved up to now
    990            double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low);
    991
    992            if (Double.isNaN(result)) {
    993                // either we have split infinite numbers or some coefficients were NaNs,
    994                // just rely on the naive implementation and let IEEE754 handle this
    995                result = a1 * b1 + a2 * b2 + a3 * b3;
    996            }
    997
    998            return result;
    999        }
    1000
    1001        /**
    1002         * Compute a linear combination accurately.
    1003         * <p>
    1004         * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
    1005         * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
    1006         * a<sub>4</sub>&times;b<sub>4</sub>
    1007         * to high accuracy. It does so by using specific multiplication and
    1008         * addition algorithms to preserve accuracy and reduce cancellation effects.
    1009         * It is based on the 2005 paper <a
    1010         * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
    1011         * Accurate Sum and Dot Product</a> by Takeshi Ogita,
    1012         * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
    1013         * </p>
    1014         * @param a1 first factor of the first term
    1015         * @param b1 second factor of the first term
    1016         * @param a2 first factor of the second term
    1017         * @param b2 second factor of the second term
    1018         * @param a3 first factor of the third term
    1019         * @param b3 second factor of the third term
    1020         * @param a4 first factor of the third term
    1021         * @param b4 second factor of the third term
    1022         * @return a<sub>1</sub>&times;b<sub>1</sub> +
    1023         * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
    1024         * a<sub>4</sub>&times;b<sub>4</sub>
    1025         * @see #linearCombination(double, double, double, double)
    1026         * @see #linearCombination(double, double, double, double, double, double)
    1027         */
    1028        public static double linearCombination(final double a1, final double b1,
    1029                                               final double a2, final double b2,
    1030                                               final double a3, final double b3,
    1031                                               final double a4, final double b4) {
    1032
    1033            // the code below is split in many additions/subtractions that may
    1034            // appear redundant. However, they should NOT be simplified, as they
    1035            // do use IEEE754 floating point arithmetic rounding properties.
    1036            // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same as "a1"
    1037            // The variables naming conventions are that xyzHigh contains the most significant
    1038            // bits of xyz and xyzLow contains its least significant bits. So theoretically
    1039            // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
    1040            // be represented in only one double precision number so we preserve two numbers
    1041            // to hold it as long as we can, combining the high and low order bits together
    1042            // only at the end, after cancellation may have occurred on high order bits
    1043
    1044            // split a1 and b1 as two 26 bits numbers
    1045            final double ca1        = SPLIT_FACTOR * a1;
    1046            final double a1High     = ca1 - (ca1 - a1);
    1047            final double a1Low      = a1 - a1High;
    1048            final double cb1        = SPLIT_FACTOR * b1;
    1049            final double b1High     = cb1 - (cb1 - b1);
    1050            final double b1Low      = b1 - b1High;
    1051
    1052            // accurate multiplication a1 * b1
    1053            final double prod1High  = a1 * b1;
    1054            final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);
    1055
    1056            // split a2 and b2 as two 26 bits numbers
    1057            final double ca2        = SPLIT_FACTOR * a2;
    1058            final double a2High     = ca2 - (ca2 - a2);
    1059            final double a2Low      = a2 - a2High;
    1060            final double cb2        = SPLIT_FACTOR * b2;
    1061            final double b2High     = cb2 - (cb2 - b2);
    1062            final double b2Low      = b2 - b2High;
    1063
    1064            // accurate multiplication a2 * b2
    1065            final double prod2High  = a2 * b2;
    1066            final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);
    1067
    1068            // split a3 and b3 as two 26 bits numbers
    1069            final double ca3        = SPLIT_FACTOR * a3;
    1070            final double a3High     = ca3 - (ca3 - a3);
    1071            final double a3Low      = a3 - a3High;
    1072            final double cb3        = SPLIT_FACTOR * b3;
    1073            final double b3High     = cb3 - (cb3 - b3);
    1074            final double b3Low      = b3 - b3High;
    1075
    1076            // accurate multiplication a3 * b3
    1077            final double prod3High  = a3 * b3;
    1078            final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);
    1079
    1080            // split a4 and b4 as two 26 bits numbers
    1081            final double ca4        = SPLIT_FACTOR * a4;
    1082            final double a4High     = ca4 - (ca4 - a4);
    1083            final double a4Low      = a4 - a4High;
    1084            final double cb4        = SPLIT_FACTOR * b4;
    1085            final double b4High     = cb4 - (cb4 - b4);
    1086            final double b4Low      = b4 - b4High;
    1087
    1088            // accurate multiplication a4 * b4
    1089            final double prod4High  = a4 * b4;
    1090            final double prod4Low   = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low);
    1091
    1092            // accurate addition a1 * b1 + a2 * b2
    1093            final double s12High    = prod1High + prod2High;
    1094            final double s12Prime   = s12High - prod2High;
    1095            final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);
    1096
    1097            // accurate addition a1 * b1 + a2 * b2 + a3 * b3
    1098            final double s123High   = s12High + prod3High;
    1099            final double s123Prime  = s123High - prod3High;
    1100            final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);
    1101
    1102            // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4
    1103            final double s1234High  = s123High + prod4High;
    1104            final double s1234Prime = s1234High - prod4High;
    1105            final double s1234Low   = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime);
    1106
    1107            // final rounding, s1234 may have suffered many cancellations, we try
    1108            // to recover some bits from the extra words we have saved up to now
    1109            double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low);
    1110
    1111            if (Double.isNaN(result)) {
    1112                // either we have split infinite numbers or some coefficients were NaNs,
    1113                // just rely on the naive implementation and let IEEE754 handle this
    1114                result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
    1115            }
    1116
    1117            return result;
    1118        }
    1119
    1120        /**
    1121         * Returns true iff both arguments are null or have same dimensions and all
    1122         * their elements are equal as defined by
    1123         * {@link Precision#equals(float,float)}.
    1124         *
    1125         * @param x first array
    1126         * @param y second array
    1127         * @return true if the values are both null or have same dimension
    1128         * and equal elements.
    1129         */
    1130        public static boolean equals(float[] x, float[] y) {
    1131            if ((x == null) || (y == null)) {
    1132                return !((x == null) ^ (y == null));
    1133            }
    1134            if (x.length != y.length) {
    1135                return false;
    1136            }
    1137            for (int i = 0; i < x.length; ++i) {
    1138                if (!Precision.equals(x[i], y[i])) {
    1139                    return false;
    1140                }
    1141            }
    1142            return true;
    1143        }
    1144
    1145        /**
    1146         * Returns true iff both arguments are null or have same dimensions and all
    1147         * their elements are equal as defined by
    1148         * {@link Precision#equalsIncludingNaN(double,double) this method}.
    1149         *
    1150         * @param x first array
    1151         * @param y second array
    1152         * @return true if the values are both null or have same dimension and
    1153         * equal elements
    1154         * @since 2.2
    1155         */
    1156        public static boolean equalsIncludingNaN(float[] x, float[] y) {
    1157            if ((x == null) || (y == null)) {
    1158                return !((x == null) ^ (y == null));
    1159            }
    1160            if (x.length != y.length) {
    1161                return false;
    1162            }
    1163            for (int i = 0; i < x.length; ++i) {
    1164                if (!Precision.equalsIncludingNaN(x[i], y[i])) {
    1165                    return false;
    1166                }
    1167            }
    1168            return true;
    1169        }
    1170
    1171        /**
    1172         * Returns {@code true} iff both arguments are {@code null} or have same
    1173         * dimensions and all their elements are equal as defined by
    1174         * {@link Precision#equals(double,double)}.
    1175         *
    1176         * @param x First array.
    1177         * @param y Second array.
    1178         * @return {@code true} if the values are both {@code null} or have same
    1179         * dimension and equal elements.
    1180         */
    1181        public static boolean equals(double[] x, double[] y) {
    1182            if ((x == null) || (y == null)) {
    1183                return !((x == null) ^ (y == null));
    1184            }
    1185            if (x.length != y.length) {
    1186                return false;
    1187            }
    1188            for (int i = 0; i < x.length; ++i) {
    1189                if (!Precision.equals(x[i], y[i])) {
    1190                    return false;
    1191                }
    1192            }
    1193            return true;
    1194        }
    1195
    1196        /**
    1197         * Returns {@code true} iff both arguments are {@code null} or have same
    1198         * dimensions and all their elements are equal as defined by
    1199         * {@link Precision#equalsIncludingNaN(double,double) this method}.
    1200         *
    1201         * @param x First array.
    1202         * @param y Second array.
    1203         * @return {@code true} if the values are both {@code null} or have same
    1204         * dimension and equal elements.
    1205         * @since 2.2
    1206         */
    1207        public static boolean equalsIncludingNaN(double[] x, double[] y) {
    1208            if ((x == null) || (y == null)) {
    1209                return !((x == null) ^ (y == null));
    1210            }
    1211            if (x.length != y.length) {
    1212                return false;
    1213            }
    1214            for (int i = 0; i < x.length; ++i) {
    1215                if (!Precision.equalsIncludingNaN(x[i], y[i])) {
    1216                    return false;
    1217                }
    1218            }
    1219            return true;
    1220        }
    1221
    1222         /**
    1223          * Normalizes an array to make it sum to a specified value.
    1224          * Returns the result of the transformation <pre>
    1225          *    x |-> x * normalizedSum / sum
    1226          * </pre>
    1227          * applied to each non-NaN element x of the input array, where sum is the
    1228          * sum of the non-NaN entries in the input array.</p>
    1229          *
    1230          * <p>Throws IllegalArgumentException if {@code normalizedSum} is infinite
    1231          * or NaN and ArithmeticException if the input array contains any infinite elements
    1232          * or sums to 0.</p>
    1233          *
    1234          * <p>Ignores (i.e., copies unchanged to the output array) NaNs in the input array.</p>
    1235          *
    1236          * @param values Input array to be normalized
    1237          * @param normalizedSum Target sum for the normalized array
    1238          * @return the normalized array.
    1239          * @throws MathArithmeticException if the input array contains infinite
    1240          * elements or sums to zero.
    1241          * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}.
    1242          * @since 2.1
    1243          */
    1244         public static double[] normalizeArray(double[] values, double normalizedSum)
    1245             throws MathIllegalArgumentException, MathArithmeticException {
    1246             if (Double.isInfinite(normalizedSum)) {
    1247                 throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_INFINITE);
    1248             }
    1249             if (Double.isNaN(normalizedSum)) {
    1250                 throw new MathIllegalArgumentException(LocalizedFormats.NORMALIZE_NAN);
    1251             }
    1252             double sum = 0d;
    1253             final int len = values.length;
    1254             double[] out = new double[len];
    1255             for (int i = 0; i < len; i++) {
    1256                 if (Double.isInfinite(values[i])) {
    1257                     throw new MathIllegalArgumentException(LocalizedFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
    1258                 }
    1259                 if (!Double.isNaN(values[i])) {
    1260                     sum += values[i];
    1261                 }
    1262             }
    1263             if (sum == 0) {
    1264                 throw new MathArithmeticException(LocalizedFormats.ARRAY_SUMS_TO_ZERO);
    1265             }
    1266             for (int i = 0; i < len; i++) {
    1267                 if (Double.isNaN(values[i])) {
    1268                     out[i] = Double.NaN;
    1269                 } else {
    1270                     out[i] = values[i] * normalizedSum / sum;
    1271                 }
    1272             }
    1273             return out;
    1274         }
    1275    }