MathArrays { public static double linearCombination(final double[] a, final double[] b) throws DimensionMismatchException { final int len = a.length; if (len != b.length) { throw new DimensionMismatchException(len, b.length); } if (len == 1) { return a[0] * b[0]; } final double[] prodHigh = new double[len]; double prodLowSum = 0; for (int i = 0; i < len; i++) { final double ai = a[i]; final double ca = SPLIT_FACTOR * ai; final double aHigh = ca - (ca - ai); final double aLow = ai - aHigh; final double bi = b[i]; final double cb = SPLIT_FACTOR * bi; final double bHigh = cb - (cb - bi); final double bLow = bi - bHigh; prodHigh[i] = ai * bi; final double prodLow = aLow * bLow - (((prodHigh[i] - aHigh * bHigh) - aLow * bHigh) - aHigh * bLow); prodLowSum += prodLow; } final double prodHighCur = prodHigh[0]; double prodHighNext = prodHigh[1]; double sHighPrev = prodHighCur + prodHighNext; double sPrime = sHighPrev - prodHighNext; double sLowSum = (prodHighNext - (sHighPrev - sPrime)) + (prodHighCur - sPrime); final int lenMinusOne = len - 1; for (int i = 1; i < lenMinusOne; i++) { prodHighNext = prodHigh[i + 1]; final double sHighCur = sHighPrev + prodHighNext; sPrime = sHighCur - prodHighNext; sLowSum += (prodHighNext - (sHighCur - sPrime)) + (sHighPrev - sPrime); sHighPrev = sHighCur; } double result = sHighPrev + (prodLowSum + sLowSum); if (Double.isNaN(result)) { result = 0; for (int i = 0; i < len; ++i) { result += a[i] * b[i]; } } return result; } private  MathArrays(); static double[] scale(double val, final double[] arr); static void scaleInPlace(double val, final double[] arr); static double[] ebeAdd(double[] a, double[] b); static double[] ebeSubtract(double[] a, double[] b); static double[] ebeMultiply(double[] a, double[] b); static double[] ebeDivide(double[] a, double[] b); static double distance1(double[] p1, double[] p2); static int distance1(int[] p1, int[] p2); static double distance(double[] p1, double[] p2); static double distance(int[] p1, int[] p2); static double distanceInf(double[] p1, double[] p2); static int distanceInf(int[] p1, int[] p2); static boolean isMonotonic(T[] val,
                                      OrderDirection dir,
                                      boolean strict); static boolean isMonotonic(double[] val, OrderDirection dir, boolean strict); static boolean checkOrder(double[] val, OrderDirection dir,
                                     boolean strict, boolean abort); static void checkOrder(double[] val, OrderDirection dir,
                                  boolean strict); static void checkOrder(double[] val); static void checkRectangular(final long[][] in); static void checkPositive(final double[] in); static void checkNotNaN(final double[] in); static void checkNonNegative(final long[] in); static void checkNonNegative(final long[][] in); static double safeNorm(double[] v); static void sortInPlace(double[] x, double[] ... yList); static void sortInPlace(double[] x,
                                   final OrderDirection dir,
                                   double[] ... yList); static int[] copyOf(int[] source); static double[] copyOf(double[] source); static int[] copyOf(int[] source, int len); static double[] copyOf(double[] source, int len); static double[] copyOfRange(double[] source, int from, int to); static double linearCombination(final double[] a, final double[] b); static double linearCombination(final double a1, final double b1,
                                           final double a2, final double b2); static double linearCombination(final double a1, final double b1,
                                           final double a2, final double b2,
                                           final double a3, final double b3); static double linearCombination(final double a1, final double b1,
                                           final double a2, final double b2,
                                           final double a3, final double b3,
                                           final double a4, final double b4); static boolean equals(float[] x, float[] y); static boolean equalsIncludingNaN(float[] x, float[] y); static boolean equals(double[] x, double[] y); static boolean equalsIncludingNaN(double[] x, double[] y); static double[] normalizeArray(double[] values, double normalizedSum); static T[] buildArray(final Field<T> field, final int length); @SuppressWarnings("unchecked") static T[][] buildArray(final Field<T> field, final int rows, final int columns); static double[] convolve(double[] x, double[] h); static void shuffle(int[] list,
                               int start,
                               Position pos); static void shuffle(int[] list,
                               int start,
                               Position pos,
                               RandomGenerator rng); static void shuffle(int[] list,
                               RandomGenerator rng); static void shuffle(int[] list); static int[] natural(int n); static int[] sequence(int size,
                                 int start,
                                 int stride); static boolean verifyValues(final double[] values, final int begin, final int length); static boolean verifyValues(final double[] values, final int begin,
            final int length, final boolean allowEmpty); static boolean verifyValues(
        final double[] values,
        final double[] weights,
        final int begin,
        final int length); static boolean verifyValues(final double[] values, final double[] weights,
            final int begin, final int length, final boolean allowEmpty);  }
