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