From 1b8ded4702a9e4352e13235d3f646299750bb2b1 Mon Sep 17 00:00:00 2001 From: Thomas Fischer <thomas.fischer@ufz.de> Date: Thu, 25 Jul 2024 14:31:58 +0200 Subject: [PATCH] [MaL] Implement KahanSum algorithm --- MathLib/KahanSum.h | 49 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 MathLib/KahanSum.h diff --git a/MathLib/KahanSum.h b/MathLib/KahanSum.h new file mode 100644 index 00000000000..0a360f355dd --- /dev/null +++ b/MathLib/KahanSum.h @@ -0,0 +1,49 @@ +/** + * \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 <iomanip> +#include <ostream> + +namespace MathLib +{ + +class KahanSum +{ +public: + explicit constexpr KahanSum(double const value = 0) : value_(value) {} + + 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 -- GitLab