Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-05 08:35:21

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