@Test public void testLinearCombinationWithSingleElementArray() { final double[] a = { 1.23456789 }; final double[] b = { 98765432.1 }; Assert.assertEquals(a[0] * b[0], MathArrays.linearCombination(a, b), 0d); }
@Test public void testLinearCombination1() { final double[] a = new double[] { -1321008684645961.0 / 268435456.0, -5774608829631843.0 / 268435456.0, -7645843051051357.0 / 8589934592.0 }; final double[] b = new double[] { -5712344449280879.0 / 2097152.0, -4550117129121957.0 / 2097152.0, 8846951984510141.0 / 131072.0 }; final double abSumInline = MathArrays.linearCombination(a[0], b[0], a[1], b[1], a[2], b[2]); final double abSumArray = MathArrays.linearCombination(a, b); Assert.assertEquals(abSumInline, abSumArray, 0); Assert.assertEquals(-1.8551294182586248737720779899, abSumInline, 1.0e-15); final double naive = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; Assert.assertTrue(FastMath.abs(naive - abSumInline) > 1.5); }
@Test public void testLinearCombination2() { Well1024a random = new Well1024a(553267312521321234l); for (int i = 0; i < 10000; ++i) { final double ux = 1e17 * random.nextDouble(); final double uy = 1e17 * random.nextDouble(); final double uz = 1e17 * random.nextDouble(); final double vx = 1e17 * random.nextDouble(); final double vy = 1e17 * random.nextDouble(); final double vz = 1e17 * random.nextDouble(); final double sInline = MathArrays.linearCombination(ux, vx, uy, vy, uz, vz); final double sArray = MathArrays.linearCombination(new double[] {ux, uy, uz}, new double[] {vx, vy, vz}); Assert.assertEquals(sInline, sArray, 0); } }
@Test public void testLinearCombinationInfinite() { final double[][] a = new double[][] { { 1, 2, 3, 4}, { 1, Double.POSITIVE_INFINITY, 3, 4}, { 1, 2, Double.POSITIVE_INFINITY, 4}, { 1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY}, { 1, 2, 3, 4}, { 1, 2, 3, 4}, { 1, 2, 3, 4}, { 1, 2, 3, 4} }; final double[][] b = new double[][] { { 1, -2, 3, 4}, { 1, -2, 3, 4}, { 1, -2, 3, 4}, { 1, -2, 3, 4}, { 1, Double.POSITIVE_INFINITY, 3, 4}, { 1, -2, Double.POSITIVE_INFINITY, 4}, { 1, Double.POSITIVE_INFINITY, 3, Double.NEGATIVE_INFINITY}, { Double.NaN, -2, 3, 4} }; Assert.assertEquals(-3, MathArrays.linearCombination(a[0][0], b[0][0], a[0][1], b[0][1]), 1.0e-10); Assert.assertEquals(6, MathArrays.linearCombination(a[0][0], b[0][0], a[0][1], b[0][1], a[0][2], b[0][2]), 1.0e-10); Assert.assertEquals(22, MathArrays.linearCombination(a[0][0], b[0][0], a[0][1], b[0][1], a[0][2], b[0][2], a[0][3], b[0][3]), 1.0e-10); Assert.assertEquals(22, MathArrays.linearCombination(a[0], b[0]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[1][0], b[1][0], a[1][1], b[1][1]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[1][0], b[1][0], a[1][1], b[1][1], a[1][2], b[1][2]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[1][0], b[1][0], a[1][1], b[1][1], a[1][2], b[1][2], a[1][3], b[1][3]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[1], b[1]), 1.0e-10); Assert.assertEquals(-3, MathArrays.linearCombination(a[2][0], b[2][0], a[2][1], b[2][1]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[2][0], b[2][0], a[2][1], b[2][1], a[2][2], b[2][2]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[2][0], b[2][0], a[2][1], b[2][1], a[2][2], b[2][2], a[2][3], b[2][3]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[2], b[2]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[3][0], b[3][0], a[3][1], b[3][1]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[3][0], b[3][0], a[3][1], b[3][1], a[3][2], b[3][2]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[3][0], b[3][0], a[3][1], b[3][1], a[3][2], b[3][2], a[3][3], b[3][3]), 1.0e-10); Assert.assertEquals(Double.NEGATIVE_INFINITY, MathArrays.linearCombination(a[3], b[3]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[4][0], b[4][0], a[4][1], b[4][1]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[4][0], b[4][0], a[4][1], b[4][1], a[4][2], b[4][2]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[4][0], b[4][0], a[4][1], b[4][1], a[4][2], b[4][2], a[4][3], b[4][3]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[4], b[4]), 1.0e-10); Assert.assertEquals(-3, MathArrays.linearCombination(a[5][0], b[5][0], a[5][1], b[5][1]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[5][0], b[5][0], a[5][1], b[5][1], a[5][2], b[5][2]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[5][0], b[5][0], a[5][1], b[5][1], a[5][2], b[5][2], a[5][3], b[5][3]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[5], b[5]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[6][0], b[6][0], a[6][1], b[6][1]), 1.0e-10); Assert.assertEquals(Double.POSITIVE_INFINITY, MathArrays.linearCombination(a[6][0], b[6][0], a[6][1], b[6][1], a[6][2], b[6][2]), 1.0e-10); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[6][0], b[6][0], a[6][1], b[6][1], a[6][2], b[6][2], a[6][3], b[6][3]))); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[6], b[6]))); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[7][0], b[7][0], a[7][1], b[7][1]))); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[7][0], b[7][0], a[7][1], b[7][1], a[7][2], b[7][2]))); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[7][0], b[7][0], a[7][1], b[7][1], a[7][2], b[7][2], a[7][3], b[7][3]))); Assert.assertTrue(Double.isNaN(MathArrays.linearCombination(a[7], b[7]))); }