I am +1 on the idea here and am wondering if we should think a little more radically even. Greg has - correctly, IMO - complained that the current structure of the linear package makes it awkward to implement optimized operations for some classes of matrices and operations. I a wonder whether it might not make sense to encapsulate high-performance, array-based implementations of core operations in utils classes that can be used like the methods below. I remember back in the very early days of [math] suggesting this for statistical algorithms and being talked out of it in favor of a more classical OO approach. I wonder if by organizing things correctly, we might not be able to have the best of both worlds.
Phil On 8/5/11 8:01 AM, l...@apache.org wrote: > Author: luc > Date: Fri Aug 5 15:01:49 2011 > New Revision: 1154250 > > URL: http://svn.apache.org/viewvc?rev=1154250&view=rev > Log: > Added a few linearCombination utility methods in MathUtils to compute > accurately > linear combinations a1.b1 + a2.b2 + ... + an.bn taking great care to > compensate > for cancellation effects. This both improves and simplify several methods in > euclidean geometry classes, including linear constructors, dot product and > cross > product. > > Modified: > > commons/proper/math/trunk/src/main/java/org/apache/commons/math/util/MathUtils.java > commons/proper/math/trunk/src/site/xdoc/changes.xml > > 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=1154250&r1=1154249&r2=1154250&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 > Fri Aug 5 15:01:49 2011 > @@ -19,22 +19,22 @@ package org.apache.commons.math.util; > > import java.math.BigDecimal; > import java.math.BigInteger; > -import java.util.Arrays; > -import java.util.List; > import java.util.ArrayList; > -import java.util.Comparator; > +import java.util.Arrays; > import java.util.Collections; > +import java.util.Comparator; > +import java.util.List; > > -import org.apache.commons.math.exception.util.Localizable; > -import org.apache.commons.math.exception.util.LocalizedFormats; > -import org.apache.commons.math.exception.NonMonotonousSequenceException; > import org.apache.commons.math.exception.DimensionMismatchException; > -import org.apache.commons.math.exception.NullArgumentException; > -import org.apache.commons.math.exception.NotPositiveException; > import org.apache.commons.math.exception.MathArithmeticException; > import org.apache.commons.math.exception.MathIllegalArgumentException; > -import org.apache.commons.math.exception.NumberIsTooLargeException; > +import org.apache.commons.math.exception.NonMonotonousSequenceException; > import org.apache.commons.math.exception.NotFiniteNumberException; > +import org.apache.commons.math.exception.NotPositiveException; > +import org.apache.commons.math.exception.NullArgumentException; > +import org.apache.commons.math.exception.NumberIsTooLargeException; > +import org.apache.commons.math.exception.util.Localizable; > +import org.apache.commons.math.exception.util.LocalizedFormats; > > /** > * Some useful additions to the built-in functions in {@link Math}. > @@ -91,6 +91,9 @@ public final class MathUtils { > 1307674368000l, 20922789888000l, 355687428096000l, > 6402373705728000l, 121645100408832000l, 2432902008176640000l }; > > + /** Factor used for splitting double numbers: n = 2^27 + 1. */ > + private static final int SPLIT_FACTOR = 0x8000001; > + > /** > * Private Constructor > */ > @@ -2332,4 +2335,277 @@ public final class MathUtils { > throw new NullArgumentException(); > } > } > + > + /** > + * Compute a linear combination accurately. > + * <p> > + * This method computes a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> 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 <a > + * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547";> > + * Accurate Sum and Dot Product</a> by Takeshi Ogita, > + * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. > Comput. > + * </p> > + * @param a1 first factor of the first term > + * @param b1 second factor of the first term > + * @param a2 first factor of the second term > + * @param b2 second factor of the second term > + * @return a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> > + * @see #linearCombination(double, double, double, double, double, > double) > + * @see #linearCombination(double, double, double, double, double, > double, double, double) > + */ > + public static double linearCombination(final double a1, final double b1, > + final double a2, final double b2) > { > + > + // the code below is split in many additions/subtractions that may > + // appear redundant. However, they shoud NOT be simplified, as they > + // do use IEEE753 floating point arithmetic rouding properties. > + // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same > as "a1" > + // The variables naming conventions are that xyzHigh contains the > most significant > + // bits of xyz and xyzLow contains its least significant bits. So > theoretically > + // xyz is the sum xyzHigh + xyzLow, but in many cases below, this > sum cannot > + // be represented in only one double precision number so we preserve > two numbers > + // to hold it as long as we can, combining the high and low order > bits together > + // only at the end, after cancellation may have occurred on high > order bits > + > + // split a1 and b1 as two 26 bits numbers > + final double ca1 = SPLIT_FACTOR * a1; > + final double a1High = ca1 - (ca1 - a1); > + final double a1Low = a1 - a1High; > + final double cb1 = SPLIT_FACTOR * b1; > + final double b1High = cb1 - (cb1 - b1); > + final double b1Low = b1 - b1High; > + > + // accurate multiplication a1 * b1 > + final double prod1High = a1 * b1; > + final double prod1Low = a1Low * b1Low - (((prod1High - a1High * > b1High) - a1Low * b1High) - a1High * b1Low); > + > + // split a2 and b2 as two 26 bits numbers > + final double ca2 = SPLIT_FACTOR * a2; > + final double a2High = ca2 - (ca2 - a2); > + final double a2Low = a2 - a2High; > + final double cb2 = SPLIT_FACTOR * b2; > + final double b2High = cb2 - (cb2 - b2); > + final double b2Low = b2 - b2High; > + > + // accurate multiplication a2 * b2 > + final double prod2High = a2 * b2; > + final double prod2Low = a2Low * b2Low - (((prod2High - a2High * > b2High) - a2Low * b2High) - a2High * b2Low); > + > + // accurate addition a1 * b1 + a2 * b2 > + final double s12High = prod1High + prod2High; > + final double s12Prime = s12High - prod2High; > + final double s12Low = (prod2High - (s12High - s12Prime)) + > (prod1High - s12Prime); > + > + // final rounding, s12 may have suffered many cancellations, we try > + // to recover some bits from the extra words we have saved up to now > + return s12High + (prod1Low + prod2Low + s12Low); > + > + } > + > + /** > + * Compute a linear combination accurately. > + * <p> > + * This method computes a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> > + * 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 <a > + * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547";> > + * Accurate Sum and Dot Product</a> by Takeshi Ogita, > + * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. > Comput. > + * </p> > + * @param a1 first factor of the first term > + * @param b1 second factor of the first term > + * @param a2 first factor of the second term > + * @param b2 second factor of the second term > + * @param a3 first factor of the third term > + * @param b3 second factor of the third term > + * @return a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> > + * @see #linearCombination(double, double, double, double) > + * @see #linearCombination(double, double, double, double, double, > double, double, double) > + */ > + public static double linearCombination(final double a1, final double b1, > + final double a2, final double b2, > + final double a3, final double b3) > { > + > + // the code below is split in many additions/subtractions that may > + // appear redundant. However, they shoud NOT be simplified, as they > + // do use IEEE753 floating point arithmetic rouding properties. > + // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same > as "a1" > + // The variables naming conventions are that xyzHigh contains the > most significant > + // bits of xyz and xyzLow contains its least significant bits. So > theoretically > + // xyz is the sum xyzHigh + xyzLow, but in many cases below, this > sum cannot > + // be represented in only one double precision number so we preserve > two numbers > + // to hold it as long as we can, combining the high and low order > bits together > + // only at the end, after cancellation may have occurred on high > order bits > + > + // split a1 and b1 as two 26 bits numbers > + final double ca1 = SPLIT_FACTOR * a1; > + final double a1High = ca1 - (ca1 - a1); > + final double a1Low = a1 - a1High; > + final double cb1 = SPLIT_FACTOR * b1; > + final double b1High = cb1 - (cb1 - b1); > + final double b1Low = b1 - b1High; > + > + // accurate multiplication a1 * b1 > + final double prod1High = a1 * b1; > + final double prod1Low = a1Low * b1Low - (((prod1High - a1High * > b1High) - a1Low * b1High) - a1High * b1Low); > + > + // split a2 and b2 as two 26 bits numbers > + final double ca2 = SPLIT_FACTOR * a2; > + final double a2High = ca2 - (ca2 - a2); > + final double a2Low = a2 - a2High; > + final double cb2 = SPLIT_FACTOR * b2; > + final double b2High = cb2 - (cb2 - b2); > + final double b2Low = b2 - b2High; > + > + // accurate multiplication a2 * b2 > + final double prod2High = a2 * b2; > + final double prod2Low = a2Low * b2Low - (((prod2High - a2High * > b2High) - a2Low * b2High) - a2High * b2Low); > + > + // split a3 and b3 as two 26 bits numbers > + final double ca3 = SPLIT_FACTOR * a3; > + final double a3High = ca3 - (ca3 - a3); > + final double a3Low = a3 - a3High; > + final double cb3 = SPLIT_FACTOR * b3; > + final double b3High = cb3 - (cb3 - b3); > + final double b3Low = b3 - b3High; > + > + // accurate multiplication a3 * b3 > + final double prod3High = a3 * b3; > + final double prod3Low = a3Low * b3Low - (((prod3High - a3High * > b3High) - a3Low * b3High) - a3High * b3Low); > + > + // accurate addition a1 * b1 + a2 * b2 > + final double s12High = prod1High + prod2High; > + final double s12Prime = s12High - prod2High; > + final double s12Low = (prod2High - (s12High - s12Prime)) + > (prod1High - s12Prime); > + > + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 > + final double s123High = s12High + prod3High; > + final double s123Prime = s123High - prod3High; > + final double s123Low = (prod3High - (s123High - s123Prime)) + > (s12High - s123Prime); > + > + // final rounding, s123 may have suffered many cancellations, we try > + // to recover some bits from the extra words we have saved up to now > + return s123High + (prod1Low + prod2Low + prod3Low + s12Low + > s123Low); > + > + } > + > + /** > + * Compute a linear combination accurately. > + * <p> > + * This method computes a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> > + > + * a<sub>4</sub>×b<sub>4</sub> > + * 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 <a > + * href="http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547";> > + * Accurate Sum and Dot Product</a> by Takeshi Ogita, > + * Siegfried M. Rump, and Shin'ichi Oishi published in SIAM J. Sci. > Comput. > + * </p> > + * @param a1 first factor of the first term > + * @param b1 second factor of the first term > + * @param a2 first factor of the second term > + * @param b2 second factor of the second term > + * @param a3 first factor of the third term > + * @param b3 second factor of the third term > + * @param a4 first factor of the third term > + * @param b4 second factor of the third term > + * @return a<sub>1</sub>×b<sub>1</sub> + > + * a<sub>2</sub>×b<sub>2</sub> + a<sub>3</sub>×b<sub>3</sub> > + > + * a<sub>4</sub>×b<sub>4</sub> > + * @see #linearCombination(double, double, double, double) > + * @see #linearCombination(double, double, double, double, double, > double) > + */ > + public 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) > { > + > + // the code below is split in many additions/subtractions that may > + // appear redundant. However, they shoud NOT be simplified, as they > + // do use IEEE753 floating point arithmetic rouding properties. > + // as an example, the expression "ca1 - (ca1 - a1)" is NOT the same > as "a1" > + // The variables naming conventions are that xyzHigh contains the > most significant > + // bits of xyz and xyzLow contains its least significant bits. So > theoretically > + // xyz is the sum xyzHigh + xyzLow, but in many cases below, this > sum cannot > + // be represented in only one double precision number so we preserve > two numbers > + // to hold it as long as we can, combining the high and low order > bits together > + // only at the end, after cancellation may have occurred on high > order bits > + > + // split a1 and b1 as two 26 bits numbers > + final double ca1 = SPLIT_FACTOR * a1; > + final double a1High = ca1 - (ca1 - a1); > + final double a1Low = a1 - a1High; > + final double cb1 = SPLIT_FACTOR * b1; > + final double b1High = cb1 - (cb1 - b1); > + final double b1Low = b1 - b1High; > + > + // accurate multiplication a1 * b1 > + final double prod1High = a1 * b1; > + final double prod1Low = a1Low * b1Low - (((prod1High - a1High * > b1High) - a1Low * b1High) - a1High * b1Low); > + > + // split a2 and b2 as two 26 bits numbers > + final double ca2 = SPLIT_FACTOR * a2; > + final double a2High = ca2 - (ca2 - a2); > + final double a2Low = a2 - a2High; > + final double cb2 = SPLIT_FACTOR * b2; > + final double b2High = cb2 - (cb2 - b2); > + final double b2Low = b2 - b2High; > + > + // accurate multiplication a2 * b2 > + final double prod2High = a2 * b2; > + final double prod2Low = a2Low * b2Low - (((prod2High - a2High * > b2High) - a2Low * b2High) - a2High * b2Low); > + > + // split a3 and b3 as two 26 bits numbers > + final double ca3 = SPLIT_FACTOR * a3; > + final double a3High = ca3 - (ca3 - a3); > + final double a3Low = a3 - a3High; > + final double cb3 = SPLIT_FACTOR * b3; > + final double b3High = cb3 - (cb3 - b3); > + final double b3Low = b3 - b3High; > + > + // accurate multiplication a3 * b3 > + final double prod3High = a3 * b3; > + final double prod3Low = a3Low * b3Low - (((prod3High - a3High * > b3High) - a3Low * b3High) - a3High * b3Low); > + > + // split a4 and b4 as two 26 bits numbers > + final double ca4 = SPLIT_FACTOR * a4; > + final double a4High = ca4 - (ca4 - a4); > + final double a4Low = a4 - a4High; > + final double cb4 = SPLIT_FACTOR * b4; > + final double b4High = cb4 - (cb4 - b4); > + final double b4Low = b4 - b4High; > + > + // accurate multiplication a4 * b4 > + final double prod4High = a4 * b4; > + final double prod4Low = a4Low * b4Low - (((prod4High - a4High * > b4High) - a4Low * b4High) - a4High * b4Low); > + > + // accurate addition a1 * b1 + a2 * b2 > + final double s12High = prod1High + prod2High; > + final double s12Prime = s12High - prod2High; > + final double s12Low = (prod2High - (s12High - s12Prime)) + > (prod1High - s12Prime); > + > + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 > + final double s123High = s12High + prod3High; > + final double s123Prime = s123High - prod3High; > + final double s123Low = (prod3High - (s123High - s123Prime)) + > (s12High - s123Prime); > + > + // accurate addition a1 * b1 + a2 * b2 + a3 * b3 + a4 * b4 > + final double s1234High = s123High + prod4High; > + final double s1234Prime = s1234High - prod4High; > + final double s1234Low = (prod4High - (s1234High - s1234Prime)) + > (s123High - s1234Prime); > + > + // final rounding, s1234 may have suffered many cancellations, we try > + // to recover some bits from the extra words we have saved up to now > + return s1234High + (prod1Low + prod2Low + prod3Low + prod4Low + > s12Low + s123Low + s1234Low); > + > + } > + > } > > Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml > URL: > http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1154250&r1=1154249&r2=1154250&view=diff > ============================================================================== > --- commons/proper/math/trunk/src/site/xdoc/changes.xml (original) > +++ commons/proper/math/trunk/src/site/xdoc/changes.xml Fri Aug 5 15:01:49 > 2011 > @@ -52,6 +52,13 @@ The <action> type attribute can be add,u > If the output is not quite correct, check for invisible trailing spaces! > --> > <release version="3.0" date="TBD" description="TBD"> > + <action dev="luc" type="add" > > + Added a few linearCombination utility methods in MathUtils to > compute accurately > + linear combinations a1.b1 + a2.b2 + ... + an.bn taking great care to > compensate > + for cancellation effects. This both improves and simplify several > methods in > + euclidean geometry classes, including linear constructors, dot > product and cross > + product. > + </action> > <action dev="psteitz" type="fix" issue="MATH-640"> > Fixed bugs in AbstractRandomGenerator nextInt() and nextLong() > default > implementations. Prior to the fix for this issue, these methods > > > --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
