Return-Path: X-Original-To: apmail-commons-commits-archive@minotaur.apache.org Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 536CD7E9D for ; Sat, 6 Aug 2011 00:14:12 +0000 (UTC) Received: (qmail 23990 invoked by uid 500); 6 Aug 2011 00:14:11 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 23728 invoked by uid 500); 6 Aug 2011 00:14:10 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 23721 invoked by uid 99); 6 Aug 2011 00:14:10 -0000 Received: from athena.apache.org (HELO athena.apache.org) (140.211.11.136) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 06 Aug 2011 00:14:10 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=5.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 06 Aug 2011 00:14:09 +0000 Received: from eris.apache.org (localhost [127.0.0.1]) by eris.apache.org (Postfix) with ESMTP id 81D6723889B1 for ; Sat, 6 Aug 2011 00:13:49 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r1154416 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/util/MathUtils.java test/java/org/apache/commons/math/util/MathUtilsTest.java Date: Sat, 06 Aug 2011 00:13:49 -0000 To: commits@commons.apache.org From: erans@apache.org X-Mailer: svnmailer-1.0.8 Message-Id: <20110806001349.81D6723889B1@eris.apache.org> Author: erans Date: Sat Aug 6 00:13:49 2011 New Revision: 1154416 URL: http://svn.apache.org/viewvc?rev=1154416&view=rev Log: Array version of "linearCombination". Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java?rev=1154416&r1=1154415&r2=1154416&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java Sat Aug 6 00:13:49 2011 @@ -2608,4 +2608,64 @@ public final class MathUtils { } + /** + * Compute a linear combination accurately. + * This method computes the sum of the products + * a_i b_i to high accuracy. + * It does so by using specific multiplication and addition algorithms to + * preserve accuracy and reduce cancellation effects. + *
+ * It is based on the 2005 paper + * + * Accurate Sum and Dot Product by Takeshi Ogita, Siegfried M. Rump, + * and Shin'ichi Oishi published in SIAM J. Sci. Comput. + * + * @param a Factors. + * @param b Factors. + * @return Σ_i a_i b_i. + */ + public static double linearCombination(final double[] a, final double[] b) { + final int len = a.length; + if (len != b.length) { + throw new DimensionMismatchException(len, b.length); + } + + 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 int lenMinusOne = len - 1; + final double[] sHigh = new double[lenMinusOne]; + + sHigh[0] = prodHigh[0] + prodHigh[1]; + double sPrime = sHigh[0] - prodHigh[1]; + double sLowSum = (prodHigh[1] - (sHigh[0] - sPrime)) + (prodHigh[0] - sPrime); + + for (int i = 1; i < lenMinusOne; i++) { + final int prev = i - 1; + final int next = i + 1; + sHigh[i] = sHigh[prev] + prodHigh[next]; + sPrime = sHigh[i] - prodHigh[next]; + sLowSum += (prodHigh[next] - (sHigh[i] - sPrime)) + (sHigh[prev] - sPrime); + } + + return sHigh[lenMinusOne - 1] + (prodLowSum + sLowSum); + } } Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java?rev=1154416&r1=1154415&r2=1154416&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/util/MathUtilsTest.java Sat Aug 6 00:13:49 2011 @@ -1841,4 +1841,25 @@ public final class MathUtilsTest { // Expected. } } + + @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 = MathUtils.linearCombination(a[0], b[0], + a[1], b[1], + a[2], b[2]); + final double abSumArray = MathUtils.linearCombination(a, b); + + Assert.assertEquals(abSumInline, abSumArray, 0); + } }