View Javadoc
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    *      http://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   */
17  package org.apache.commons.numbers.core;
18  
19  import java.util.function.DoubleConsumer;
20  import java.util.function.DoubleSupplier;
21  
22  /**
23   * Class providing accurate floating-point sums and linear combinations.
24   * This class uses techniques to mitigate round off errors resulting from
25   * standard floating-point operations, increasing the overall accuracy of
26   * results at the cost of a moderate increase in the number of computations.
27   * This functionality can be viewed as filling the gap between standard
28   * floating point operations (fast but prone to round off errors) and
29   * {@link java.math.BigDecimal} (perfectly accurate but slow).
30   *
31   * <p><strong>Usage</strong>
32   * <p>This class use a builder pattern in order to maximize the flexibility
33   * of the API. Typical use involves constructing an instance from one
34   * of the factory methods, adding any number of {@link #add(double) single value terms}
35   * and/or {@link #addProduct(double, double) products}, and then extracting the
36   * computed sum. Convenience methods exist for adding multiple values or products at once.
37   * The examples below demonstrate some simple use cases.
38   *
39   * <pre>
40   * // compute the sum a1 + a2 + a3 + a4
41   * Sum sum = Sum.of(a1);
42   *      .add(a2)
43   *      .add(a3)
44   *      .add(a4);
45   * double result = sum.getAsDouble();
46   *
47   * // same as above but using the varargs factory method
48   * double result = Sum.of(a1, a2, a3, a4).getAsDouble();
49   *
50   * // compute the dot product of two arrays of the same length, a and b
51   * Sum sum = Sum.create();
52   * for (int i = 0; i &lt; a.length; ++i) {
53   *      sum.addProduct(a[i], b[i]);
54   * }
55   * double result = sum.getAsDouble();
56   *
57   * // same as above but using a convenience factory method
58   * double result = Sum.ofProducts(a, b).getAsDouble();
59   * </pre>
60   *
61   * <p>It is worth noting that this class is designed to reduce floating point errors
62   * <em>across a sequence of operations</em> and not just a single add or multiply. The
63   * standard IEEE floating point operations already produce the most accurate results
64   * possible given two arguments and this class does not improve on them. Rather, it tracks
65   * the errors inherent with each operation and uses them to reduce the error of the overall
66   * result. Therefore, this class is only beneficial in cases involving 3 or more floating point
67   * operations. Code such as {@code Sum.of(a, b).getAsDouble()} and
68   * {@code Sum.create().addProduct(a, b).getAsDouble()} only adds overhead with no benefit.
69   *
70   * <p><strong>Implementation Notes</strong>
71   * <p>This class internally uses the <em>Sum2S</em> and <em>Dot2S</em> algorithms described in
72   * <a href="https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547">
73   * Accurate Sum and Dot Product</a> by Takeshi Ogita, Siegfried M. Rump,
74   * and Shin'ichi Oishi (<em>SIAM J. Sci. Comput</em>, 2005). These are compensated
75   * summation and multiplication algorithms chosen here for their good
76   * balance of precision and performance. Future releases may choose to use
77   * different algorithms.
78   *
79   * <p>Results follow the IEEE 754 rules for addition: For example, if any
80   * input value is {@link Double#NaN}, the result is {@link Double#NaN}.
81   *
82   * <p>Instances of this class are mutable and not safe for use by multiple
83   * threads.
84   */
85  public final class Sum
86      implements DoubleSupplier,
87                 DoubleConsumer {
88      /** Standard sum. */
89      private double sum;
90      /** Compensation value. */
91      private double comp;
92  
93      /**
94       * Constructs a new instance with the given initial value.
95       *
96       * @param initialValue Initial value.
97       */
98      private Sum(final double initialValue) {
99          sum = initialValue;
100         comp = 0d;
101     }
102 
103     /**
104      * Adds a single term to this sum.
105      *
106      * @param t Value to add.
107      * @return this instance.
108      */
109     public Sum add(final double t) {
110         final double newSum = sum + t;
111         comp += ExtendedPrecision.twoSumLow(sum, t, newSum);
112         sum = newSum;
113 
114         return this;
115     }
116 
117     /**
118      * Adds values from the given array to the sum.
119      *
120      * @param terms Terms to add.
121      * @return this instance.
122      */
123     public Sum add(final double... terms) {
124         for (double t : terms) {
125             add(t);
126         }
127 
128         return this;
129     }
130 
131     /**
132      * Adds the high-accuracy product \( a b \) to this sum.
133      *
134      * @param a Factor
135      * @param b Factor.
136      * @return this instance
137      */
138     public Sum addProduct(final double a,
139                           final double b) {
140         final double ab = a * b;
141         final double pLow = ExtendedPrecision.productLow(a, b, ab);
142 
143         final double newSum = sum + ab;
144         comp += ExtendedPrecision.twoSumLow(sum, ab, newSum) + pLow;
145         sum = newSum;
146 
147         return this;
148     }
149 
150     /**
151      * Adds \( \sum_i a_i b_i \) to this sum.
152      *
153      * @param a Factors.
154      * @param b Factors.
155      * @return this instance.
156      * @throws IllegalArgumentException if the arrays do not have the same length.
157      */
158     public Sum addProducts(final double[] a,
159                            final double[] b) {
160         final int len = a.length;
161         if (len != b.length) {
162             throw new IllegalArgumentException("Dimension mismatch: " +
163                                                a.length + " != " + b.length);
164         }
165 
166         for (int i = 0; i < len; ++i) {
167             addProduct(a[i], b[i]);
168         }
169 
170         return this;
171     }
172 
173     /**
174      * Adds another sum to this sum.
175      *
176      * @param other Sum to add.
177      * @return this instance.
178      */
179     public Sum add(final Sum other) {
180         // Pull both values first to ensure there are
181         // no issues when adding a sum to itself.
182         final double s = other.sum;
183         final double c = other.comp;
184 
185         return add(s).add(c);
186     }
187 
188     /**
189      * Adds a single term to this sum.
190      * This is equivalent to {@link #add(double)}.
191      *
192      * @param value Value to add.
193      *
194      * @see #add(double)
195      */
196     @Override
197     public void accept(final double value) {
198         add(value);
199     }
200 
201     /**
202      * Gets the sum value.
203      *
204      * @return the sum value.
205      */
206     @Override
207     public double getAsDouble() {
208         // High-precision value if it is finite, standard IEEE754 result otherwise.
209         final double hpsum = sum + comp;
210         return Double.isFinite(hpsum) ?
211                 hpsum :
212                 sum;
213     }
214 
215     /**
216      * Creates a new instance with an initial value of zero.
217      *
218      * @return a new instance.
219      */
220     public static Sum create() {
221         return new Sum(0d);
222     }
223 
224     /**
225      * Creates an instance initialized to the given value.
226      *
227      * @param a Initial value.
228      * @return a new instance.
229      */
230     public static Sum of(final double a) {
231         return new Sum(a);
232     }
233 
234     /**
235      * Creates an instance containing the sum of the given values.
236      *
237      * @param values Values to add.
238      * @return a new instance.
239      */
240     public static Sum of(final double... values) {
241         return create().add(values);
242     }
243 
244     /**
245      * Creates a new instance containing \( \sum_i a_i b_i \).
246      *
247      * @param a Factors.
248      * @param b Factors.
249      * @return a new instance.
250      */
251     public static Sum ofProducts(final double[] a,
252                                  final double[] b) {
253         return create().addProducts(a, b);
254     }
255 }