File almostequal.h

File List > code_source > templat > include > TempLat > util > almostequal.h
#ifndef TEMPLAT_UTIL_ALMOSTEQUAL_H
#define TEMPLAT_UTIL_ALMOSTEQUAL_H

/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
   Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Wessel Valkenburg, Year: 2019

#include <limits> 
#include "TempLat/parallel/device.h"

namespace TempLat
{
  template <typename T1, typename T2, typename T3 = T1>
    requires(std::is_arithmetic_v<T1> && std::is_arithmetic_v<T2> && std::is_arithmetic_v<T3>)
  DEVICE_FUNCTION bool AlmostEqual(const T1 &a, const T2 &b,
                                   const T3 &epsilon = std::sqrt(std::numeric_limits<T3>::epsilon()))
  {
    if (std::isnan(a) || std::isnan(b)) return false;
    if (a == b) return true;

    // Relative error is meaningless once either side approaches zero, so below `epsilon` we fall back
    // to comparing the difference absolutely. This has to stay symmetric in a and b: comparing
    // "is the other one also below epsilon" instead would report two values straddling epsilon as
    // unequal however well they agree -- for epsilon = 1e3 * sqrt(FLT_EPSILON) = 0.345, the pair
    // (0.345338, 0.345192) took that branch and failed despite agreeing to 4e-4 relative.
    const bool aNearZero = std::abs(a) < epsilon;
    const bool bNearZero = std::abs(b) < epsilon;
    if (aNearZero || bNearZero) return (aNearZero && bNearZero) || std::abs(a - b) < epsilon;

    // Also test the absolute difference, to catch values that are effectively equal.
    return std::abs(a / b - 1) < epsilon || std::abs(a - b) < std::numeric_limits<T3>::epsilon() * 2;
  };

  template <typename T>
  DEVICE_FUNCTION bool AlmostEqual(const complex<T> &a, const complex<T> &b,
                                   const T epsilon = std::sqrt(std::numeric_limits<T>::epsilon()))
  {
    return AlmostEqual(a.real(), b.real(), epsilon) && AlmostEqual(a.imag(), b.imag(), epsilon);
  };

  template <typename T, size_t N>
  DEVICE_FUNCTION bool AlmostEqual(const std::array<T, N> &a, const std::array<T, N> &b,
                                   const T epsilon = std::sqrt(std::numeric_limits<T>::epsilon()))
  {
    bool result = true;
    for (ptrdiff_t i = 0; i < (ptrdiff_t)N; ++i) {
      result = result && AlmostEqual(a[i], b[i], epsilon);
    }
    return result;
  };

} // namespace TempLat

#endif