001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.core;
018
019import java.util.function.DoubleConsumer;
020import java.util.function.DoubleSupplier;
021
022/**
023 * Class providing accurate floating-point sums and linear combinations.
024 * This class uses techniques to mitigate round off errors resulting from
025 * standard floating-point operations, increasing the overall accuracy of
026 * results at the cost of a moderate increase in the number of computations.
027 * This functionality can be viewed as filling the gap between standard
028 * floating point operations (fast but prone to round off errors) and
029 * {@link java.math.BigDecimal} (perfectly accurate but slow).
030 *
031 * <p><strong>Usage</strong>
032 * <p>This class use a builder pattern in order to maximize the flexibility
033 * of the API. Typical use involves constructing an instance from one
034 * of the factory methods, adding any number of {@link #add(double) single value terms}
035 * and/or {@link #addProduct(double, double) products}, and then extracting the
036 * computed sum. Convenience methods exist for adding multiple values or products at once.
037 * The examples below demonstrate some simple use cases.
038 *
039 * <pre>
040 * // compute the sum a1 + a2 + a3 + a4
041 * Sum sum = Sum.of(a1);
042 *      .add(a2)
043 *      .add(a3)
044 *      .add(a4);
045 * double result = sum.getAsDouble();
046 *
047 * // same as above but using the varargs factory method
048 * double result = Sum.of(a1, a2, a3, a4).getAsDouble();
049 *
050 * // compute the dot product of two arrays of the same length, a and b
051 * Sum sum = Sum.create();
052 * for (int i = 0; i &lt; a.length; ++i) {
053 *      sum.addProduct(a[i], b[i]);
054 * }
055 * double result = sum.getAsDouble();
056 *
057 * // same as above but using a convenience factory method
058 * double result = Sum.ofProducts(a, b).getAsDouble();
059 * </pre>
060 *
061 * <p>It is worth noting that this class is designed to reduce floating point errors
062 * <em>across a sequence of operations</em> and not just a single add or multiply. The
063 * standard IEEE floating point operations already produce the most accurate results
064 * possible given two arguments and this class does not improve on them. Rather, it tracks
065 * the errors inherent with each operation and uses them to reduce the error of the overall
066 * result. Therefore, this class is only beneficial in cases involving 3 or more floating point
067 * operations. Code such as {@code Sum.of(a, b).getAsDouble()} and
068 * {@code Sum.create().addProduct(a, b).getAsDouble()} only adds overhead with no benefit.
069 *
070 * <p><strong>Implementation Notes</strong>
071 * <p>This class internally uses the <em>Sum2S</em> and <em>Dot2S</em> algorithms described in
072 * <a href="https://doi.org/10.1137/030601818">
073 * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
074 * and Shin'ichi Oishi (<em>SIAM J. Sci. Comput</em>, 2005). These are compensated
075 * summation and multiplication algorithms chosen here for their good
076 * balance of precision and performance. Future releases may choose to use
077 * different algorithms.
078 *
079 * <p>Results follow the IEEE 754 rules for addition: For example, if any
080 * input value is {@link Double#NaN}, the result is {@link Double#NaN}.
081 *
082 * <p>Instances of this class are mutable and not safe for use by multiple
083 * threads.
084 */
085public final class Sum
086    implements DoubleSupplier,
087               DoubleConsumer {
088    /** Standard sum. */
089    private double sum;
090    /** Compensation value. */
091    private double comp;
092
093    /**
094     * Constructs a new instance with the given initial value.
095     *
096     * @param initialValue Initial value.
097     */
098    private Sum(final double initialValue) {
099        sum = initialValue;
100    }
101
102    /**
103     * Adds a single term to this sum.
104     *
105     * @param t Value to add.
106     * @return this instance.
107     */
108    public Sum add(final double t) {
109        final double newSum = sum + t;
110        comp += DD.twoSumLow(sum, t, newSum);
111        sum = newSum;
112
113        return this;
114    }
115
116    /**
117     * Adds values from the given array to the sum.
118     *
119     * @param terms Terms to add.
120     * @return this instance.
121     */
122    public Sum add(final double... terms) {
123        for (final double t : terms) {
124            add(t);
125        }
126
127        return this;
128    }
129
130    /**
131     * Adds the high-accuracy product \( a b \) to this sum.
132     *
133     * @param a Factor
134     * @param b Factor.
135     * @return this instance
136     */
137    public Sum addProduct(final double a,
138                          final double b) {
139        final double ab = a * b;
140        final double pLow = ExtendedPrecision.productLow(a, b, ab);
141
142        final double newSum = sum + ab;
143        comp += DD.twoSumLow(sum, ab, newSum) + pLow;
144        sum = newSum;
145
146        return this;
147    }
148
149    /**
150     * Adds \( \sum_i a_i b_i \) to this sum.
151     *
152     * @param a Factors.
153     * @param b Factors.
154     * @return this instance.
155     * @throws IllegalArgumentException if the arrays do not have the same length.
156     */
157    public Sum addProducts(final double[] a,
158                           final double[] b) {
159        final int len = a.length;
160        if (len != b.length) {
161            throw new IllegalArgumentException("Dimension mismatch: " +
162                                               a.length + " != " + b.length);
163        }
164
165        for (int i = 0; i < len; ++i) {
166            addProduct(a[i], b[i]);
167        }
168
169        return this;
170    }
171
172    /**
173     * Adds another sum to this sum.
174     *
175     * @param other Sum to add.
176     * @return this instance.
177     */
178    public Sum add(final Sum other) {
179        return add(other.sum, other.comp);
180    }
181
182    /**
183     * Subtracts another sum from this sum.
184     *
185     * @param other Sum to subtract.
186     * @return this instance.
187     * @since 1.2
188     */
189    public Sum subtract(final Sum other) {
190        return add(-other.sum, -other.comp);
191    }
192
193    /**
194     * Adds the sum and compensation to this sum.
195     *
196     * <p>This is a utility method to extract both values from a sum
197     * before addition to ensure there are no issues when adding a sum
198     * to itself.
199     *
200     * <p>This method sums the values in double-double precision.
201     * This enforces symmetry when combining sums and maximises the available precision.
202     *
203     * @param s Sum.
204     * @param c Compensation.
205     * @return this instance.
206     */
207    private Sum add(double s, double c) {
208        // Re-normalise the sums and combine in double-double precision.
209        // Note: The conversion to a DD is lossless.
210        final DD result = DD.ofSum(sum, comp).add(DD.ofSum(s, c));
211        if (result.isFinite()) {
212            sum = result.hi();
213            comp = result.lo();
214        } else {
215            // compensation can be NaN from accumulating one or more same-signed infinite values.
216            // Do not pollute the regular IEEE754 sum with a spurious NaN.
217            add(s);
218            if (!Double.isNaN(c)) {
219                add(c);
220            }
221        }
222        return this;
223    }
224
225    /**
226     * Adds a single term to this sum.
227     * This is equivalent to {@link #add(double)}.
228     *
229     * @param value Value to add.
230     *
231     * @see #add(double)
232     */
233    @Override
234    public void accept(final double value) {
235        add(value);
236    }
237
238    /**
239     * Gets the sum value.
240     *
241     * @return the sum value.
242     */
243    @Override
244    public double getAsDouble() {
245        // High-precision value if it is finite, standard IEEE754 result otherwise.
246        final double hpsum = sum + comp;
247        return Double.isFinite(hpsum) ?
248                hpsum :
249                sum;
250    }
251
252    /**
253     * Creates a new instance with an initial value of zero.
254     *
255     * @return a new instance.
256     */
257    public static Sum create() {
258        return new Sum(0d);
259    }
260
261    /**
262     * Creates an instance initialized to the given value.
263     *
264     * @param a Initial value.
265     * @return a new instance.
266     */
267    public static Sum of(final double a) {
268        return new Sum(a);
269    }
270
271    /**
272     * Creates an instance containing the sum of the given values.
273     *
274     * @param values Values to add.
275     * @return a new instance.
276     */
277    public static Sum of(final double... values) {
278        return create().add(values);
279    }
280
281    /**
282     * Creates a new instance containing \( \sum_i a_i b_i \).
283     *
284     * @param a Factors.
285     * @param b Factors.
286     * @return a new instance.
287     */
288    public static Sum ofProducts(final double[] a,
289                                 final double[] b) {
290        return create().addProducts(a, b);
291    }
292}