LSSTApplications  20.0.0
LSSTDataManagementBasePackage
SafeSum.h
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * Developed for the LSST Data Management System.
4  * This product includes software developed by the LSST Project
5  * (https://www.lsst.org).
6  * See the COPYRIGHT file at the top-level directory of this distribution
7  * for details of code ownership.
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <https://www.gnu.org/licenses/>.
21  */
22 #ifndef LSST_AFW_MATH_POLYNOMIALS_SafeSum_h_INCLUDED
23 #define LSST_AFW_MATH_POLYNOMIALS_SafeSum_h_INCLUDED
24 
25 #include <cmath>
26 
27 namespace lsst { namespace geom { namespace polynomials {
28 
32 enum class SumMode {
34  FAST,
36  SAFE
37 };
38 
61 template <typename T>
62 class SafeSum {
63 public:
64 
65  explicit SafeSum(T initial=static_cast<T>(0)) noexcept :
66  _sum(initial),
67  _correction(static_cast<T>(0))
68  {}
69 
70  SafeSum & operator=(T value) noexcept {
71  _sum = value;
72  _correction = static_cast<T>(0);
73  return *this;
74  }
75 
76  SafeSum & operator+=(T value) noexcept {
77  T t = _sum + value;
78  // update _correction to account for lost low-order bits of the
79  // greater of _sum and value
80  if (std::abs(_sum) >= std::abs(value)) {
81  _correction += (_sum - t) + value;
82  } else {
83  _correction += (value - t) + _sum;
84  }
85  _sum = t;
86  return *this;
87  }
88 
89  SafeSum & operator-=(T value) noexcept {
90  return operator+=(-value);
91  }
92 
93  explicit operator T() const noexcept {
94  return _sum + _correction;
95  }
96 
97 private:
98  T _sum;
99  T _correction;
100 };
101 
102 }}} // namespace lsst::geom::polynomials
103 
104 #endif // !LSST_AFW_MATH_POLYNOMIALS_SafeSum_h_INCLUDED
lsst::sphgeom::abs
Angle abs(Angle const &a)
Definition: Angle.h:106
lsst::geom::polynomials::SafeSum::SafeSum
SafeSum(T initial=static_cast< T >(0)) noexcept
Definition: SafeSum.h:65
lsst::geom::polynomials::SumMode
SumMode
Enum used to control how to sum polynomial terms.
Definition: SafeSum.h:32
lsst::geom::polynomials::SafeSum::operator+=
SafeSum & operator+=(T value) noexcept
Definition: SafeSum.h:76
lsst::geom::polynomials::SumMode::SAFE
@ SAFE
Compensated summation using SafeSum. Involves ~4x as many floating point operations.
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::geom
Definition: geomOperators.dox:4
lsst::geom::polynomials::SafeSum
A numerically stable summation algorithm for floating-point numbers.
Definition: SafeSum.h:62
lsst::geom::polynomials::SumMode::FAST
@ FAST
Summation using regular floating-point addition.
lsst::geom::polynomials::SafeSum::operator-=
SafeSum & operator-=(T value) noexcept
Definition: SafeSum.h:89
lsst::geom::polynomials::SafeSum::operator=
SafeSum & operator=(T value) noexcept
Definition: SafeSum.h:70