numerical_utils.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <iterator>
#include <limits>
#include <stdexcept>
 
namespace statcpp {
 
constexpr double epsilon = std::numeric_limits<double>::epsilon();
 
constexpr double default_rel_tol = 1e-9;
 
constexpr double default_abs_tol = 1e-12;
 
inline bool approx_equal(double a, double b,
                         double rel_tol = default_rel_tol,
                         double abs_tol = default_abs_tol)
{
    // Handle exact equality (including infinities)
    if (a == b) {
        return true;
    }
 
    // Handle NaN
    if (std::isnan(a) || std::isnan(b)) {
        return false;
    }
 
    // Compute absolute difference
    double diff = std::abs(a - b);
 
    // Check absolute tolerance
    if (diff <= abs_tol) {
        return true;
    }
 
    // Check relative tolerance
    double max_abs = std::max(std::abs(a), std::abs(b));
    return diff <= rel_tol * max_abs;
}
 
inline bool is_zero(double x, double tol = default_abs_tol)
{
    return std::abs(x) <= tol;
}
 
inline bool is_finite(double x)
{
    return std::isfinite(x);
}
 
template <typename Iterator>
bool all_finite(Iterator first, Iterator last)
{
    for (auto it = first; it != last; ++it) {
        if (!std::isfinite(static_cast<double>(*it))) {
            return false;
        }
    }
    return true;
}
 
inline bool has_converged_abs(double x_new, double x_old,
                               double tol = 1e-6)
{
    return std::abs(x_new - x_old) <= tol;
}
 
inline bool has_converged_rel(double x_new, double x_old,
                               double tol = 1e-6)
{
    double denominator = std::max(std::abs(x_old), epsilon);
    return std::abs(x_new - x_old) / denominator <= tol;
}
 
inline bool has_converged(double x_new, double x_old,
                          double abs_tol = 1e-8,
                          double rel_tol = 1e-6)
{
    double threshold = abs_tol + rel_tol * std::abs(x_old);
    return std::abs(x_new - x_old) <= threshold;
}
 
inline double log1p_safe(double x)
{
    return std::log1p(x);
}
 
inline double expm1_safe(double x)
{
    return std::expm1(x);
}
 
inline double clamp(double x, double min_val, double max_val)
{
    if (min_val > max_val) {
        throw std::invalid_argument("statcpp::clamp: min_val must be <= max_val");
    }
    return std::max(min_val, std::min(x, max_val));
}
 
inline bool in_range(double x, double min_val, double max_val)
{
    return x >= min_val && x <= max_val;
}
 
inline double relative_error(double x, double x_ref)
{
    double denominator = std::max(std::abs(x_ref), epsilon);
    return std::abs(x - x_ref) / denominator;
}
 
inline double safe_divide(double numerator, double denominator,
                          double default_value = std::numeric_limits<double>::quiet_NaN())
{
    if (is_zero(denominator)) {
        return default_value;
    }
    return numerator / denominator;
}
 
template <typename Iterator>
double kahan_sum(Iterator first, Iterator last)
{
    double sum = 0.0;
    double compensation = 0.0;  // Running compensation for lost low-order bits
 
    for (auto it = first; it != last; ++it) {
        double value = static_cast<double>(*it);
        double y = value - compensation;
        double t = sum + y;
        compensation = (t - sum) - y;
        sum = t;
    }
 
    return sum;
}
 
template <typename Iterator, typename Proj>
double kahan_sum(Iterator first, Iterator last, Proj proj)
{
    double sum = 0.0;
    double compensation = 0.0;
 
    for (auto it = first; it != last; ++it) {
        double value = static_cast<double>(proj(*it));
        double y = value - compensation;
        double t = sum + y;
        compensation = (t - sum) - y;
        sum = t;
    }
 
    return sum;
}
 
template <typename Iterator1, typename Iterator2>
bool approx_equal_range(Iterator1 first1, Iterator1 last1,
                        Iterator2 first2, Iterator2 last2,
                        double rel_tol = default_rel_tol,
                        double abs_tol = default_abs_tol)
{
    auto n1 = std::distance(first1, last1);
    auto n2 = std::distance(first2, last2);
 
    if (n1 != n2) {
        return false;
    }
 
    auto it1 = first1;
    auto it2 = first2;
 
    while (it1 != last1) {
        double v1 = static_cast<double>(*it1);
        double v2 = static_cast<double>(*it2);
 
        if (!approx_equal(v1, v2, rel_tol, abs_tol)) {
            return false;
        }
 
        ++it1;
        ++it2;
    }
 
    return true;
}
 
} // namespace statcpp