Merge pull request #5798 from InsightSoftwareConsortium/add-constexpr-for-math-operations

add constexpr for math operations

Author: hjmjohnson

Modules/Core/Common/include/itkMath.h

@@ -104,7 +104,7 @@ static constexpr float float_sqrteps = vnl_math::float_sqrteps;
  *  bit, the 64 bit or the vanilla version */
 #define itkTemplateFloatingToIntegerMacro(name)                             \
   template <typename TReturn, typename TInput>                              \
-  inline TReturn name(TInput x)                                             \
+  inline constexpr auto name(TInput x)                                      \
   {                                                                         \
     if constexpr (sizeof(TReturn) <= 4)                                     \
     {                                                                       \
@@ -173,7 +173,7 @@ itkTemplateFloatingToIntegerMacro(RoundHalfIntegerUp)
  *  \sa RoundHalfIntegerUp<TReturn, TInput>()
  */
 template <typename TReturn, typename TInput>
-inline TReturn
+inline constexpr auto
 Round(TInput x)
 {
   return RoundHalfIntegerUp<TReturn, TInput>(x);
@@ -208,7 +208,7 @@ itkTemplateFloatingToIntegerMacro(Ceil)
 #undef itkTemplateFloatingToIntegerMacro
 
 template <typename TReturn, typename TInput>
-inline TReturn
+inline auto
 CastWithRangeCheck(TInput x)
 {
   itkConceptMacro(OnlyDefinedForIntegerTypes1, (itk::Concept::IsInteger<TReturn>));
@@ -247,7 +247,7 @@ CastWithRangeCheck(TInput x)
  * \sa FloatAddULP
  */
 template <typename T>
-inline typename Detail::FloatIEEE<T>::IntType
+inline constexpr typename Detail::FloatIEEE<T>::IntType
 FloatDifferenceULP(T x1, T x2)
 {
   const Detail::FloatIEEE<T> x1f(x1);
@@ -263,7 +263,7 @@ FloatDifferenceULP(T x1, T x2)
  * \sa FloatDifferenceULP
  */
 template <typename T>
-inline T
+inline constexpr T
 FloatAddULP(T x, typename Detail::FloatIEEE<T>::IntType ulps)
 {
   const Detail::FloatIEEE<T> representInput(x);
@@ -714,7 +714,7 @@ NotAlmostEquals(T1 x1, T2 x2)
 
 // The ExactlyEquals function
 template <typename TInput1, typename TInput2>
-inline bool
+inline constexpr bool
 ExactlyEquals(const TInput1 & x1, const TInput2 & x2)
 {
   ITK_GCC_PRAGMA_PUSH
@@ -725,7 +725,7 @@ ExactlyEquals(const TInput1 & x1, const TInput2 & x2)
 
 // The NotExactlyEquals function
 template <typename TInput1, typename TInput2>
-inline bool
+inline constexpr bool
 NotExactlyEquals(const TInput1 & x1, const TInput2 & x2)
 {
   return !ExactlyEquals(x1, x2);
@@ -737,33 +737,64 @@ NotExactlyEquals(const TInput1 & x1, const TInput2 & x2)
  * \note Negative numbers cannot be prime.
  */
 /** @ITKStartGrouping */
-ITKCommon_EXPORT bool
-IsPrime(unsigned short n);
-ITKCommon_EXPORT bool
-IsPrime(unsigned int n);
-ITKCommon_EXPORT bool
-IsPrime(unsigned long n);
-ITKCommon_EXPORT bool
-IsPrime(unsigned long long n);
+template <typename T, typename = std::enable_if_t<std::is_integral_v<T>>>
+constexpr bool
+IsPrime(T n)
+{
+  //   An Introduction to the Theory of Numbers — G.H. Hardy & E.M. Wright, 6th ed., Oxford University Press.
+  //  See Chapter 1 (Elementary Results on Primes). The observation that primes > 3 lie in residue classes
+  //  +/1 mod 6 follows directly from eliminating multiples of 2 and 3.
+  if (n <= 1)
+  {
+    return false;
+  }
+  if (n == 2 || n == 3)
+  {
+    return true;
+  }
+  if (n % 2 == 0 || n % 3 == 0)
+  {
+    return false;
+  }
+  // 6k +/- 1, and avoid overflow via i <= n / i
+  for (T x = 5; x <= n / x; x += 6)
+  {
+    if (n % x == 0 || n % (x + 2) == 0)
+    {
+      return false;
+    }
+  }
+  return true;
+}
 /** @ITKEndGrouping */
 
 /** Return the greatest factor of the decomposition in prime numbers. */
 /** @ITKStartGrouping */
-ITKCommon_EXPORT unsigned short
-GreatestPrimeFactor(unsigned short n);
-ITKCommon_EXPORT unsigned int
-GreatestPrimeFactor(unsigned int n);
-ITKCommon_EXPORT unsigned long
-GreatestPrimeFactor(unsigned long n);
-ITKCommon_EXPORT unsigned long long
-GreatestPrimeFactor(unsigned long long n);
+template <typename T, typename = std::enable_if_t<std::is_integral_v<T>>>
+constexpr T
+GreatestPrimeFactor(T n)
+{
+  T v = 2;
+  while (v <= n)
+  {
+    if (n % v == 0 && IsPrime(v))
+    {
+      n /= v;
+    }
+    else
+    {
+      v += 1;
+    }
+  }
+  return v;
+}
 /** @ITKEndGrouping */
 
 /**  Returns `a * b`. Numeric overflow triggers a compilation error in
  * "constexpr context" and a debug assert failure at run-time.
  */
 template <typename TReturnType = uintmax_t>
-constexpr TReturnType
+constexpr auto
 UnsignedProduct(const uintmax_t a, const uintmax_t b) noexcept
 {
   static_assert(std::is_unsigned_v<TReturnType>, "UnsignedProduct only supports unsigned return types");
@@ -785,8 +816,8 @@ UnsignedProduct(const uintmax_t a, const uintmax_t b) noexcept
  * no agreed-upon value: https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero
  */
 template <typename TReturnType = uintmax_t>
-constexpr TReturnType
-UnsignedPower(const uintmax_t base, const uintmax_t exponent) noexcept
+constexpr auto
+UnsignedPower(const uintmax_t base, const uintmax_t exponent) noexcept -> TReturnType
 {
   static_assert(std::is_unsigned_v<TReturnType>, "UnsignedPower only supports unsigned return types");
 

Modules/Core/Common/include/itkMathDetail.h

@@ -89,7 +89,7 @@ ITK_GCC_PRAGMA_PUSH
 ITK_GCC_SUPPRESS_Wfloat_equal
 
 template <typename TReturn, typename TInput>
-inline TReturn
+constexpr TReturn
 RoundHalfIntegerToEven_base(TInput x)
 {
   if (NumericTraits<TInput>::IsNonnegative(x))
@@ -106,7 +106,7 @@ RoundHalfIntegerToEven_base(TInput x)
 }
 
 template <typename TReturn, typename TInput>
-inline TReturn
+constexpr TReturn
 RoundHalfIntegerUp_base(TInput x)
 {
   x += static_cast<TInput>(0.5);
@@ -116,7 +116,7 @@ RoundHalfIntegerUp_base(TInput x)
 }
 
 template <typename TReturn, typename TInput>
-inline TReturn
+constexpr TReturn
 Floor_base(TInput x)
 {
   const auto r = static_cast<TReturn>(x);
@@ -126,7 +126,7 @@ Floor_base(TInput x)
 }
 
 template <typename TReturn, typename TInput>
-inline TReturn
+constexpr TReturn
 Ceil_base(TInput x)
 {
   const auto r = static_cast<TReturn>(x);
@@ -217,15 +217,11 @@ RoundHalfIntegerToEven_32(float x)
 
 #else // Base implementation
 
-inline int32_t
-RoundHalfIntegerToEven_32(double x)
+template <typename TInput>
+constexpr int32_t
+RoundHalfIntegerToEven_32(TInput x)
 {
-  return RoundHalfIntegerToEven_base<int32_t, double>(x);
-}
-inline int32_t
-RoundHalfIntegerToEven_32(float x)
-{
-  return RoundHalfIntegerToEven_base<int32_t, float>(x);
+  return RoundHalfIntegerToEven_base<int32_t, TInput>(x);
 }
 
 #endif
@@ -267,37 +263,25 @@ Ceil_32(float x)
 
 #else // Base implementation
 
-inline int32_t
-RoundHalfIntegerUp_32(double x)
-{
-  return RoundHalfIntegerUp_base<int32_t, double>(x);
-}
-inline int32_t
-RoundHalfIntegerUp_32(float x)
+template <typename TInput>
+constexpr int32_t
+RoundHalfIntegerUp_32(TInput x)
 {
-  return RoundHalfIntegerUp_base<int32_t, float>(x);
+  return RoundHalfIntegerUp_base<int32_t, TInput>(x);
 }
 
-inline int32_t
-Floor_32(double x)
-{
-  return Floor_base<int32_t, double>(x);
-}
-inline int32_t
-Floor_32(float x)
+template <typename TInput>
+constexpr int32_t
+Floor_32(TInput x)
 {
-  return Floor_base<int32_t, float>(x);
+  return Floor_base<int32_t, TInput>(x);
 }
 
-inline int32_t
-Ceil_32(double x)
-{
-  return Ceil_base<int32_t, double>(x);
-}
-inline int32_t
-Ceil_32(float x)
+template <typename TInput>
+constexpr int32_t
+Ceil_32(TInput x)
 {
-  return Ceil_base<int32_t, float>(x);
+  return Ceil_base<int32_t, TInput>(x);
 }
 
 #endif // USE_SSE2_32IMPL || GCC_USE_ASM_32IMPL || VC_USE_ASM_32IMPL
@@ -383,15 +367,11 @@ RoundHalfIntegerToEven_64(float x)
 
 #else // Base implementation
 
-inline int64_t
-RoundHalfIntegerToEven_64(double x)
-{
-  return RoundHalfIntegerToEven_base<int64_t, double>(x);
-}
-inline int64_t
-RoundHalfIntegerToEven_64(float x)
+template <typename TInput>
+constexpr int64_t
+RoundHalfIntegerToEven_64(TInput x)
 {
-  return RoundHalfIntegerToEven_base<int64_t, float>(x);
+  return RoundHalfIntegerToEven_base<int64_t, TInput>(x);
 }
 
 #endif
@@ -433,37 +413,25 @@ Ceil_64(float x)
 
 #else // Base implementation
 
-inline int64_t
-RoundHalfIntegerUp_64(double x)
-{
-  return RoundHalfIntegerUp_base<int64_t, double>(x);
-}
-inline int64_t
-RoundHalfIntegerUp_64(float x)
+template <typename TInput>
+constexpr int64_t
+RoundHalfIntegerUp_64(TInput x)
 {
-  return RoundHalfIntegerUp_base<int64_t, float>(x);
+  return RoundHalfIntegerUp_base<int64_t, TInput>(x);
 }
 
-inline int64_t
-Floor_64(double x)
-{
-  return Floor_base<int64_t, double>(x);
-}
-inline int64_t
-Floor_64(float x)
+template <typename TInput>
+constexpr int64_t
+Floor_64(TInput x)
 {
-  return Floor_base<int64_t, float>(x);
+  return Floor_base<int64_t, TInput>(x);
 }
 
-inline int64_t
-Ceil_64(double x)
-{
-  return Ceil_base<int64_t, double>(x);
-}
-inline int64_t
-Ceil_64(float x)
+template <typename TInput>
+constexpr int64_t
+Ceil_64(TInput x)
 {
-  return Ceil_base<int64_t, float>(x);
+  return Ceil_base<int64_t, TInput>(x);
 }
 
 #endif // USE_SSE2_64IMPL || GCC_USE_ASM_64IMPL || VC_USE_ASM_64IMPL
@@ -496,18 +464,21 @@ union FloatIEEE
   IntType   asInt;
   UIntType  asUInt;
 
-  FloatIEEE(FloatType f)
+  constexpr FloatIEEE(FloatType f)
     : asFloat(f)
   {}
-  FloatIEEE(IntType i)
+  constexpr FloatIEEE(IntType i)
     : asInt(i)
   {}
-  [[nodiscard]] bool
+  constexpr FloatIEEE(UIntType ui)
+    : asUInt(ui)
+  {}
+  [[nodiscard]] constexpr bool
   Sign() const
   {
     return (asUInt >> (sizeof(asUInt) * 8 - 1)) != 0;
   }
-  [[nodiscard]] IntType
+  [[nodiscard]] constexpr IntType
   AsULP() const
   {
     return this->Sign() ? IntType(~(~UIntType(0) >> 1) - asUInt) : asInt;

Modules/Core/Common/src/CMakeLists.txt

@@ -147,7 +147,6 @@ set(
   itkArrayOutputSpecialization.cxx
   itkNumberToString.cxx
   itkRandomVariateGeneratorBase.cxx
-  itkMath.cxx
   itkProgressTransformer.cxx
   itkSymmetricEigenAnalysis.cxx
   itkExtractImageFilter.cxx