/** * \file * \copyright * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org) * Distributed under a Modified BSD License. * See accompanying file LICENSE.txt or * http://www.opengeosys.org/project/license */ #pragma once #include <spdlog/fmt/bundled/ostream.h> #include <iomanip> #include <ostream> #include <range/v3/range/concepts.hpp> #include <range/v3/range/traits.hpp> namespace MathLib { class KahanSum; } template <> struct fmt::formatter<MathLib::KahanSum> : fmt::ostream_formatter { }; namespace MathLib { class KahanSum { public: explicit constexpr KahanSum(double const value = 0) : value_(value) {} explicit constexpr KahanSum(ranges::range auto const& range) : value_(0) { for (auto const v : range) { *this += v; } } constexpr KahanSum operator+(double const increment) const { KahanSum result = *this; return result += increment; } constexpr KahanSum operator-(double const increment) const { KahanSum result = *this; return result += -increment; } constexpr KahanSum& operator-=(double const increment) { return *this += -increment; } constexpr KahanSum& operator+=(double const increment) { double const y = increment - correction_; double const t = value_ + y; correction_ = (t - value_) - y; value_ = t; return *this; } constexpr double value() const { return value_; } constexpr double operator()() const { return value_; } friend inline std::ostream& operator<<(std::ostream& os, KahanSum const& x) { auto const precision = os.precision(); return os << std::setprecision( std::numeric_limits<double>::max_digits10) << x.value() << " (± " << x.correction_ << ')' << std::setprecision(precision); } private: double value_; double correction_ = 0.; }; } // namespace MathLib