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 < 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 }