VecGeom/base/Math.h

0001 #ifndef VECGEOM_MATH_H
0002 #define VECGEOM_MATH_H
0003
0004 #include <cmath>
0005 #include <limits>
0006 #include "VecGeom/base/Config.h"
0007 #include "VecCore/Limits.h"
0008
0009 namespace vecgeom {
0010
0011 #ifdef VECGEOM_SINGLE_PRECISION
0012 #define VECCORE_SINGLE_PRECISION
0013 using Precision = float;
0014 #else
0015 using Precision = double;
0016 #endif
0017
0018 inline namespace VECGEOM_IMPL_NAMESPACE {
0019
0020 #ifdef __CUDA_ARCH__
0021 #define VECGEOM_CONST static __constant__ const
0022 #else
0023 #define VECGEOM_CONST static constexpr
0024 #endif
0025
0026 #ifdef VECGEOM_SINGLE_PRECISION
0027 VECGEOM_CONST Precision kTolerance     = 1e-3;
0028 VECGEOM_CONST Precision kInvTolerance  = 1e3;
0029 VECGEOM_CONST Precision kPushTolerance = 1e-3;
0030 VECGEOM_CONST Precision kSqrtTolerance = 3.1622777e-2;
0031 VECGEOM_CONST Precision kAngTolerance  = 1e-2;
0032 VECGEOM_CONST Precision kConeTolerance = 1e-3;
0033 VECGEOM_CONST Precision kFarAway       = 1e5;
0034 #else
0035 VECGEOM_CONST Precision kTolerance     = 1e-9;
0036 VECGEOM_CONST Precision kInvTolerance  = 1e9;
0037 VECGEOM_CONST Precision kPushTolerance = 1e-6;
0038 VECGEOM_CONST Precision kSqrtTolerance = 3.1622777e-5;
0039 VECGEOM_CONST Precision kAngTolerance  = 1e-9;
0040 VECGEOM_CONST Precision kConeTolerance = 1e-7;
0041 VECGEOM_CONST Precision kFarAway       = 1e10;
0042 #endif
0043
0044 // Tolerace distance constant specializations
0045 template <typename Real_t>
0046 constexpr Real_t kToleranceDist = Real_t(0);
0047 template <>
0048 inline constexpr double kToleranceDist<double> = double(1e-9);
0049 template <>
0050 inline constexpr float kToleranceDist<float> = float(1e-3);
0051
0052 template <typename Real_t>
0053 constexpr Real_t kToleranceStrict = Real_t(0);
0054 template <>
0055 inline constexpr double kToleranceStrict<double> = double(1e-9);
0056 template <>
0057 inline constexpr float kToleranceStrict<float> = float(1e-6);
0058
0059 template <typename Real_t>
0060 constexpr Real_t kToleranceDistSquared = Real_t(0);
0061 template <>
0062 inline constexpr double kToleranceDistSquared<double> = kToleranceDist<double> * kToleranceDist<double>;
0063 template <>
0064 inline constexpr float kToleranceDistSquared<float> = kToleranceDist<float> * kToleranceDist<float>;
0065
0066 template <typename Real_t>
0067 constexpr Real_t kToleranceCone = Real_t(0);
0068 template <>
0069 inline constexpr double kToleranceCone<double> = double(1e-7);
0070 template <>
0071 inline constexpr float kToleranceCone<float> = float(1e-3);
0072
0073 template <typename Real_t>
0074 constexpr Real_t kToleranceArb4 = Real_t(0);
0075 template <>
0076 inline constexpr double kToleranceArb4<double> = double(1e-6);
0077 template <>
0078 inline constexpr float kToleranceArb4<float> = float(1e-3);
0079
0080 #ifdef VECGEOM_BVH_SINGLE
0081 constexpr double kToleranceBVH = double(1e-4);
0082 #else
0083 constexpr float kToleranceBVH = float(1e-9);
0084 #endif
0085
0086 using namespace vecCore::math;
0087
0088 VECGEOM_CONST Precision kAvogadro = 6.02214085774e23;
0089 VECGEOM_CONST Precision kEpsilon  = std::numeric_limits<Precision>::epsilon();
0090 template <typename Real_t>
0091 constexpr Real_t kEpsilonT        = std::numeric_limits<Real_t>::epsilon();
0092 VECGEOM_CONST double kInfinityDbl = std::numeric_limits<double>::infinity();
0093 // a function to estimate ULP *unit in the last place for a number
0094 // Compute ULP of a given number x (templated on precision type)
0095 template <typename T>
0096 VECCORE_ATT_HOST_DEVICE T ULP(T x)
0097 {
0098   static_assert(std::is_floating_point<T>::value, "T must be a floating point type");
0099   // Use nextafter to find the next representable value greater than x
0100   T next = std::nextafter(x, vecCore::NumericLimits<T>::Infinity());
0101   return vecCore::math::Abs(next - x);
0102 }
0103
0104 // a special constant to indicate a "miss" length
0105 VECGEOM_CONST Precision kInfLength        = vecCore::NumericLimits<Precision>::Max();
0106 VECGEOM_CONST Precision kMaximum          = vecCore::NumericLimits<Precision>::Max();
0107 VECGEOM_CONST int kMaximumInt             = vecCore::NumericLimits<int>::Max();
0108 VECGEOM_CONST Precision kMinimum          = vecCore::NumericLimits<Precision>::Min();
0109 VECGEOM_CONST Precision kPi               = 3.14159265358979323846;
0110 VECGEOM_CONST Precision kHalfPi           = 0.5 * kPi;
0111 VECGEOM_CONST Precision kTwoPi            = 2. * kPi;
0112 VECGEOM_CONST Precision kTwoPiInv         = 1. / kTwoPi;
0113 VECGEOM_CONST Precision kDegToRad         = kPi / 180.;
0114 VECGEOM_CONST Precision kRadToDeg         = 180. / kPi;
0115 VECGEOM_CONST Precision kRadTolerance     = 1e-9;
0116 VECGEOM_CONST Precision kTiny             = 1e-30;
0117 VECGEOM_CONST Precision kHalfTolerance    = 0.5 * kTolerance;
0118 VECGEOM_CONST Precision kToleranceSquared = kTolerance * kTolerance;
0119
0120 template <typename T>
0121 struct Tiny {
0122   static constexpr T kValue = 1.e-30;
0123 };
0124
0125 template <template <typename, typename> class ImplementationType, typename T, typename Q>
0126 struct Tiny<ImplementationType<T, Q>> {
0127   static constexpr typename ImplementationType<T, Q>::value_type kValue = 1.e-30;
0128 };
0129 } // namespace VECGEOM_IMPL_NAMESPACE
0130 } // namespace vecgeom
0131
0132 #endif