MathArrays.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.util;

  25. import java.lang.reflect.Array;
  26. import java.util.ArrayList;
  27. import java.util.Arrays;
  28. import java.util.Comparator;
  29. import java.util.Iterator;
  30. import java.util.List;
  31. import java.util.NavigableSet;
  32. import java.util.TreeSet;

  34. import org.hipparchus.Field;
  35. import org.hipparchus.FieldElement;
  36. import org.hipparchus.CalculusFieldElement;
  37. import org.hipparchus.exception.LocalizedCoreFormats;
  38. import org.hipparchus.exception.MathIllegalArgumentException;
  39. import org.hipparchus.exception.MathRuntimeException;
  40. import org.hipparchus.exception.NullArgumentException;
  41. import org.hipparchus.random.RandomGenerator;
  42. import org.hipparchus.random.Well19937c;

  44. /**
  45.  * Arrays utilities.
  46.  */
  47. public class MathArrays {

  49.     /**
  50.      * Private constructor.
  51.      */
  52.     private MathArrays() {}

  54.     /**
  55.      * Real-valued function that operates on an array or a part of it.
  56.      */
  57.     public interface Function {

  59.         /** Operates on an entire array.
  60.          * @param array Array to operate on.
  61.          * @return the result of the operation.
  62.          */
  63.         double evaluate(double[] array);

  65.         /** Operates on a sub-array.
  66.          * @param array Array to operate on.
  67.          * @param startIndex Index of the first element to take into account.
  68.          * @param numElements Number of elements to take into account.
  69.          * @return the result of the operation.
  70.          */
  71.         double evaluate(double[] array,
  72.                         int startIndex,
  73.                         int numElements);

  75.     }

  77.     /**
  78.      * Create a copy of an array scaled by a value.
  79.      *
  80.      * @param arr Array to scale.
  81.      * @param val Scalar.
  82.      * @return scaled copy of array with each entry multiplied by val.
  83.      */
  84.     public static double[] scale(double val, final double[] arr) {
  85.         double[] newArr = new double[arr.length];
  86.         for (int i = 0; i < arr.length; i++) {
  87.             newArr[i] = arr[i] * val;
  88.         }
  89.         return newArr;
  90.     }

  92.     /**
  93.      * Multiply each element of an array by a value.
  94.      * <p>
  95.      * The array is modified in place (no copy is created).
  96.      *
  97.      * @param arr Array to scale
  98.      * @param val Scalar
  99.      */
  100.     public static void scaleInPlace(double val, final double[] arr) {
  101.         for (int i = 0; i < arr.length; i++) {
  102.             arr[i] *= val;
  103.         }
  104.     }

  106.     /**
  107.      * Creates an array whose contents will be the element-by-element
  108.      * addition of the arguments.
  109.      *
  110.      * @param a First term of the addition.
  111.      * @param b Second term of the addition.
  112.      * @return a new array {@code r} where {@code r[i] = a[i] + b[i]}.
  113.      * @throws MathIllegalArgumentException if the array lengths differ.
  114.      */
  115.     public static double[] ebeAdd(double[] a, double[] b)
  116.         throws MathIllegalArgumentException {
  117.         checkEqualLength(a, b);

  119.         final double[] result = a.clone();
  120.         for (int i = 0; i < a.length; i++) {
  121.             result[i] += b[i];
  122.         }
  123.         return result;
  124.     }
  125.     /**
  126.      * Creates an array whose contents will be the element-by-element
  127.      * subtraction of the second argument from the first.
  128.      *
  129.      * @param a First term.
  130.      * @param b Element to be subtracted.
  131.      * @return a new array {@code r} where {@code r[i] = a[i] - b[i]}.
  132.      * @throws MathIllegalArgumentException if the array lengths differ.
  133.      */
  134.     public static double[] ebeSubtract(double[] a, double[] b)
  135.         throws MathIllegalArgumentException {
  136.         checkEqualLength(a, b);

  138.         final double[] result = a.clone();
  139.         for (int i = 0; i < a.length; i++) {
  140.             result[i] -= b[i];
  141.         }
  142.         return result;
  143.     }
  144.     /**
  145.      * Creates an array whose contents will be the element-by-element
  146.      * multiplication of the arguments.
  147.      *
  148.      * @param a First factor of the multiplication.
  149.      * @param b Second factor of the multiplication.
  150.      * @return a new array {@code r} where {@code r[i] = a[i] * b[i]}.
  151.      * @throws MathIllegalArgumentException if the array lengths differ.
  152.      */
  153.     public static double[] ebeMultiply(double[] a, double[] b)
  154.         throws MathIllegalArgumentException {
  155.         checkEqualLength(a, b);

  157.         final double[] result = a.clone();
  158.         for (int i = 0; i < a.length; i++) {
  159.             result[i] *= b[i];
  160.         }
  161.         return result;
  162.     }
  163.     /**
  164.      * Creates an array whose contents will be the element-by-element
  165.      * division of the first argument by the second.
  166.      *
  167.      * @param a Numerator of the division.
  168.      * @param b Denominator of the division.
  169.      * @return a new array {@code r} where {@code r[i] = a[i] / b[i]}.
  170.      * @throws MathIllegalArgumentException if the array lengths differ.
  171.      */
  172.     public static double[] ebeDivide(double[] a, double[] b)
  173.         throws MathIllegalArgumentException {
  174.         checkEqualLength(a, b);

  176.         final double[] result = a.clone();
  177.         for (int i = 0; i < a.length; i++) {
  178.             result[i] /= b[i];
  179.         }
  180.         return result;
  181.     }

  183.     /**
  184.      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
  185.      *
  186.      * @param p1 the first point
  187.      * @param p2 the second point
  188.      * @return the L<sub>1</sub> distance between the two points
  189.      * @throws MathIllegalArgumentException if the array lengths differ.
  190.      */
  191.     public static double distance1(double[] p1, double[] p2)
  192.         throws MathIllegalArgumentException {
  193.         checkEqualLength(p1, p2);
  194.         double sum = 0;
  195.         for (int i = 0; i < p1.length; i++) {
  196.             sum += FastMath.abs(p1[i] - p2[i]);
  197.         }
  198.         return sum;
  199.     }

  201.     /**
  202.      * Calculates the L<sub>1</sub> (sum of abs) distance between two points.
  203.      *
  204.      * @param p1 the first point
  205.      * @param p2 the second point
  206.      * @return the L<sub>1</sub> distance between the two points
  207.      * @throws MathIllegalArgumentException if the array lengths differ.
  208.      */
  209.     public static int distance1(int[] p1, int[] p2)
  210.         throws MathIllegalArgumentException {
  211.         checkEqualLength(p1, p2);
  212.         int sum = 0;
  213.         for (int i = 0; i < p1.length; i++) {
  214.             sum += FastMath.abs(p1[i] - p2[i]);
  215.         }
  216.         return sum;
  217.     }

  219.     /**
  220.      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
  221.      *
  222.      * @param p1 the first point
  223.      * @param p2 the second point
  224.      * @return the L<sub>2</sub> distance between the two points
  225.      * @throws MathIllegalArgumentException if the array lengths differ.
  226.      */
  227.     public static double distance(double[] p1, double[] p2)
  228.         throws MathIllegalArgumentException {
  229.         checkEqualLength(p1, p2);
  230.         double sum = 0;
  231.         for (int i = 0; i < p1.length; i++) {
  232.             final double dp = p1[i] - p2[i];
  233.             sum += dp * dp;
  234.         }
  235.         return FastMath.sqrt(sum);
  236.     }

  238.     /**
  239.      * Calculates the cosine of the angle between two vectors.
  240.      *
  241.      * @param v1 Cartesian coordinates of the first vector.
  242.      * @param v2 Cartesian coordinates of the second vector.
  243.      * @return the cosine of the angle between the vectors.
  244.      */
  245.     public static double cosAngle(double[] v1, double[] v2) {
  246.         return linearCombination(v1, v2) / (safeNorm(v1) * safeNorm(v2));
  247.     }

  249.     /**
  250.      * Calculates the L<sub>2</sub> (Euclidean) distance between two points.
  251.      *
  252.      * @param p1 the first point
  253.      * @param p2 the second point
  254.      * @return the L<sub>2</sub> distance between the two points
  255.      * @throws MathIllegalArgumentException if the array lengths differ.
  256.      */
  257.     public static double distance(int[] p1, int[] p2)
  258.         throws MathIllegalArgumentException {
  259.         checkEqualLength(p1, p2);
  260.         double sum = 0;
  261.         for (int i = 0; i < p1.length; i++) {
  262.             final double dp = p1[i] - p2[i];
  263.             sum += dp * dp;
  264.         }
  265.         return FastMath.sqrt(sum);
  266.     }

  268.     /**
  269.      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
  270.      *
  271.      * @param p1 the first point
  272.      * @param p2 the second point
  273.      * @return the L<sub>&infin;</sub> distance between the two points
  274.      * @throws MathIllegalArgumentException if the array lengths differ.
  275.      */
  276.     public static double distanceInf(double[] p1, double[] p2)
  277.         throws MathIllegalArgumentException {
  278.         checkEqualLength(p1, p2);
  279.         double max = 0;
  280.         for (int i = 0; i < p1.length; i++) {
  281.             max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
  282.         }
  283.         return max;
  284.     }

  286.     /**
  287.      * Calculates the L<sub>&infin;</sub> (max of abs) distance between two points.
  288.      *
  289.      * @param p1 the first point
  290.      * @param p2 the second point
  291.      * @return the L<sub>&infin;</sub> distance between the two points
  292.      * @throws MathIllegalArgumentException if the array lengths differ.
  293.      */
  294.     public static int distanceInf(int[] p1, int[] p2)
  295.         throws MathIllegalArgumentException {
  296.         checkEqualLength(p1, p2);
  297.         int max = 0;
  298.         for (int i = 0; i < p1.length; i++) {
  299.             max = FastMath.max(max, FastMath.abs(p1[i] - p2[i]));
  300.         }
  301.         return max;
  302.     }

  304.     /**
  305.      * Specification of ordering direction.
  306.      */
  307.     public enum OrderDirection {
  308.         /** Constant for increasing direction. */
  309.         INCREASING,
  310.         /** Constant for decreasing direction. */
  311.         DECREASING
  312.     }

  314.     /**
  315.      * Check that an array is monotonically increasing or decreasing.
  316.      *
  317.      * @param <T> the type of the elements in the specified array
  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 <T extends Comparable<? super T>> boolean isMonotonic(T[] val,
  324.                                                                         OrderDirection dir,
  325.                                                                         boolean strict) {
  326.         T previous = val[0];
  327.         final int max = val.length;
  328.         for (int i = 1; i < max; i++) {
  329.             final int comp;
  330.             switch (dir) {
  331.             case INCREASING:
  332.                 comp = previous.compareTo(val[i]);
  333.                 if (strict) {
  334.                     if (comp >= 0) {
  335.                         return false;
  336.                     }
  337.                 } else {
  338.                     if (comp > 0) {
  339.                         return false;
  340.                     }
  341.                 }
  342.                 break;
  343.             case DECREASING:
  344.                 comp = val[i].compareTo(previous);
  345.                 if (strict) {
  346.                     if (comp >= 0) {
  347.                         return false;
  348.                     }
  349.                 } else {
  350.                     if (comp > 0) {
  351.                        return false;
  352.                     }
  353.                 }
  354.                 break;
  355.             default:
  356.                 // Should never happen.
  357.                 throw MathRuntimeException.createInternalError();
  358.             }

  360.             previous = val[i];
  361.         }
  362.         return true;
  363.     }

  365.     /**
  366.      * Check that an array is monotonically increasing or decreasing.
  367.      *
  368.      * @param val Values.
  369.      * @param dir Ordering direction.
  370.      * @param strict Whether the order should be strict.
  371.      * @return {@code true} if sorted, {@code false} otherwise.
  372.      */
  373.     public static boolean isMonotonic(double[] val, OrderDirection dir, boolean strict) {
  374.         return checkOrder(val, dir, strict, false);
  375.     }

  377.     /**
  378.      * Check that both arrays have the same length.
  379.      *
  380.      * @param a Array.
  381.      * @param b Array.
  382.      * @param abort Whether to throw an exception if the check fails.
  383.      * @return {@code true} if the arrays have the same length.
  384.      * @throws MathIllegalArgumentException if the lengths differ and
  385.      * {@code abort} is {@code true}.
  386.      */
  387.     public static boolean checkEqualLength(double[] a,
  388.                                            double[] b,
  389.                                            boolean abort) {
  390.         if (a.length == b.length) {
  391.             return true;
  392.         } else {
  393.             if (abort) {
  394.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  395.                                                        a.length, b.length);
  396.             }
  397.             return false;
  398.         }
  399.     }

  401.     /**
  402.      * Check that both arrays have the same length.
  403.      *
  404.      * @param a Array.
  405.      * @param b Array.
  406.      * @throws MathIllegalArgumentException if the lengths differ.
  407.      */
  408.     public static void checkEqualLength(double[] a,
  409.                                         double[] b) {
  410.         checkEqualLength(a, b, true);
  411.     }

  413.     /**
  414.      * Check that both arrays have the same length.
  415.      *
  416.      * @param a Array.
  417.      * @param b Array.
  418.      * @param abort Whether to throw an exception if the check fails.
  419.      * @return {@code true} if the arrays have the same length.
  420.      * @throws MathIllegalArgumentException if the lengths differ and
  421.      * {@code abort} is {@code true}.
  422.      * @param <T> the type of the field elements
  423.      * @since 1.5
  424.      */
  425.     public static <T extends CalculusFieldElement<T>> boolean checkEqualLength(final T[] a,
  426.                                                                            final T[] b,
  427.                                            boolean abort) {
  428.         if (a.length == b.length) {
  429.             return true;
  430.         } else {
  431.             if (abort) {
  432.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  433.                                                        a.length, b.length);
  434.             }
  435.             return false;
  436.         }
  437.     }

  439.     /**
  440.      * Check that both arrays have the same length.
  441.      *
  442.      * @param a Array.
  443.      * @param b Array.
  444.      * @throws MathIllegalArgumentException if the lengths differ.
  445.      * @param <T> the type of the field elements
  446.      * @since 1.5
  447.      */
  448.     public static <T extends CalculusFieldElement<T>> void checkEqualLength(final T[] a, final T[] b) {
  449.         checkEqualLength(a, b, true);
  450.     }

  452.     /**
  453.      * Check that both arrays have the same length.
  454.      *
  455.      * @param a Array.
  456.      * @param b Array.
  457.      * @param abort Whether to throw an exception if the check fails.
  458.      * @return {@code true} if the arrays have the same length.
  459.      * @throws MathIllegalArgumentException if the lengths differ and
  460.      * {@code abort} is {@code true}.
  461.      */
  462.     public static boolean checkEqualLength(int[] a,
  463.                                            int[] b,
  464.                                            boolean abort) {
  465.         if (a.length == b.length) {
  466.             return true;
  467.         } else {
  468.             if (abort) {
  469.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  470.                                                        a.length, b.length);
  471.             }
  472.             return false;
  473.         }
  474.     }

  476.     /**
  477.      * Check that both arrays have the same length.
  478.      *
  479.      * @param a Array.
  480.      * @param b Array.
  481.      * @throws MathIllegalArgumentException if the lengths differ.
  482.      */
  483.     public static void checkEqualLength(int[] a,
  484.                                         int[] b) {
  485.         checkEqualLength(a, b, true);
  486.     }

  488.     /**
  489.      * Check that the given array is sorted.
  490.      *
  491.      * @param val Values.
  492.      * @param dir Ordering direction.
  493.      * @param strict Whether the order should be strict.
  494.      * @param abort Whether to throw an exception if the check fails.
  495.      * @return {@code true} if the array is sorted.
  496.      * @throws MathIllegalArgumentException if the array is not sorted
  497.      * and {@code abort} is {@code true}.
  498.      */
  499.     public static boolean checkOrder(double[] val, OrderDirection dir,
  500.                                      boolean strict, boolean abort)
  501.         throws MathIllegalArgumentException {
  502.         double previous = val[0];
  503.         final int max = val.length;

  505.         int index;
  506.         ITEM:
  507.         for (index = 1; index < max; index++) {
  508.             switch (dir) {
  509.             case INCREASING:
  510.                 if (strict) {
  511.                     if (val[index] <= previous) {
  512.                         break ITEM;
  513.                     }
  514.                 } else {
  515.                     if (val[index] < previous) {
  516.                         break ITEM;
  517.                     }
  518.                 }
  519.                 break;
  520.             case DECREASING:
  521.                 if (strict) {
  522.                     if (val[index] >= previous) {
  523.                         break ITEM;
  524.                     }
  525.                 } else {
  526.                     if (val[index] > previous) {
  527.                         break ITEM;
  528.                     }
  529.                 }
  530.                 break;
  531.             default:
  532.                 // Should never happen.
  533.                 throw MathRuntimeException.createInternalError();
  534.             }

  536.             previous = val[index];
  537.         }

  539.         if (index == max) {
  540.             // Loop completed.
  541.             return true;
  542.         }

  544.         // Loop early exit means wrong ordering.
  545.         if (abort) {
  546.             throw new MathIllegalArgumentException(dir == MathArrays.OrderDirection.INCREASING ?
  547.                                                     (strict ?
  548.                                                      LocalizedCoreFormats.NOT_STRICTLY_INCREASING_SEQUENCE :
  549.                                                      LocalizedCoreFormats.NOT_INCREASING_SEQUENCE) :
  550.                                                     (strict ?
  551.                                                      LocalizedCoreFormats.NOT_STRICTLY_DECREASING_SEQUENCE :
  552.                                                      LocalizedCoreFormats.NOT_DECREASING_SEQUENCE),
  553.                                                     val[index], previous, index, index - 1);
  554.         } else {
  555.             return false;
  556.         }
  557.     }

  559.     /**
  560.      * Check that the given array is sorted.
  561.      *
  562.      * @param val Values.
  563.      * @param dir Ordering direction.
  564.      * @param strict Whether the order should be strict.
  565.      * @throws MathIllegalArgumentException if the array is not sorted.
  566.      */
  567.     public static void checkOrder(double[] val, OrderDirection dir,
  568.                                   boolean strict) throws MathIllegalArgumentException {
  569.         checkOrder(val, dir, strict, true);
  570.     }

  572.     /**
  573.      * Check that the given array is sorted in strictly increasing order.
  574.      *
  575.      * @param val Values.
  576.      * @throws MathIllegalArgumentException if the array is not sorted.
  577.      */
  578.     public static void checkOrder(double[] val) throws MathIllegalArgumentException {
  579.         checkOrder(val, OrderDirection.INCREASING, true);
  580.     }

  582.     /**
  583.      * Check that the given array is sorted.
  584.      *
  585.      * @param val Values.
  586.      * @param dir Ordering direction.
  587.      * @param strict Whether the order should be strict.
  588.      * @param abort Whether to throw an exception if the check fails.
  589.      * @return {@code true} if the array is sorted.
  590.      * @throws MathIllegalArgumentException if the array is not sorted
  591.      * and {@code abort} is {@code true}.
  592.      * @param <T> the type of the field elements
  593.      * @since 1.5
  594.      */
  595.     public static <T extends CalculusFieldElement<T>>boolean checkOrder(T[] val, OrderDirection dir,
  596.                                                                     boolean strict, boolean abort)
  597.         throws MathIllegalArgumentException {
  598.         double previous = val[0].getReal();
  599.         final int max = val.length;

  601.         int index;
  602.         ITEM:
  603.         for (index = 1; index < max; index++) {
  604.             switch (dir) {
  605.             case INCREASING:
  606.                 if (strict) {
  607.                     if (val[index].getReal() <= previous) {
  608.                         break ITEM;
  609.                     }
  610.                 } else {
  611.                     if (val[index].getReal() < previous) {
  612.                         break ITEM;
  613.                     }
  614.                 }
  615.                 break;
  616.             case DECREASING:
  617.                 if (strict) {
  618.                     if (val[index].getReal() >= previous) {
  619.                         break ITEM;
  620.                     }
  621.                 } else {
  622.                     if (val[index].getReal() > previous) {
  623.                         break ITEM;
  624.                     }
  625.                 }
  626.                 break;
  627.             default:
  628.                 // Should never happen.
  629.                 throw MathRuntimeException.createInternalError();
  630.             }

  632.             previous = val[index].getReal();
  633.         }

  635.         if (index == max) {
  636.             // Loop completed.
  637.             return true;
  638.         }

  640.         // Loop early exit means wrong ordering.
  641.         if (abort) {
  642.             throw new MathIllegalArgumentException(dir == MathArrays.OrderDirection.INCREASING ?
  643.                                                     (strict ?
  644.                                                      LocalizedCoreFormats.NOT_STRICTLY_INCREASING_SEQUENCE :
  645.                                                      LocalizedCoreFormats.NOT_INCREASING_SEQUENCE) :
  646.                                                     (strict ?
  647.                                                      LocalizedCoreFormats.NOT_STRICTLY_DECREASING_SEQUENCE :
  648.                                                      LocalizedCoreFormats.NOT_DECREASING_SEQUENCE),
  649.                                                     val[index], previous, index, index - 1);
  650.         } else {
  651.             return false;
  652.         }
  653.     }

  655.     /**
  656.      * Check that the given array is sorted.
  657.      *
  658.      * @param val Values.
  659.      * @param dir Ordering direction.
  660.      * @param strict Whether the order should be strict.
  661.      * @throws MathIllegalArgumentException if the array is not sorted.
  662.      * @param <T> the type of the field elements
  663.      * @since 1.5
  664.      */
  665.     public static <T extends CalculusFieldElement<T>> void checkOrder(T[] val, OrderDirection dir,
  666.                                                                   boolean strict) throws MathIllegalArgumentException {
  667.         checkOrder(val, dir, strict, true);
  668.     }

  670.     /**
  671.      * Check that the given array is sorted in strictly increasing order.
  672.      *
  673.      * @param val Values.
  674.      * @throws MathIllegalArgumentException if the array is not sorted.
  675.      * @param <T> the type of the field elements
  676.      * @since 1.5
  677.      */
  678.     public static <T extends CalculusFieldElement<T>> void checkOrder(T[] val) throws MathIllegalArgumentException {
  679.         checkOrder(val, OrderDirection.INCREASING, true);
  680.     }

  682.     /**
  683.      * Throws MathIllegalArgumentException if the input array is not rectangular.
  684.      *
  685.      * @param in array to be tested
  686.      * @throws NullArgumentException if input array is null
  687.      * @throws MathIllegalArgumentException if input array is not rectangular
  688.      */
  689.     public static void checkRectangular(final long[][] in)
  690.         throws MathIllegalArgumentException, NullArgumentException {
  691.         MathUtils.checkNotNull(in);
  692.         for (int i = 1; i < in.length; i++) {
  693.             if (in[i].length != in[0].length) {
  694.                 throw new MathIllegalArgumentException(
  695.                         LocalizedCoreFormats.DIFFERENT_ROWS_LENGTHS,
  696.                         in[i].length, in[0].length);
  697.             }
  698.         }
  699.     }

  701.     /**
  702.      * Check that all entries of the input array are strictly positive.
  703.      *
  704.      * @param in Array to be tested
  705.      * @throws MathIllegalArgumentException if any entries of the array are not
  706.      * strictly positive.
  707.      */
  708.     public static void checkPositive(final double[] in)
  709.         throws MathIllegalArgumentException {
  710.         for (double v : in) {
  711.             if (v <= 0) {
  712.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
  713.                         v, 0);
  714.             }
  715.         }
  716.     }

  718.     /**
  719.      * Check that no entry of the input array is {@code NaN}.
  720.      *
  721.      * @param in Array to be tested.
  722.      * @throws MathIllegalArgumentException if an entry is {@code NaN}.
  723.      */
  724.     public static void checkNotNaN(final double[] in)
  725.         throws MathIllegalArgumentException {
  726.         for(int i = 0; i < in.length; i++) {
  727.             if (Double.isNaN(in[i])) {
  728.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NAN_ELEMENT_AT_INDEX, i);
  729.             }
  730.         }
  731.     }

  733.     /**
  734.      * Check that all entries of the input array are &gt;= 0.
  735.      *
  736.      * @param in Array to be tested
  737.      * @throws MathIllegalArgumentException if any array entries are less than 0.
  738.      */
  739.     public static void checkNonNegative(final long[] in)
  740.         throws MathIllegalArgumentException {
  741.         for (long l : in) {
  742.             if (l < 0) {
  743.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL, l, 0);
  744.             }
  745.         }
  746.     }

  748.     /**
  749.      * Check all entries of the input array are &gt;= 0.
  750.      *
  751.      * @param in Array to be tested
  752.      * @throws MathIllegalArgumentException if any array entries are less than 0.
  753.      */
  754.     public static void checkNonNegative(final long[][] in)
  755.         throws MathIllegalArgumentException {
  756.         for (long[] longs : in) {
  757.             for (int j = 0; j < longs.length; j++) {
  758.                 if (longs[j] < 0) {
  759.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL, longs[j], 0);
  760.                 }
  761.             }
  762.         }
  763.     }

  765.     /**
  766.      * Returns the Cartesian norm (2-norm), handling both overflow and underflow.
  767.      * Translation of the minpack enorm subroutine.
  768.      *
  769.      * <p>
  770.      * The redistribution policy for MINPACK is available
  771.      * <a href="http://www.netlib.org/minpack/disclaimer">here</a>, for
  772.      * convenience, it is reproduced below.
  773.      * </p>
  774.      *
  775.      * <blockquote>
  776.      * <p>
  777.      *    Minpack Copyright Notice (1999) University of Chicago.
  778.      *    All rights reserved
  779.      * </p>
  780.      * <p>
  781.      * Redistribution and use in source and binary forms, with or without
  782.      * modification, are permitted provided that the following conditions
  783.      * are met:
  784.      * </p>
  785.      * <ol>
  786.      *  <li>Redistributions of source code must retain the above copyright
  787.      *      notice, this list of conditions and the following disclaimer.</li>
  788.      * <li>Redistributions in binary form must reproduce the above
  789.      *     copyright notice, this list of conditions and the following
  790.      *     disclaimer in the documentation and/or other materials provided
  791.      *     with the distribution.</li>
  792.      * <li>The end-user documentation included with the redistribution, if any,
  793.      *     must include the following acknowledgment:
  794.      *     {@code This product includes software developed by the University of
  795.      *           Chicago, as Operator of Argonne National Laboratory.}
  796.      *     Alternately, this acknowledgment may appear in the software itself,
  797.      *     if and wherever such third-party acknowledgments normally appear.</li>
  798.      * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
  799.      *     WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
  800.      *     UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
  801.      *     THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
  802.      *     IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
  803.      *     OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
  804.      *     OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
  805.      *     OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
  806.      *     USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
  807.      *     THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
  808.      *     DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
  809.      *     UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
  810.      *     BE CORRECTED.</strong></li>
  811.      * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
  812.      *     HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
  813.      *     ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
  814.      *     INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
  815.      *     ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
  816.      *     PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
  817.      *     SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
  818.      *     (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
  819.      *     EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
  820.      *     POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li>
  821.      * </ol>
  822.      * </blockquote>
  823.      *
  824.      * @param v Vector of doubles.
  825.      * @return the 2-norm of the vector.
  826.      */
  827.     public static double safeNorm(double[] v) {
  828.         double rdwarf = 3.834e-20;
  829.         double rgiant = 1.304e+19;
  830.         double s1 = 0;
  831.         double s2 = 0;
  832.         double s3 = 0;
  833.         double x1max = 0;
  834.         double x3max = 0;
  835.         double floatn = v.length;
  836.         double agiant = rgiant / floatn;
  837.         for (double value : v) {
  838.             double xabs = FastMath.abs(value);
  839.             if (xabs < rdwarf || xabs > agiant) {
  840.                 if (xabs > rdwarf) {
  841.                     if (xabs > x1max) {
  842.                         double r = x1max / xabs;
  843.                         s1 = 1 + s1 * r * r;
  844.                         x1max = xabs;
  845.                     } else {
  846.                         double r = xabs / x1max;
  847.                         s1 += r * r;
  848.                     }
  849.                 } else {
  850.                     if (xabs > x3max) {
  851.                         double r = x3max / xabs;
  852.                         s3 = 1 + s3 * r * r;
  853.                         x3max = xabs;
  854.                     } else {
  855.                         if (xabs != 0) {
  856.                             double r = xabs / x3max;
  857.                             s3 += r * r;
  858.                         }
  859.                     }
  860.                 }
  861.             } else {
  862.                 s2 += xabs * xabs;
  863.             }
  864.         }
  865.         double norm;
  866.         if (s1 != 0) {
  867.             norm = x1max * Math.sqrt(s1 + (s2 / x1max) / x1max);
  868.         } else {
  869.             if (s2 == 0) {
  870.                 norm = x3max * Math.sqrt(s3);
  871.             } else {
  872.                 if (s2 >= x3max) {
  873.                     norm = Math.sqrt(s2 * (1 + (x3max / s2) * (x3max * s3)));
  874.                 } else {
  875.                     norm = Math.sqrt(x3max * ((s2 / x3max) + (x3max * s3)));
  876.                 }
  877.             }
  878.         }
  879.         return norm;
  880.     }

  882.     /**
  883.      * Sort an array in ascending order in place and perform the same reordering
  884.      * of entries on other arrays. For example, if
  885.      * {@code x = [3, 1, 2], y = [1, 2, 3]} and {@code z = [0, 5, 7]}, then
  886.      * {@code sortInPlace(x, y, z)} will update {@code x} to {@code [1, 2, 3]},
  887.      * {@code y} to {@code [2, 3, 1]} and {@code z} to {@code [5, 7, 0]}.
  888.      *
  889.      * @param x Array to be sorted and used as a pattern for permutation
  890.      * of the other arrays.
  891.      * @param yList Set of arrays whose permutations of entries will follow
  892.      * those performed on {@code x}.
  893.      * @throws MathIllegalArgumentException if any {@code y} is not the same
  894.      * size as {@code x}.
  895.      * @throws NullArgumentException if {@code x} or any {@code y} is null.
  896.      */
  897.     public static void sortInPlace(double[] x, double[]... yList)
  898.         throws MathIllegalArgumentException, NullArgumentException {
  899.         sortInPlace(x, OrderDirection.INCREASING, yList);
  900.     }

  902.     /**
  903.      * Helper data structure holding a (double, integer) pair.
  904.      */
  905.     private static class PairDoubleInteger {
  906.         /** Key */
  907.         private final double key;
  908.         /** Value */
  909.         private final int value;

  911.         /**
  912.          * @param key Key.
  913.          * @param value Value.
  914.          */
  915.         PairDoubleInteger(double key, int value) {
  916.             this.key = key;
  917.             this.value = value;
  918.         }

  920.         /** @return the key. */
  921.         public double getKey() {
  922.             return key;
  923.         }

  925.         /** @return the value. */
  926.         public int getValue() {
  927.             return value;
  928.         }
  929.     }

  931.     /**
  932.      * Sort an array in place and perform the same reordering of entries on
  933.      * other arrays.  This method works the same as the other
  934.      * {@link #sortInPlace(double[], double[][]) sortInPlace} method, but
  935.      * allows the order of the sort to be provided in the {@code dir}
  936.      * parameter.
  937.      *
  938.      * @param x Array to be sorted and used as a pattern for permutation
  939.      * of the other arrays.
  940.      * @param dir Order direction.
  941.      * @param yList Set of arrays whose permutations of entries will follow
  942.      * those performed on {@code x}.
  943.      * @throws MathIllegalArgumentException if any {@code y} is not the same
  944.      * size as {@code x}.
  945.      * @throws NullArgumentException if {@code x} or any {@code y} is null
  946.      */
  947.     public static void sortInPlace(double[] x,
  948.                                    final OrderDirection dir,
  949.                                    double[]... yList)
  950.         throws MathIllegalArgumentException, NullArgumentException {

  952.         // Consistency checks.
  953.         if (x == null) {
  954.             throw new NullArgumentException();
  955.         }

  957.         final int len = x.length;

  959.         for (final double[] y : yList) {
  960.             if (y == null) {
  961.                 throw new NullArgumentException();
  962.             }
  963.             if (y.length != len) {
  964.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  965.                         y.length, len);
  966.             }
  967.         }

  969.         // Associate each abscissa "x[i]" with its index "i".
  970.         final List<PairDoubleInteger> list = new ArrayList<>(len);
  971.         for (int i = 0; i < len; i++) {
  972.             list.add(new PairDoubleInteger(x[i], i));
  973.         }

  975.         // Create comparators for increasing and decreasing orders.
  976.         final Comparator<PairDoubleInteger> comp =
  977.             dir == MathArrays.OrderDirection.INCREASING ?
  978.             new Comparator<PairDoubleInteger>() {
  979.                 /** {@inheritDoc} */
  980.                 @Override
  981.                 public int compare(PairDoubleInteger o1,
  982.                                    PairDoubleInteger o2) {
  983.                     return Double.compare(o1.getKey(), o2.getKey());
  984.                 }
  985.             } :
  986.             new Comparator<PairDoubleInteger>() {
  987.                 /** {@inheritDoc} */
  988.                 @Override
  989.                 public int compare(PairDoubleInteger o1,
  990.                                    PairDoubleInteger o2) {
  991.                     return Double.compare(o2.getKey(), o1.getKey());
  992.                 }
  993.             };

  995.         // Sort.
  996.         list.sort(comp);

  998.         // Modify the original array so that its elements are in
  999.         // the prescribed order.
  1000.         // Retrieve indices of original locations.
  1001.         final int[] indices = new int[len];
  1002.         for (int i = 0; i < len; i++) {
  1003.             final PairDoubleInteger e = list.get(i);
  1004.             x[i] = e.getKey();
  1005.             indices[i] = e.getValue();
  1006.         }

  1008.         // In each of the associated arrays, move the
  1009.         // elements to their new location.
  1010.         for (final double[] yInPlace : yList) {
  1011.             // Input array will be modified in place.
  1012.             final double[] yOrig = yInPlace.clone();

  1014.             for (int i = 0; i < len; i++) {
  1015.                 yInPlace[i] = yOrig[indices[i]];
  1016.             }
  1017.         }
  1018.     }

  1020.     /**
  1021.      * Compute a linear combination accurately.
  1022.      * This method computes the sum of the products
  1023.      * <code>a<sub>i</sub> b<sub>i</sub></code> to high accuracy.
  1024.      * It does so by using specific multiplication and addition algorithms to
  1025.      * preserve accuracy and reduce cancellation effects.
  1026.      * <br>
  1027.      * It is based on the 2005 paper
  1028.      * <a href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
  1029.      * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
  1030.      * and Shin'ichi Oishi published in SIAM J. Sci. Comput.
  1031.      *
  1032.      * @param a Factors.
  1033.      * @param b Factors.
  1034.      * @return <code>&Sigma;<sub>i</sub> a<sub>i</sub> b<sub>i</sub></code>.
  1035.      * @throws MathIllegalArgumentException if arrays dimensions don't match
  1036.      */
  1037.     public static double linearCombination(final double[] a, final double[] b)
  1038.         throws MathIllegalArgumentException {
  1039.         checkEqualLength(a, b);
  1040.         final int len = a.length;

  1042.         if (len == 1) {
  1043.             // Revert to scalar multiplication.
  1044.             return a[0] * b[0];
  1045.         }

  1047.         final double[] prodHigh = new double[len];
  1048.         double prodLowSum = 0;

  1050.         for (int i = 0; i < len; i++) {
  1051.             final double ai    = a[i];
  1052.             final double aHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(ai) & ((-1L) << 27));
  1053.             final double aLow  = ai - aHigh;

  1055.             final double bi    = b[i];
  1056.             final double bHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(bi) & ((-1L) << 27));
  1057.             final double bLow  = bi - bHigh;
  1058.             prodHigh[i] = ai * bi;
  1059.             final double prodLow = aLow * bLow - (((prodHigh[i] -
  1060.                                                     aHigh * bHigh) -
  1061.                                                    aLow * bHigh) -
  1062.                                                   aHigh * bLow);
  1063.             prodLowSum += prodLow;
  1064.         }


  1067.         final double prodHighCur = prodHigh[0];
  1068.         double prodHighNext = prodHigh[1];
  1069.         double sHighPrev = prodHighCur + prodHighNext;
  1070.         double sPrime = sHighPrev - prodHighNext;
  1071.         double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime);

  1073.         final int lenMinusOne = len - 1;
  1074.         for (int i = 1; i < lenMinusOne; i++) {
  1075.             prodHighNext = prodHigh[i + 1];
  1076.             final double sHighCur = sHighPrev + prodHighNext;
  1077.             sPrime = sHighCur - prodHighNext;
  1078.             sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime);
  1079.             sHighPrev = sHighCur;
  1080.         }

  1082.         double result = sHighPrev + (prodLowSum + sLowSum);

  1084.         if (Double.isNaN(result) || result == 0.0) {
  1085.             // either we have split infinite numbers or some coefficients were NaNs or signed zeros,
  1086.             // just rely on the naive implementation and let IEEE754 handle this
  1087.             // we do this for zeros too as we want to preserve the sign of zero (see issue #76)
  1088.             result = a[0] * b[0];
  1089.             for (int i = 1; i < len; ++i) {
  1090.                 result += a[i] * b[i];
  1091.             }
  1092.         }

  1094.         return result;
  1095.     }

  1097.     /**
  1098.      * Compute a linear combination accurately.
  1099.      * <p>
  1100.      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
  1101.      * a<sub>2</sub>&times;b<sub>2</sub> to high accuracy. It does
  1102.      * so by using specific multiplication and addition algorithms to
  1103.      * preserve accuracy and reduce cancellation effects. It is based
  1104.      * on the 2005 paper <a
  1105.      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
  1106.      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
  1107.      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
  1108.      * </p>
  1109.      * @param a1 first factor of the first term
  1110.      * @param b1 second factor of the first term
  1111.      * @param a2 first factor of the second term
  1112.      * @param b2 second factor of the second term
  1113.      * @return a<sub>1</sub>&times;b<sub>1</sub> +
  1114.      * a<sub>2</sub>&times;b<sub>2</sub>
  1115.      * @see #linearCombination(double, double, double, double, double, double)
  1116.      * @see #linearCombination(double, double, double, double, double, double, double, double)
  1117.      */
  1118.     public static double linearCombination(final double a1, final double b1,
  1119.                                            final double a2, final double b2) {

  1121.         // the code below is split in many additions/subtractions that may
  1122.         // appear redundant. However, they should NOT be simplified, as they
  1123.         // use IEEE754 floating point arithmetic rounding properties.
  1124.         // The variable naming conventions are that xyzHigh contains the most significant
  1125.         // bits of xyz and xyzLow contains its least significant bits. So theoretically
  1126.         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
  1127.         // be represented in only one double precision number so we preserve two numbers
  1128.         // to hold it as long as we can, combining the high and low order bits together
  1129.         // only at the end, after cancellation may have occurred on high order bits

  1131.         // split a1 and b1 as one 26 bits number and one 27 bits number
  1132.         final double a1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
  1133.         final double a1Low      = a1 - a1High;
  1134.         final double b1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
  1135.         final double b1Low      = b1 - b1High;

  1137.         // accurate multiplication a1 * b1
  1138.         final double prod1High  = a1 * b1;
  1139.         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);

  1141.         // split a2 and b2 as one 26 bits number and one 27 bits number
  1142.         final double a2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
  1143.         final double a2Low      = a2 - a2High;
  1144.         final double b2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
  1145.         final double b2Low      = b2 - b2High;

  1147.         // accurate multiplication a2 * b2
  1148.         final double prod2High  = a2 * b2;
  1149.         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);

  1151.         // accurate addition a1 * b1 + a2 * b2
  1152.         final double s12High    = prod1High + prod2High;
  1153.         final double s12Prime   = s12High - prod2High;
  1154.         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);

  1156.         // final rounding, s12 may have suffered many cancellations, we try
  1157.         // to recover some bits from the extra words we have saved up to now
  1158.         double result = s12High + (prod1Low + prod2Low + s12Low);

  1160.         if (Double.isNaN(result) || result == 0.0) {
  1161.             // either we have split infinite numbers or some coefficients were NaNs or signed zeros,
  1162.             // just rely on the naive implementation and let IEEE754 handle this
  1163.             // we do this for zeros too as we want to preserve the sign of zero (see issue #76)
  1164.             result = a1 * b1 + a2 * b2;
  1165.         }

  1167.         return result;
  1168.     }

  1170.     /**
  1171.      * Compute a linear combination accurately.
  1172.      * <p>
  1173.      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
  1174.      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
  1175.      * to high accuracy. It does so by using specific multiplication and
  1176.      * addition algorithms to preserve accuracy and reduce cancellation effects.
  1177.      * It is based on the 2005 paper <a
  1178.      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
  1179.      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
  1180.      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
  1181.      * </p>
  1182.      * @param a1 first factor of the first term
  1183.      * @param b1 second factor of the first term
  1184.      * @param a2 first factor of the second term
  1185.      * @param b2 second factor of the second term
  1186.      * @param a3 first factor of the third term
  1187.      * @param b3 second factor of the third term
  1188.      * @return a<sub>1</sub>&times;b<sub>1</sub> +
  1189.      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub>
  1190.      * @see #linearCombination(double, double, double, double)
  1191.      * @see #linearCombination(double, double, double, double, double, double, double, double)
  1192.      */
  1193.     public static double linearCombination(final double a1, final double b1,
  1194.                                            final double a2, final double b2,
  1195.                                            final double a3, final double b3) {

  1197.         // the code below is split in many additions/subtractions that may
  1198.         // appear redundant. However, they should NOT be simplified, as they
  1199.         // do use IEEE754 floating point arithmetic rounding properties.
  1200.         // The variables naming conventions are that xyzHigh contains the most significant
  1201.         // bits of xyz and xyzLow contains its least significant bits. So theoretically
  1202.         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
  1203.         // be represented in only one double precision number so we preserve two numbers
  1204.         // to hold it as long as we can, combining the high and low order bits together
  1205.         // only at the end, after cancellation may have occurred on high order bits

  1207.         // split a1 and b1 as one 26 bits number and one 27 bits number
  1208.         final double a1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
  1209.         final double a1Low      = a1 - a1High;
  1210.         final double b1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
  1211.         final double b1Low      = b1 - b1High;

  1213.         // accurate multiplication a1 * b1
  1214.         final double prod1High  = a1 * b1;
  1215.         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);

  1217.         // split a2 and b2 as one 26 bits number and one 27 bits number
  1218.         final double a2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
  1219.         final double a2Low      = a2 - a2High;
  1220.         final double b2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
  1221.         final double b2Low      = b2 - b2High;

  1223.         // accurate multiplication a2 * b2
  1224.         final double prod2High  = a2 * b2;
  1225.         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);

  1227.         // split a3 and b3 as one 26 bits number and one 27 bits number
  1228.         final double a3High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
  1229.         final double a3Low      = a3 - a3High;
  1230.         final double b3High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
  1231.         final double b3Low      = b3 - b3High;

  1233.         // accurate multiplication a3 * b3
  1234.         final double prod3High  = a3 * b3;
  1235.         final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);

  1237.         // accurate addition a1 * b1 + a2 * b2
  1238.         final double s12High    = prod1High + prod2High;
  1239.         final double s12Prime   = s12High - prod2High;
  1240.         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);

  1242.         // accurate addition a1 * b1 + a2 * b2 + a3 * b3
  1243.         final double s123High   = s12High + prod3High;
  1244.         final double s123Prime  = s123High - prod3High;
  1245.         final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);

  1247.         // final rounding, s123 may have suffered many cancellations, we try
  1248.         // to recover some bits from the extra words we have saved up to now
  1249.         double result = s123High + (prod1Low + prod2Low + prod3Low + s12Low + s123Low);

  1251.         if (Double.isNaN(result) || result == 0.0) {
  1252.             // either we have split infinite numbers or some coefficients were NaNs or signed zeros,
  1253.             // just rely on the naive implementation and let IEEE754 handle this
  1254.             // we do this for zeros too as we want to preserve the sign of zero (see issue #76)
  1255.             result = a1 * b1 + a2 * b2 + a3 * b3;
  1256.         }

  1258.         return result;
  1259.     }

  1261.     /**
  1262.      * Compute a linear combination accurately.
  1263.      * <p>
  1264.      * This method computes a<sub>1</sub>&times;b<sub>1</sub> +
  1265.      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
  1266.      * a<sub>4</sub>&times;b<sub>4</sub>
  1267.      * to high accuracy. It does so by using specific multiplication and
  1268.      * addition algorithms to preserve accuracy and reduce cancellation effects.
  1269.      * It is based on the 2005 paper <a
  1270.      * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
  1271.      * Accurate Sum and Dot Product</a> by Takeshi Ogita,
  1272.      * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. Comput.
  1273.      * </p>
  1274.      * @param a1 first factor of the first term
  1275.      * @param b1 second factor of the first term
  1276.      * @param a2 first factor of the second term
  1277.      * @param b2 second factor of the second term
  1278.      * @param a3 first factor of the third term
  1279.      * @param b3 second factor of the third term
  1280.      * @param a4 first factor of the third term
  1281.      * @param b4 second factor of the third term
  1282.      * @return a<sub>1</sub>&times;b<sub>1</sub> +
  1283.      * a<sub>2</sub>&times;b<sub>2</sub> + a<sub>3</sub>&times;b<sub>3</sub> +
  1284.      * a<sub>4</sub>&times;b<sub>4</sub>
  1285.      * @see #linearCombination(double, double, double, double)
  1286.      * @see #linearCombination(double, double, double, double, double, double)
  1287.      */
  1288.     public static double linearCombination(final double a1, final double b1,
  1289.                                            final double a2, final double b2,
  1290.                                            final double a3, final double b3,
  1291.                                            final double a4, final double b4) {

  1293.         // the code below is split in many additions/subtractions that may
  1294.         // appear redundant. However, they should NOT be simplified, as they
  1295.         // do use IEEE754 floating point arithmetic rounding properties.
  1296.         // The variables naming conventions are that xyzHigh contains the most significant
  1297.         // bits of xyz and xyzLow contains its least significant bits. So theoretically
  1298.         // xyz is the sum xyzHigh + xyzLow, but in many cases below, this sum cannot
  1299.         // be represented in only one double precision number so we preserve two numbers
  1300.         // to hold it as long as we can, combining the high and low order bits together
  1301.         // only at the end, after cancellation may have occurred on high order bits

  1303.         // split a1 and b1 as one 26 bits number and one 27 bits number
  1304.         final double a1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a1) & ((-1L) << 27));
  1305.         final double a1Low      = a1 - a1High;
  1306.         final double b1High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b1) & ((-1L) << 27));
  1307.         final double b1Low      = b1 - b1High;

  1309.         // accurate multiplication a1 * b1
  1310.         final double prod1High  = a1 * b1;
  1311.         final double prod1Low   = a1Low * b1Low - (((prod1High - a1High * b1High) - a1Low * b1High) - a1High * b1Low);

  1313.         // split a2 and b2 as one 26 bits number and one 27 bits number
  1314.         final double a2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a2) & ((-1L) << 27));
  1315.         final double a2Low      = a2 - a2High;
  1316.         final double b2High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b2) & ((-1L) << 27));
  1317.         final double b2Low      = b2 - b2High;

  1319.         // accurate multiplication a2 * b2
  1320.         final double prod2High  = a2 * b2;
  1321.         final double prod2Low   = a2Low * b2Low - (((prod2High - a2High * b2High) - a2Low * b2High) - a2High * b2Low);

  1323.         // split a3 and b3 as one 26 bits number and one 27 bits number
  1324.         final double a3High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a3) & ((-1L) << 27));
  1325.         final double a3Low      = a3 - a3High;
  1326.         final double b3High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b3) & ((-1L) << 27));
  1327.         final double b3Low      = b3 - b3High;

  1329.         // accurate multiplication a3 * b3
  1330.         final double prod3High  = a3 * b3;
  1331.         final double prod3Low   = a3Low * b3Low - (((prod3High - a3High * b3High) - a3Low * b3High) - a3High * b3Low);

  1333.         // split a4 and b4 as one 26 bits number and one 27 bits number
  1334.         final double a4High     = Double.longBitsToDouble(Double.doubleToRawLongBits(a4) & ((-1L) << 27));
  1335.         final double a4Low      = a4 - a4High;
  1336.         final double b4High     = Double.longBitsToDouble(Double.doubleToRawLongBits(b4) & ((-1L) << 27));
  1337.         final double b4Low      = b4 - b4High;

  1339.         // accurate multiplication a4 * b4
  1340.         final double prod4High  = a4 * b4;
  1341.         final double prod4Low   = a4Low * b4Low - (((prod4High - a4High * b4High) - a4Low * b4High) - a4High * b4Low);

  1343.         // accurate addition a1 * b1 + a2 * b2
  1344.         final double s12High    = prod1High + prod2High;
  1345.         final double s12Prime   = s12High - prod2High;
  1346.         final double s12Low     = (prod2High - (s12High - s12Prime)) + (prod1High - s12Prime);

  1348.         // accurate addition a1 * b1 + a2 * b2 + a3 * b3
  1349.         final double s123High   = s12High + prod3High;
  1350.         final double s123Prime  = s123High - prod3High;
  1351.         final double s123Low    = (prod3High - (s123High - s123Prime)) + (s12High - s123Prime);

  1353.         // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4
  1354.         final double s1234High  = s123High + prod4High;
  1355.         final double s1234Prime = s1234High - prod4High;
  1356.         final double s1234Low   = (prod4High - (s1234High - s1234Prime)) + (s123High - s1234Prime);

  1358.         // final rounding, s1234 may have suffered many cancellations, we try
  1359.         // to recover some bits from the extra words we have saved up to now
  1360.         double result = s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + s12Low + s123Low + s1234Low);

  1362.         if (Double.isNaN(result) || result == 0.0) {
  1363.             // either we have split infinite numbers or some coefficients were NaNs or signed zeros,
  1364.             // just rely on the naive implementation and let IEEE754 handle this
  1365.             // we do this for zeros too as we want to preserve the sign of zero (see issue #76)
  1366.             result = a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4;
  1367.         }

  1369.         return result;
  1370.     }

  1372.     /**
  1373.      * Returns true iff both arguments are null or have same dimensions and all
  1374.      * their elements are equal as defined by
  1375.      * {@link Precision#equals(float,float)}.
  1376.      *
  1377.      * @param x first array
  1378.      * @param y second array
  1379.      * @return true if the values are both null or have same dimension
  1380.      * and equal elements.
  1381.      */
  1382.     public static boolean equals(float[] x, float[] y) {
  1383.         if ((x == null) || (y == null)) {
  1384.             return !((x == null) ^ (y == null));
  1385.         }
  1386.         if (x.length != y.length) {
  1387.             return false;
  1388.         }
  1389.         for (int i = 0; i < x.length; ++i) {
  1390.             if (!Precision.equals(x[i], y[i])) {
  1391.                 return false;
  1392.             }
  1393.         }
  1394.         return true;
  1395.     }

  1397.     /**
  1398.      * Returns true iff both arguments are null or have same dimensions and all
  1399.      * their elements are equal as defined by
  1400.      * {@link Precision#equalsIncludingNaN(double,double) this method}.
  1401.      *
  1402.      * @param x first array
  1403.      * @param y second array
  1404.      * @return true if the values are both null or have same dimension and
  1405.      * equal elements
  1406.      */
  1407.     public static boolean equalsIncludingNaN(float[] x, float[] y) {
  1408.         if ((x == null) || (y == null)) {
  1409.             return !((x == null) ^ (y == null));
  1410.         }
  1411.         if (x.length != y.length) {
  1412.             return false;
  1413.         }
  1414.         for (int i = 0; i < x.length; ++i) {
  1415.             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
  1416.                 return false;
  1417.             }
  1418.         }
  1419.         return true;
  1420.     }

  1422.     /**
  1423.      * Returns {@code true} iff both arguments are {@code null} or have same
  1424.      * dimensions and all their elements are equal as defined by
  1425.      * {@link Precision#equals(double,double)}.
  1426.      *
  1427.      * @param x First array.
  1428.      * @param y Second array.
  1429.      * @return {@code true} if the values are both {@code null} or have same
  1430.      * dimension and equal elements.
  1431.      */
  1432.     public static boolean equals(double[] x, double[] y) {
  1433.         if ((x == null) || (y == null)) {
  1434.             return !((x == null) ^ (y == null));
  1435.         }
  1436.         if (x.length != y.length) {
  1437.             return false;
  1438.         }
  1439.         for (int i = 0; i < x.length; ++i) {
  1440.             if (!Precision.equals(x[i], y[i])) {
  1441.                 return false;
  1442.             }
  1443.         }
  1444.         return true;
  1445.     }

  1447.     /**
  1448.      * Returns {@code true} iff both arguments are {@code null} or have same
  1449.      * dimensions and all their elements are equal as defined by
  1450.      * {@link Object#equals(Object)}.
  1451.      *
  1452.      * @param <T> type of the field elements
  1453.      * @param x First array.
  1454.      * @param y Second array.
  1455.      * @return {@code true} if the values are both {@code null} or have same
  1456.      * dimension and equal elements.
  1457.      * @since 4.0
  1458.      */
  1459.     public static <T extends FieldElement<T>> boolean equals(T[] x, T[] y) {
  1460.         if ((x == null) || (y == null)) {
  1461.             return !((x == null) ^ (y == null));
  1462.         }
  1463.         if (x.length != y.length) {
  1464.             return false;
  1465.         }
  1466.         for (int i = 0; i < x.length; ++i) {
  1467.             if (!x[i].equals(y[i])) {
  1468.                 return false;
  1469.             }
  1470.         }
  1471.         return true;
  1472.     }

  1474.     /**
  1475.      * Returns {@code true} iff both arguments are {@code null} or have same
  1476.      * dimensions and all their elements are equal as defined by
  1477.      * {@link Precision#equalsIncludingNaN(double,double) this method}.
  1478.      *
  1479.      * @param x First array.
  1480.      * @param y Second array.
  1481.      * @return {@code true} if the values are both {@code null} or have same
  1482.      * dimension and equal elements.
  1483.      */
  1484.     public static boolean equalsIncludingNaN(double[] x, double[] y) {
  1485.         if ((x == null) || (y == null)) {
  1486.             return !((x == null) ^ (y == null));
  1487.         }
  1488.         if (x.length != y.length) {
  1489.             return false;
  1490.         }
  1491.         for (int i = 0; i < x.length; ++i) {
  1492.             if (!Precision.equalsIncludingNaN(x[i], y[i])) {
  1493.                 return false;
  1494.             }
  1495.         }
  1496.         return true;
  1497.     }

  1499.     /**
  1500.      * Returns {@code true} if both arguments are {@code null} or have same
  1501.      * dimensions and all their elements are equals.
  1502.      *
  1503.      * @param x First array.
  1504.      * @param y Second array.
  1505.      * @return {@code true} if the values are both {@code null} or have same
  1506.      * dimension and equal elements.
  1507.      */
  1508.     public static boolean equals(long[] x, long[] y) {
  1509.         if ((x == null) || (y == null)) {
  1510.             return !((x == null) ^ (y == null));
  1511.         }
  1512.         if (x.length != y.length) {
  1513.             return false;
  1514.         }
  1515.         for (int i = 0; i < x.length; ++i) {
  1516.             if (x[i] != y[i]) {
  1517.                 return false;
  1518.             }
  1519.         }
  1520.         return true;
  1521.     }

  1523.     /**
  1524.      * Returns {@code true} if both arguments are {@code null} or have same
  1525.      * dimensions and all their elements are equals.
  1526.      *
  1527.      * @param x First array.
  1528.      * @param y Second array.
  1529.      * @return {@code true} if the values are both {@code null} or have same
  1530.      * dimension and equal elements.
  1531.      */
  1532.     public static boolean equals(int[] x, int[] y) {
  1533.         if ((x == null) || (y == null)) {
  1534.             return !((x == null) ^ (y == null));
  1535.         }
  1536.         if (x.length != y.length) {
  1537.             return false;
  1538.         }
  1539.         for (int i = 0; i < x.length; ++i) {
  1540.             if (x[i] != y[i]) {
  1541.                 return false;
  1542.             }
  1543.         }
  1544.         return true;
  1545.     }

  1547.     /**
  1548.      * Returns {@code true} if both arguments are {@code null} or have same
  1549.      * dimensions and all their elements are equals.
  1550.      *
  1551.      * @param x First array.
  1552.      * @param y Second array.
  1553.      * @return {@code true} if the values are both {@code null} or have same
  1554.      * dimension and equal elements.
  1555.      */
  1556.     public static boolean equals(byte[] x, byte[] y) {
  1557.         if ((x == null) || (y == null)) {
  1558.             return !((x == null) ^ (y == null));
  1559.         }
  1560.         if (x.length != y.length) {
  1561.             return false;
  1562.         }
  1563.         for (int i = 0; i < x.length; ++i) {
  1564.             if (x[i] != y[i]) {
  1565.                 return false;
  1566.             }
  1567.         }
  1568.         return true;
  1569.     }

  1571.     /**
  1572.      * Returns {@code true} if both arguments are {@code null} or have same
  1573.      * dimensions and all their elements are equals.
  1574.      *
  1575.      * @param x First array.
  1576.      * @param y Second array.
  1577.      * @return {@code true} if the values are both {@code null} or have same
  1578.      * dimension and equal elements.
  1579.      */
  1580.     public static boolean equals(short[] x, short[] y) {
  1581.         if ((x == null) || (y == null)) {
  1582.             return !((x == null) ^ (y == null));
  1583.         }
  1584.         if (x.length != y.length) {
  1585.             return false;
  1586.         }
  1587.         for (int i = 0; i < x.length; ++i) {
  1588.             if (x[i] != y[i]) {
  1589.                 return false;
  1590.             }
  1591.         }
  1592.         return true;
  1593.     }

  1595.     /**
  1596.      * Normalizes an array to make it sum to a specified value.
  1597.      * Returns the result of the transformation
  1598.      * <pre>
  1599.      *    x ↦ x * normalizedSum / sum
  1600.      * </pre>
  1601.      * applied to each non-NaN element x of the input array, where sum is the
  1602.      * sum of the non-NaN entries in the input array.
  1603.      * <p>
  1604.      * Throws IllegalArgumentException if {@code normalizedSum} is infinite
  1605.      * or NaN and ArithmeticException if the input array contains any infinite elements
  1606.      * or sums to 0.
  1607.      * <p>
  1608.      * Ignores (i.e., copies unchanged to the output array) NaNs in the input array.
  1609.      * The input array is unchanged by this method.
  1610.      *
  1611.      * @param values Input array to be normalized
  1612.      * @param normalizedSum Target sum for the normalized array
  1613.      * @return the normalized array
  1614.      * @throws MathRuntimeException if the input array contains infinite
  1615.      * elements or sums to zero
  1616.      * @throws MathIllegalArgumentException if the target sum is infinite or {@code NaN}
  1617.      */
  1618.     public static double[] normalizeArray(double[] values, double normalizedSum)
  1619.         throws MathIllegalArgumentException, MathRuntimeException {
  1620.         if (Double.isInfinite(normalizedSum)) {
  1621.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NORMALIZE_INFINITE);
  1622.         }
  1623.         if (Double.isNaN(normalizedSum)) {
  1624.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NORMALIZE_NAN);
  1625.         }
  1626.         double sum = 0d;
  1627.         final int len = values.length;
  1628.         double[] out = new double[len];
  1629.         for (int i = 0; i < len; i++) {
  1630.             if (Double.isInfinite(values[i])) {
  1631.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.INFINITE_ARRAY_ELEMENT, values[i], i);
  1632.             }
  1633.             if (!Double.isNaN(values[i])) {
  1634.                 sum += values[i];
  1635.             }
  1636.         }
  1637.         if (sum == 0) {
  1638.             throw new MathRuntimeException(LocalizedCoreFormats.ARRAY_SUMS_TO_ZERO);
  1639.         }
  1640.         for (int i = 0; i < len; i++) {
  1641.             if (Double.isNaN(values[i])) {
  1642.                 out[i] = Double.NaN;
  1643.             } else {
  1644.                 out[i] = values[i] * normalizedSum / sum;
  1645.             }
  1646.         }
  1647.         return out;
  1648.     }

  1650.     /** Build an array of elements.
  1651.      * <p>
  1652.      * Arrays are filled with {@code field.getZero()}
  1653.      *
  1654.      * @param <T> the type of the field elements
  1655.      * @param field field to which array elements belong
  1656.      * @param length of the array
  1657.      * @return a new array
  1658.      */
  1659.     public static <T extends FieldElement<T>> T[] buildArray(final Field<T> field, final int length) {
  1660.         @SuppressWarnings("unchecked") // OK because field must be correct class
  1661.         T[] array = (T[]) Array.newInstance(field.getRuntimeClass(), length);
  1662.         Arrays.fill(array, field.getZero());
  1663.         return array;
  1664.     }

  1666.     /** Build a double dimension array of elements.
  1667.      * <p>
  1668.      * Arrays are filled with {@code field.getZero()}
  1669.      *
  1670.      * @param <T> the type of the field elements
  1671.      * @param field field to which array elements belong
  1672.      * @param rows number of rows in the array
  1673.      * @param columns number of columns (may be negative to build partial
  1674.      * arrays in the same way {@code new Field[rows][]} works)
  1675.      * @return a new array
  1676.      */
  1677.     @SuppressWarnings("unchecked")
  1678.     public static <T extends FieldElement<T>> T[][] buildArray(final Field<T> field, final int rows, final int columns) {
  1679.         final T[][] array;
  1680.         if (columns < 0) {
  1681.             T[] dummyRow = buildArray(field, 0);
  1682.             array = (T[][]) Array.newInstance(dummyRow.getClass(), rows);
  1683.         } else {
  1684.             array = (T[][]) Array.newInstance(field.getRuntimeClass(), rows, columns);
  1685.             for (int i = 0; i < rows; ++i) {
  1686.                 Arrays.fill(array[i], field.getZero());
  1687.             }
  1688.         }
  1689.         return array;
  1690.     }

  1692.     /** Build a triple dimension array of elements.
  1693.      * <p>
  1694.      * Arrays are filled with {@code field.getZero()}
  1695.      *
  1696.      * @param <T> the type of the field elements
  1697.      * @param field field to which array elements belong
  1698.      * @param l1 number of elements along first dimension
  1699.      * @param l2 number of elements along second dimension
  1700.      * @param l3 number of elements along third dimension (may be negative to build partial
  1701.      * arrays in the same way {@code new Field[l1][l2][]} works)
  1702.      * @return a new array
  1703.      * @since 1.4
  1704.      */
  1705.     @SuppressWarnings("unchecked")
  1706.     public static <T extends FieldElement<T>> T[][][] buildArray(final Field<T> field, final int l1, final int l2, final int l3) {
  1707.         final T[][][] array;
  1708.         if (l3 < 0) {
  1709.             T[] dummyRow = buildArray(field, 0);
  1710.             array = (T[][][]) Array.newInstance(dummyRow.getClass(), l1, l2);
  1711.         } else {
  1712.             array = (T[][][]) Array.newInstance(field.getRuntimeClass(), l1, l2, l3);
  1713.             for (int i = 0; i < l1; ++i) {
  1714.                 for (int j = 0; j < l2; ++j) {
  1715.                     Arrays.fill(array[i][j], field.getZero());
  1716.                 }
  1717.             }
  1718.         }
  1719.         return array;
  1720.     }

  1722.     /**
  1723.      * Calculates the <a href="http://en.wikipedia.org/wiki/Convolution">
  1724.      * convolution</a> between two sequences.
  1725.      * <p>
  1726.      * The solution is obtained via straightforward computation of the
  1727.      * convolution sum (and not via FFT). Whenever the computation needs
  1728.      * an element that would be located at an index outside the input arrays,
  1729.      * the value is assumed to be zero.
  1730.      *
  1731.      * @param x First sequence.
  1732.      * Typically, this sequence will represent an input signal to a system.
  1733.      * @param h Second sequence.
  1734.      * Typically, this sequence will represent the impulse response of the system.
  1735.      * @return the convolution of {@code x} and {@code h}.
  1736.      * This array's length will be {@code x.length + h.length - 1}.
  1737.      * @throws NullArgumentException if either {@code x} or {@code h} is {@code null}.
  1738.      * @throws MathIllegalArgumentException if either {@code x} or {@code h} is empty.
  1739.      */
  1740.     public static double[] convolve(double[] x, double[] h)
  1741.         throws MathIllegalArgumentException, NullArgumentException {
  1742.         MathUtils.checkNotNull(x);
  1743.         MathUtils.checkNotNull(h);

  1745.         final int xLen = x.length;
  1746.         final int hLen = h.length;

  1748.         if (xLen == 0 || hLen == 0) {
  1749.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  1750.         }

  1752.         // initialize the output array
  1753.         final int totalLength = xLen + hLen - 1;
  1754.         final double[] y = new double[totalLength];

  1756.         // straightforward implementation of the convolution sum
  1757.         for (int n = 0; n < totalLength; n++) {
  1758.             double yn = 0;
  1759.             int k = FastMath.max(0, n + 1 - xLen);
  1760.             int j = n - k;
  1761.             while (k < hLen && j >= 0) {
  1762.                 yn += x[j--] * h[k++];
  1763.             }
  1764.             y[n] = yn;
  1765.         }

  1767.         return y;
  1768.     }

  1770.     /**
  1771.      * Specification for indicating that some operation applies
  1772.      * before or after a given index.
  1773.      */
  1774.     public enum Position {
  1775.         /** Designates the beginning of the array (near index 0). */
  1776.         HEAD,
  1777.         /** Designates the end of the array. */
  1778.         TAIL
  1779.     }

  1781.     /**
  1782.      * Shuffle the entries of the given array.
  1783.      * The {@code start} and {@code pos} parameters select which portion
  1784.      * of the array is randomized and which is left untouched.
  1785.      *
  1786.      * @see #shuffle(int[],int,Position,RandomGenerator)
  1787.      *
  1788.      * @param list Array whose entries will be shuffled (in-place).
  1789.      * @param start Index at which shuffling begins.
  1790.      * @param pos Shuffling is performed for index positions between
  1791.      * {@code start} and either the end (if {@link Position#TAIL})
  1792.      * or the beginning (if {@link Position#HEAD}) of the array.
  1793.      */
  1794.     public static void shuffle(int[] list,
  1795.                                int start,
  1796.                                Position pos) {
  1797.         shuffle(list, start, pos, new Well19937c());
  1798.     }

  1800.     /**
  1801.      * Shuffle the entries of the given array, using the
  1802.      * <a href="https://en.wikipedia.org/wiki/Fisher-Yates_shuffle#The_modern_algorithm">
  1803.      * Fisher–Yates</a> algorithm.
  1804.      * The {@code start} and {@code pos} parameters select which portion
  1805.      * of the array is randomized and which is left untouched.
  1806.      *
  1807.      * @param list Array whose entries will be shuffled (in-place).
  1808.      * @param start Index at which shuffling begins.
  1809.      * @param pos Shuffling is performed for index positions between
  1810.      * {@code start} and either the end (if {@link Position#TAIL})
  1811.      * or the beginning (if {@link Position#HEAD}) of the array.
  1812.      * @param rng Random number generator.
  1813.      */
  1814.     public static void shuffle(int[] list,
  1815.                                int start,
  1816.                                Position pos,
  1817.                                RandomGenerator rng) {
  1818.         switch (pos) {
  1819.         case TAIL:
  1820.             for (int i = list.length - 1; i > start; i--) {
  1821.                 final int target = start + rng.nextInt(i - start + 1);
  1822.                 final int temp = list[target];
  1823.                 list[target] = list[i];
  1824.                 list[i] = temp;
  1825.             }
  1826.             break;

  1828.         case HEAD:
  1829.             for (int i = 0; i < start; i++) {
  1830.                 final int target = i + rng.nextInt(start - i + 1);
  1831.                 final int temp = list[target];
  1832.                 list[target] = list[i];
  1833.                 list[i] = temp;
  1834.             }
  1835.             break;

  1837.         default:
  1838.             throw MathRuntimeException.createInternalError(); // Should never happen.
  1839.         }
  1840.     }

  1842.     /**
  1843.      * Shuffle the entries of the given array.
  1844.      *
  1845.      * @see #shuffle(int[],int,Position,RandomGenerator)
  1846.      *
  1847.      * @param list Array whose entries will be shuffled (in-place).
  1848.      * @param rng Random number generator.
  1849.      */
  1850.     public static void shuffle(int[] list,
  1851.                                RandomGenerator rng) {
  1852.         shuffle(list, 0, Position.TAIL, rng);
  1853.     }

  1855.     /**
  1856.      * Shuffle the entries of the given array.
  1857.      *
  1858.      * @see #shuffle(int[],int,Position,RandomGenerator)
  1859.      *
  1860.      * @param list Array whose entries will be shuffled (in-place).
  1861.      */
  1862.     public static void shuffle(int[] list) {
  1863.         shuffle(list, new Well19937c());
  1864.     }

  1866.     /**
  1867.      * Returns an array representing the natural number {@code n}.
  1868.      *
  1869.      * @param n Natural number.
  1870.      * @return an array whose entries are the numbers 0, 1, ..., {@code n}-1.
  1871.      * If {@code n == 0}, the returned array is empty.
  1872.      */
  1873.     public static int[] natural(int n) {
  1874.         return sequence(n, 0, 1);
  1875.     }

  1877.     /**
  1878.      * Returns an array of {@code size} integers starting at {@code start},
  1879.      * skipping {@code stride} numbers.
  1880.      *
  1881.      * @param size Natural number.
  1882.      * @param start Natural number.
  1883.      * @param stride Natural number.
  1884.      * @return an array whose entries are the numbers
  1885.      * {@code start, start + stride, ..., start + (size - 1) * stride}.
  1886.      * If {@code size == 0}, the returned array is empty.
  1887.      */
  1888.     public static int[] sequence(int size,
  1889.                                  int start,
  1890.                                  int stride) {
  1891.         final int[] a = new int[size];
  1892.         for (int i = 0; i < size; i++) {
  1893.             a[i] = start + i * stride;
  1894.         }
  1895.         return a;
  1896.     }

  1898.     /**
  1899.      * This method is used
  1900.      * to verify that the input parameters designate a subarray of positive length.
  1901.      * <ul>
  1902.      * <li>returns {@code true} iff the parameters designate a subarray of
  1903.      * positive length</li>
  1904.      * <li>throws {@code MathIllegalArgumentException} if the array is null or
  1905.      * or the indices are invalid</li>
  1906.      * <li>returns {@code false} if the array is non-null, but
  1907.      * {@code length} is 0.</li>
  1908.      * </ul>
  1909.      *
  1910.      * @param values the input array
  1911.      * @param begin index of the first array element to include
  1912.      * @param length the number of elements to include
  1913.      * @return true if the parameters are valid and designate a subarray of positive length
  1914.      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
  1915.      */
  1916.     public static boolean verifyValues(final double[] values, final int begin, final int length)
  1917.             throws MathIllegalArgumentException {
  1918.         return verifyValues(values, begin, length, false);
  1919.     }

  1921.     /**
  1922.      * This method is used
  1923.      * to verify that the input parameters designate a subarray of positive length.
  1924.      * <ul>
  1925.      * <li>returns {@code true} iff the parameters designate a subarray of
  1926.      * non-negative length</li>
  1927.      * <li>throws {@code IllegalArgumentException} if the array is null or
  1928.      * or the indices are invalid</li>
  1929.      * <li>returns {@code false} if the array is non-null, but
  1930.      * {@code length} is 0 unless {@code allowEmpty} is {@code true}</li>
  1931.      * </ul>
  1932.      *
  1933.      * @param values the input array
  1934.      * @param begin index of the first array element to include
  1935.      * @param length the number of elements to include
  1936.      * @param allowEmpty if <code>true</code> then zero length arrays are allowed
  1937.      * @return true if the parameters are valid
  1938.      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
  1939.      */
  1940.     public static boolean verifyValues(final double[] values, final int begin,
  1941.             final int length, final boolean allowEmpty) throws MathIllegalArgumentException {

  1943.         MathUtils.checkNotNull(values, LocalizedCoreFormats.INPUT_ARRAY);

  1945.         if (begin < 0) {
  1946.             throw new MathIllegalArgumentException(LocalizedCoreFormats.START_POSITION, begin);
  1947.         }

  1949.         if (length < 0) {
  1950.             throw new MathIllegalArgumentException(LocalizedCoreFormats.LENGTH, length);
  1951.         }

  1953.         if (begin + length > values.length) {
  1954.             throw new MathIllegalArgumentException(LocalizedCoreFormats.SUBARRAY_ENDS_AFTER_ARRAY_END,
  1955.                     begin + length, values.length, true);
  1956.         }

  1958.         if (length == 0 && !allowEmpty) {
  1959.             return false;
  1960.         }

  1962.         return true;
  1963.     }

  1965.     /**
  1966.      * This method is used
  1967.      * to verify that the begin and length parameters designate a subarray of positive length
  1968.      * and the weights are all non-negative, non-NaN, finite, and not all zero.
  1969.      * <ul>
  1970.      * <li>returns {@code true} iff the parameters designate a subarray of
  1971.      * positive length and the weights array contains legitimate values.</li>
  1972.      * <li>throws {@code IllegalArgumentException} if any of the following are true:
  1973.      * <ul><li>the values array is null</li>
  1974.      *     <li>the weights array is null</li>
  1975.      *     <li>the weights array does not have the same length as the values array</li>
  1976.      *     <li>the weights array contains one or more infinite values</li>
  1977.      *     <li>the weights array contains one or more NaN values</li>
  1978.      *     <li>the weights array contains negative values</li>
  1979.      *     <li>the start and length arguments do not determine a valid array</li></ul>
  1980.      * </li>
  1981.      * <li>returns {@code false} if the array is non-null, but {@code length} is 0.</li>
  1982.      * </ul>
  1983.      *
  1984.      * @param values the input array
  1985.      * @param weights the weights array
  1986.      * @param begin index of the first array element to include
  1987.      * @param length the number of elements to include
  1988.      * @return true if the parameters are valid and designate a subarray of positive length
  1989.      * @throws MathIllegalArgumentException if the indices are invalid or the array is null
  1990.      */
  1991.     public static boolean verifyValues(
  1992.         final double[] values,
  1993.         final double[] weights,
  1994.         final int begin,
  1995.         final int length) throws MathIllegalArgumentException {
  1996.         return verifyValues(values, weights, begin, length, false);
  1997.     }

  1999.     /**
  2000.      * This method is used
  2001.      * to verify that the begin and length parameters designate a subarray of positive length
  2002.      * and the weights are all non-negative, non-NaN, finite, and not all zero.
  2003.      * <ul>
  2004.      * <li>returns {@code true} iff the parameters designate a subarray of
  2005.      * non-negative length and the weights array contains legitimate values.</li>
  2006.      * <li>throws {@code MathIllegalArgumentException} if any of the following are true:
  2007.      * <ul><li>the values array is null</li>
  2008.      *     <li>the weights array is null</li>
  2009.      *     <li>the weights array does not have the same length as the values array</li>
  2010.      *     <li>the weights array contains one or more infinite values</li>
  2011.      *     <li>the weights array contains one or more NaN values</li>
  2012.      *     <li>the weights array contains negative values</li>
  2013.      *     <li>the start and length arguments do not determine a valid array</li></ul>
  2014.      * </li>
  2015.      * <li>returns {@code false} if the array is non-null, but
  2016.      * {@code length} is 0 unless {@code allowEmpty} is {@code true}</li>
  2017.      * </ul>
  2018.      *
  2019.      * @param values the input array.
  2020.      * @param weights the weights array.
  2021.      * @param begin index of the first array element to include.
  2022.      * @param length the number of elements to include.
  2023.      * @param allowEmpty if {@code true} than allow zero length arrays to pass.
  2024.      * @return {@code true} if the parameters are valid.
  2025.      * @throws NullArgumentException if either of the arrays are null
  2026.      * @throws MathIllegalArgumentException if the array indices are not valid,
  2027.      * the weights array contains NaN, infinite or negative elements, or there
  2028.      * are no positive weights.
  2029.      */
  2030.     public static boolean verifyValues(final double[] values, final double[] weights,
  2031.             final int begin, final int length, final boolean allowEmpty) throws MathIllegalArgumentException {

  2033.         MathUtils.checkNotNull(weights, LocalizedCoreFormats.INPUT_ARRAY);
  2034.         MathUtils.checkNotNull(values, LocalizedCoreFormats.INPUT_ARRAY);

  2036.         checkEqualLength(weights, values);

  2038.         boolean containsPositiveWeight = false;
  2039.         for (int i = begin; i < begin + length; i++) {
  2040.             final double weight = weights[i];
  2041.             if (Double.isNaN(weight)) {
  2042.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NAN_ELEMENT_AT_INDEX, i);
  2043.             }
  2044.             if (Double.isInfinite(weight)) {
  2045.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.INFINITE_ARRAY_ELEMENT, weight, i);
  2046.             }
  2047.             if (weight < 0) {
  2048.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NEGATIVE_ELEMENT_AT_INDEX, i, weight);
  2049.             }
  2050.             if (!containsPositiveWeight && weight > 0.0) {
  2051.                 containsPositiveWeight = true;
  2052.             }
  2053.         }

  2055.         if (!containsPositiveWeight) {
  2056.             throw new MathIllegalArgumentException(LocalizedCoreFormats.WEIGHT_AT_LEAST_ONE_NON_ZERO);
  2057.         }

  2059.         return verifyValues(values, begin, length, allowEmpty);
  2060.     }

  2062.     /**
  2063.      * Concatenates a sequence of arrays. The return array consists of the
  2064.      * entries of the input arrays concatenated in the order they appear in
  2065.      * the argument list.  Null arrays cause NullPointerExceptions; zero
  2066.      * length arrays are allowed (contributing nothing to the output array).
  2067.      *
  2068.      * @param x list of double[] arrays to concatenate
  2069.      * @return a new array consisting of the entries of the argument arrays
  2070.      * @throws NullPointerException if any of the arrays are null
  2071.      */
  2072.     public static double[] concatenate(double[]... x) {
  2073.         int combinedLength = 0;
  2074.         for (double[] a : x) {
  2075.             combinedLength += a.length;
  2076.         }
  2077.         int offset = 0;
  2078.         final double[] combined = new double[combinedLength];
  2079.         for (double[] doubles : x) {
  2080.             final int curLength = doubles.length;
  2081.             System.arraycopy(doubles, 0, combined, offset, curLength);
  2082.             offset += curLength;
  2083.         }
  2084.         return combined;
  2085.     }

  2087.     /**
  2088.      * Returns an array consisting of the unique values in {@code data}.
  2089.      * The return array is sorted in descending order.  Empty arrays
  2090.      * are allowed, but null arrays result in NullPointerException.
  2091.      * Infinities are allowed.  NaN values are allowed with maximum
  2092.      * sort order - i.e., if there are NaN values in {@code data},
  2093.      * {@code Double.NaN} will be the first element of the output array,
  2094.      * even if the array also contains {@code Double.POSITIVE_INFINITY}.
  2095.      *
  2096.      * @param data array to scan
  2097.      * @return descending list of values included in the input array
  2098.      * @throws NullPointerException if data is null
  2099.      */
  2100.     public static double[] unique(double[] data) {
  2101.         NavigableSet<Double> values = new TreeSet<>();
  2102.         for (double datum : data) {
  2103.             values.add(datum);
  2104.         }
  2105.         final int count = values.size();
  2106.         final double[] out = new double[count];
  2107.         Iterator<Double> iterator = values.descendingIterator();
  2108.         int i = 0;
  2109.         while (iterator.hasNext()) {
  2110.             out[i++] = iterator.next();
  2111.         }
  2112.         return out;
  2113.     }
  2114. }
Created with JaCoCo 0.8.12.202403310830