00001 /** 00002 * \file Accumulator.hpp 00003 * \brief Header for GeographicLib::Accumulator class 00004 * 00005 * Copyright (c) Charles Karney (2010-2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_ACCUMULATOR_HPP) 00011 #define GEOGRAPHICLIB_ACCUMULATOR_HPP 1 00012 00013 #include <GeographicLib/Constants.hpp> 00014 00015 namespace GeographicLib { 00016 00017 /** 00018 * \brief An accumulator for sums 00019 * 00020 * This allow many numbers of floating point type \e T to be added together 00021 * with twice the normal precision. Thus if \e T is double, the effective 00022 * precision of the sum is 106 bits or about 32 decimal places. 00023 * 00024 * The implementation follows J. R. Shewchuk, 00025 * <a href="http://dx.doi.org/10.1007/PL00009321"> Adaptive Precision 00026 * Floating-Point Arithmetic and Fast Robust Geometric Predicates</a>, 00027 * Discrete & Computational Geometry 18(3) 305--363 (1997). 00028 * 00029 * Approximate timings (summing a vector<double>) 00030 * - double: 2ns 00031 * - Accumulator<double>: 23ns 00032 * 00033 * In the documentation of the member functions, \e sum stands for the value 00034 * currently held in the accumulator. 00035 * 00036 * Example of use: 00037 * \include example-Accumulator.cpp 00038 **********************************************************************/ 00039 template<typename T = Math::real> 00040 class GEOGRAPHICLIB_EXPORT Accumulator { 00041 private: 00042 // _s + _t accumulators for the sum. 00043 T _s, _t; 00044 // Same as Math::sum, but requires abs(u) >= abs(v). This isn't currently 00045 // used. 00046 static inline T fastsum(T u, T v, T& t) { 00047 GEOGRAPHICLIB_VOLATILE T s = u + v; 00048 GEOGRAPHICLIB_VOLATILE T vp = s - u; 00049 t = v - vp; 00050 return s; 00051 } 00052 void Add(T y) { 00053 // Here's Shewchuk's solution... 00054 T u; // hold exact sum as [s, t, u] 00055 y = Math::sum(y, _t, u); // Accumulate starting at least significant end 00056 _s = Math::sum(y, _s, _t); 00057 // Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u) 00058 // exactly with s, t, u non-adjacent and in decreasing order (except for 00059 // possible zeros). The following code tries to normalize the result. 00060 // Ideally, we want _s = round(s+t+u) and _u = round(s+t+u - _s). The 00061 // following does an approximate job (and maintains the decreasing 00062 // non-adjacent property). Here are two "failures" using 3-bit floats: 00063 // 00064 // Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp 00065 // [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3 00066 // 00067 // Case 2: _s+_t is not as close to s+t+u as it shold be 00068 // [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1) 00069 // should be [80, -7] = 73 (exact) 00070 // 00071 // "Fixing" these problems is probably not worth the expense. The 00072 // representation inevitably leads to small errors in the accumulated 00073 // values. The additional errors illustrated here amount to 1 ulp of the 00074 // less significant word during each addition to the Accumulator and an 00075 // additional possible error of 1 ulp in the reported sum. 00076 // 00077 // Incidentally, the "ideal" representation described above is not 00078 // canonical, because _s = round(_s + _t) may not be true. For example, 00079 // with 3-bit floats: 00080 // 00081 // [128, 16] + 1 -> [160, -16] -- 160 = round(145). 00082 // But [160, 0] - 16 -> [128, 16] -- 128 = round(144). 00083 // 00084 if (_s == 0) // This implies t == 0, 00085 _s = u; // so result is u 00086 else 00087 _t += u; // otherwise just accumulate u to t. 00088 } 00089 T Sum(T y) const { 00090 Accumulator a(*this); 00091 a.Add(y); 00092 return a._s; 00093 } 00094 public: 00095 /** 00096 * Construct from a \e T. This is not declared explicit, so that you can 00097 * write <code>Accumulator<double> a = 5;</code>. 00098 * 00099 * @param[in] y set \e sum = \e y. 00100 **********************************************************************/ 00101 Accumulator(T y = T(0)) : _s(y), _t(0) { 00102 GEOGRAPHICLIB_STATIC_ASSERT(!std::numeric_limits<T>::is_integer, 00103 "Accumulator type is not floating point"); 00104 } 00105 /** 00106 * Set the accumulator to a number. 00107 * 00108 * @param[in] y set \e sum = \e y. 00109 **********************************************************************/ 00110 Accumulator& operator=(T y) { _s = y; _t = 0; return *this; } 00111 /** 00112 * Return the value held in the accumulator. 00113 * 00114 * @return \e sum. 00115 **********************************************************************/ 00116 T operator()() const { return _s; } 00117 /** 00118 * Return the result of adding a number to \e sum (but don't change \e sum). 00119 * 00120 * @param[in] y the number to be added to the sum. 00121 * @return \e sum + \e y. 00122 **********************************************************************/ 00123 T operator()(T y) const { return Sum(y); } 00124 /** 00125 * Add a number to the accumulator. 00126 * 00127 * @param[in] y set \e sum += \e y. 00128 **********************************************************************/ 00129 Accumulator& operator+=(T y) { Add(y); return *this; } 00130 /** 00131 * Subtract a number from the accumulator. 00132 * 00133 * @param[in] y set \e sum -= \e y. 00134 **********************************************************************/ 00135 Accumulator& operator-=(T y) { Add(-y); return *this; } 00136 /** 00137 * Multiply accumulator by an integer. To avoid loss of accuracy, use only 00138 * integers such that \e n × \e T is exactly representable as a \e T 00139 * (i.e., ± powers of two). Use \e n = −1 to negate \e sum. 00140 * 00141 * @param[in] n set \e sum *= \e n. 00142 **********************************************************************/ 00143 Accumulator& operator*=(int n) { _s *= n; _t *= n; return *this; } 00144 /** 00145 * Test equality of an Accumulator with a number. 00146 **********************************************************************/ 00147 bool operator==(T y) const { return _s == y; } 00148 /** 00149 * Test inequality of an Accumulator with a number. 00150 **********************************************************************/ 00151 bool operator!=(T y) const { return _s != y; } 00152 /** 00153 * Less operator on an Accumulator and a number. 00154 **********************************************************************/ 00155 bool operator<(T y) const { return _s < y; } 00156 /** 00157 * Less or equal operator on an Accumulator and a number. 00158 **********************************************************************/ 00159 bool operator<=(T y) const { return _s <= y; } 00160 /** 00161 * Greater operator on an Accumulator and a number. 00162 **********************************************************************/ 00163 bool operator>(T y) const { return _s > y; } 00164 /** 00165 * Greater or equal operator on an Accumulator and a number. 00166 **********************************************************************/ 00167 bool operator>=(T y) const { return _s >= y; } 00168 }; 00169 00170 } // namespace GeographicLib 00171 00172 #endif // GEOGRAPHICLIB_ACCUMULATOR_HPP