Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-11-01 10:04:16

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file corecel/math/SoftEqual.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008 
0009 #include <cmath>
0010 
0011 #include "corecel/Assert.hh"
0012 #include "corecel/Macros.hh"
0013 #include "corecel/Types.hh"
0014 
0015 #include "detail/SoftEqualTraits.hh"
0016 
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Square root of the soft equivalence tolerance for Celeritas.
0022  *
0023  * This tolerance is needed for operations where the accuracy is limited by the
0024  * square root of machine precision.
0025  *
0026  * \todo Move orange tolerance and related operations into corecel/math
0027  * alongside this, revisit ArrayUtils soft comparisons.
0028  */
0029 template<class T = ::celeritas::real_type>
0030 constexpr inline T sqrt_tol = detail::SoftEqualTraits<T>::sqrt_prec();
0031 
0032 //---------------------------------------------------------------------------//
0033 /*!
0034  * Functor for noninfinite floating point equality.
0035  *
0036  * This function-like class considers an \em absolute tolerance for values near
0037  * zero, and a \em relative tolerance for values far from zero. It correctly
0038  * returns "false" if either value being compared is NaN.  The call operator is
0039  * \em commutative: \c eq(a,b) should always give the same as \c eq(b,a).
0040  *
0041  * The actual comparison is: \f[
0042  |a - b| < \max(\epsilon_r \max(|a|, |b|), \epsilon_a)
0043  \f]
0044  *
0045  * \note The edge case where both values are infinite (with the same sign)
0046  * returns \em false for equality, which could be considered reasonable because
0047  * relative error is meaningless. To explicitly allow infinities to compare
0048  * equal, you must test separately, e.g., `a == b || soft_eq(a, b)` using \c
0049  * celeritas::EqualOr.
0050  */
0051 template<class RealType = ::celeritas::real_type>
0052 class SoftEqual
0053 {
0054   public:
0055     //!@{
0056     //! \name Type aliases
0057     using value_type = RealType;
0058     //!@}
0059 
0060   public:
0061     // Construct with default relative/absolute precision
0062     CELER_CONSTEXPR_FUNCTION SoftEqual();
0063 
0064     // Construct with default absolute precision
0065     explicit CELER_CONSTEXPR_FUNCTION SoftEqual(value_type rel);
0066 
0067     // Construct with both relative and absolute precision
0068     CELER_CONSTEXPR_FUNCTION SoftEqual(value_type rel, value_type abs);
0069 
0070     //// COMPARISON ////
0071 
0072     // Compare two values (implicitly casting arguments)
0073     bool CELER_FUNCTION operator()(value_type a, value_type b) const;
0074 
0075     //// ACCESSORS ////
0076 
0077     //! Relative allowable error
0078     CELER_CONSTEXPR_FUNCTION value_type rel() const { return rel_; }
0079 
0080     //! Absolute tolerance
0081     CELER_CONSTEXPR_FUNCTION value_type abs() const { return abs_; }
0082 
0083   private:
0084     value_type rel_;
0085     value_type abs_;
0086 
0087     using SETraits = detail::SoftEqualTraits<value_type>;
0088 };
0089 
0090 //---------------------------------------------------------------------------//
0091 /*!
0092  * Compare for equality before checking with the given functor.
0093  *
0094  * This CRTP class allows \c SoftEqual to work for infinities.
0095  */
0096 template<class F>
0097 class EqualOr : public F
0098 {
0099   public:
0100     //! Forward arguments to parent class
0101     template<class... C>
0102     CELER_FUNCTION EqualOr(C&&... args) : F{std::forward<C>(args)...}
0103     {
0104     }
0105 
0106     //! Forward arguments to comparison operator after comparing
0107     template<class T, class U>
0108     bool CELER_FUNCTION operator()(T a, U b) const
0109     {
0110         return a == b || static_cast<F const&>(*this)(a, b);
0111     }
0112 };
0113 
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Functor for comparing values with an absolute precision.
0117  */
0118 template<class RealType = ::celeritas::real_type>
0119 class SoftClose
0120 {
0121   public:
0122     //!@{
0123     //! \name Type aliases
0124     using value_type = RealType;
0125     //!@}
0126 
0127   public:
0128     // Construct with default absolute precision
0129     CELER_CONSTEXPR_FUNCTION SoftClose();
0130 
0131     // Construct with absolute precision
0132     CELER_FUNCTION SoftClose(value_type abs);
0133 
0134     //// COMPARISON ////
0135 
0136     // Compare two values (implicitly casting arguments)
0137     bool CELER_FUNCTION operator()(value_type a, value_type b) const;
0138 
0139     //// ACCESSORS ////
0140 
0141     //! Absolute tolerance
0142     CELER_CONSTEXPR_FUNCTION value_type abs() const { return abs_; }
0143 
0144   private:
0145     value_type abs_;
0146 
0147     using SETraits = detail::SoftEqualTraits<value_type>;
0148 };
0149 
0150 //---------------------------------------------------------------------------//
0151 /*!
0152  * Functor for comparing values to zero with an absolute precision.
0153  */
0154 template<class RealType = ::celeritas::real_type>
0155 class SoftZero
0156 {
0157   public:
0158     //!@{
0159     //! \name Type aliases
0160     using argument_type = RealType;
0161     using value_type = RealType;
0162     //!@}
0163 
0164   public:
0165     // Construct with default absolute precision
0166     CELER_CONSTEXPR_FUNCTION SoftZero();
0167 
0168     // Construct with absolute precision
0169     explicit CELER_FUNCTION SoftZero(value_type abs);
0170 
0171     //// COMPARISON ////
0172 
0173     // Compare the given value to zero
0174     inline CELER_FUNCTION bool operator()(value_type actual) const;
0175 
0176     //// ACCESSORS ////
0177 
0178     //! Absolute tolerance
0179     CELER_CONSTEXPR_FUNCTION value_type abs() const { return abs_; }
0180 
0181   private:
0182     value_type abs_;
0183 
0184     using SETraits = detail::SoftEqualTraits<value_type>;
0185 };
0186 
0187 //---------------------------------------------------------------------------//
0188 // TEMPLATE DEDUCTION GUIDES
0189 //---------------------------------------------------------------------------//
0190 template<class T>
0191 CELER_CTAD_FUNCTION SoftEqual(T) -> SoftEqual<T>;
0192 template<class T>
0193 CELER_CTAD_FUNCTION SoftEqual(T, T) -> SoftEqual<T>;
0194 template<class F>
0195 CELER_CTAD_FUNCTION EqualOr(F&&) -> EqualOr<F>;
0196 
0197 //---------------------------------------------------------------------------//
0198 // INLINE DEFINITIONS
0199 //---------------------------------------------------------------------------//
0200 /*!
0201  * Construct with default relative/absolute precision.
0202  */
0203 template<class RealType>
0204 CELER_CONSTEXPR_FUNCTION SoftEqual<RealType>::SoftEqual()
0205     : rel_{SETraits::rel_prec()}, abs_{SETraits::abs_thresh()}
0206 {
0207 }
0208 
0209 //---------------------------------------------------------------------------//
0210 /*!
0211  * Construct with scaled absolute precision.
0212  */
0213 template<class RealType>
0214 CELER_CONSTEXPR_FUNCTION SoftEqual<RealType>::SoftEqual(value_type rel)
0215     : SoftEqual(rel, rel * (SETraits::abs_thresh() / SETraits::rel_prec()))
0216 {
0217     CELER_EXPECT(rel > 0);
0218 }
0219 
0220 //---------------------------------------------------------------------------//
0221 /*!
0222  * Construct with both relative and absolute precision.
0223  *
0224  * \param rel tolerance of relative error (default 1.0e-12 for doubles)
0225  *
0226  * \param abs threshold for absolute error when comparing small quantities
0227  *           (default 1.0e-14 for doubles)
0228  */
0229 template<class RealType>
0230 CELER_CONSTEXPR_FUNCTION SoftEqual<RealType>::SoftEqual(value_type rel,
0231                                                         value_type abs)
0232     : rel_{rel}, abs_{abs}
0233 {
0234     CELER_EXPECT(rel > 0);
0235     CELER_EXPECT(abs > 0);
0236 }
0237 
0238 //---------------------------------------------------------------------------//
0239 /*!
0240  * Compare two values, implicitly casting arguments.
0241  */
0242 template<class RealType>
0243 CELER_FUNCTION bool SoftEqual<RealType>::operator()(value_type a,
0244                                                     value_type b) const
0245 {
0246     real_type rel = rel_ * std::fmax(std::fabs(a), std::fabs(b));
0247     return std::fabs(a - b) < std::fmax(abs_, rel);
0248 }
0249 
0250 //---------------------------------------------------------------------------//
0251 /*!
0252  * Construct with default absolute precision.
0253  */
0254 template<class RealType>
0255 CELER_CONSTEXPR_FUNCTION SoftClose<RealType>::SoftClose()
0256     : abs_(SETraits::abs_thresh())
0257 {
0258 }
0259 
0260 //---------------------------------------------------------------------------//
0261 /*!
0262  * Construct with absolute precision.
0263  *
0264  * \param abs threshold for absolute error (default 1.0e-14 for doubles)
0265  */
0266 template<class RealType>
0267 CELER_FUNCTION SoftClose<RealType>::SoftClose(value_type abs) : abs_(abs)
0268 {
0269     CELER_EXPECT(abs > 0);
0270 }
0271 
0272 //---------------------------------------------------------------------------//
0273 /*!
0274  * Compare two values (implicitly casting arguments)
0275  */
0276 template<class RealType>
0277 CELER_FUNCTION bool SoftClose<RealType>::operator()(value_type a,
0278                                                     value_type b) const
0279 {
0280     return std::fabs(a - b) < abs_;
0281 }
0282 
0283 //---------------------------------------------------------------------------//
0284 /*!
0285  * Construct with default absolute precision.
0286  */
0287 template<class RealType>
0288 CELER_CONSTEXPR_FUNCTION SoftZero<RealType>::SoftZero()
0289     : abs_(SETraits::abs_thresh())
0290 {
0291 }
0292 
0293 //---------------------------------------------------------------------------//
0294 /*!
0295  * Construct with absolute precision.
0296  *
0297  * \param abs threshold for absolute error (default 1.0e-14 for doubles)
0298  */
0299 template<class RealType>
0300 CELER_FUNCTION SoftZero<RealType>::SoftZero(value_type abs) : abs_(abs)
0301 {
0302     CELER_EXPECT(abs > 0);
0303 }
0304 
0305 //---------------------------------------------------------------------------//
0306 /*!
0307  * See if the value is within absolute tolerance of zero.
0308  *
0309  * \param actual scalar floating point value
0310  */
0311 template<class RealType>
0312 CELER_FUNCTION bool SoftZero<RealType>::operator()(value_type actual) const
0313 {
0314     return std::fabs(actual) < abs_;
0315 }
0316 
0317 //---------------------------------------------------------------------------//
0318 //! Soft equivalence with default tolerance
0319 template<class RealType>
0320 inline CELER_FUNCTION bool soft_equal(RealType expected, RealType actual)
0321 {
0322     return SoftEqual<RealType>()(expected, actual);
0323 }
0324 
0325 //---------------------------------------------------------------------------//
0326 //! Soft equivalence with relative error
0327 template<class RealType>
0328 inline CELER_FUNCTION bool soft_near(
0329     RealType expected, RealType actual, RealType rel_error)
0330 {
0331     return SoftEqual<RealType>(rel_error)(expected, actual);
0332 }
0333 
0334 //---------------------------------------------------------------------------//
0335 //! Soft equivalence to zero, with default tolerance
0336 template<class RealType>
0337 inline CELER_FUNCTION bool soft_zero(RealType actual)
0338 {
0339     return SoftZero<RealType>()(actual);
0340 }
0341 
0342 //---------------------------------------------------------------------------//
0343 //! Soft modulo operator
0344 template<class RealType>
0345 inline CELER_FUNCTION bool soft_mod(RealType dividend, RealType divisor)
0346 {
0347     auto remainder = std::fmod(dividend, divisor);
0348 
0349     SoftEqual<RealType> seq(detail::SoftEqualTraits<RealType>::rel_prec());
0350 
0351     return seq(0, remainder) || seq(divisor, remainder);
0352 }
0353 
0354 //---------------------------------------------------------------------------//
0355 }  // namespace celeritas