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://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.2.1547"> 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 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}