itkMath.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
/*=========================================================================
 *
 *  Portions of this file are subject to the VTK Toolkit Version 3 copyright.
 *
 *  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 *
 *  For complete copyright, license and disclaimer of warranty information
 *  please refer to the NOTICE file at the top of the ITK source tree.
 *
 *=========================================================================*/
#ifndef itkMath_h
#define itkMath_h
 
#include <cassert>
#include <complex>
#include <cmath>
#include <limits>
#include <type_traits>
// _MSVC_LANG tracks /std: even when __cplusplus stays 199711L (no /Zc:__cplusplus),
// matching the STL condition under which <limits> defines __cpp_lib_math_constants.
#if (defined(__cplusplus) && __cplusplus >= 202002L) || (defined(_MSVC_LANG) && _MSVC_LANG >= 202002L)
#  include <numbers>
#endif
#include "itkMathDetail.h"
#include "itkConceptChecking.h"
#if !defined(ITK_FUTURE_LEGACY_REMOVE)
// itk::Math no longer uses vnl_math; kept transitionally for downstream vnl_math:: users.
#  include <vnl/vnl_math.h>
#endif
 
namespace itk::Math
{
// These mathematical constants originate from VXL's vnl_math.h: correctly-rounded
// double-precision literals (OEIS references in vnl_math.h), defined identically under
// every C++ standard level and verified against C++20 <numbers> below.

inline constexpr double e = 2.71828182845904523536;
inline constexpr double log2e = 1.44269504088896340736;
inline constexpr double log10e = 0.43429448190325182765;
inline constexpr double ln2 = 0.69314718055994530942;
inline constexpr double ln10 = 2.30258509299404568402;
inline constexpr double pi = 3.14159265358979323846;
inline constexpr double twopi = 6.28318530717958647693;
inline constexpr double pi_over_2 = 1.57079632679489661923;
inline constexpr double pi_over_4 = 0.78539816339744830962;
inline constexpr double pi_over_180 = 0.01745329251994329577;
inline constexpr double one_over_pi = 0.31830988618379067154;
inline constexpr double two_over_pi = 0.63661977236758134308;
inline constexpr double deg_per_rad = 57.2957795130823208768;
inline constexpr double sqrt2pi = 2.50662827463100050242;
inline constexpr double two_over_sqrtpi = 1.12837916709551257390;
inline constexpr double one_over_sqrt2pi = 0.39894228040143267794;
inline constexpr double sqrt2 = 1.41421356237309504880;
inline constexpr double sqrt1_2 = 0.70710678118654752440;
inline constexpr double sqrt1_3 = 0.57735026918962576451;
inline constexpr double euler = 0.57721566490153286061;
 
#if defined(__cpp_lib_math_constants)
// Bit-exactness proof against C++20 <numbers>. sqrt2pi and one_over_sqrt2pi have no
// exact std::numbers identity (inv_sqrtpi / sqrt2 double-rounds 1 ULP off).
static_assert(e == std::numbers::e);
static_assert(log2e == std::numbers::log2e);
static_assert(log10e == std::numbers::log10e);
static_assert(ln2 == std::numbers::ln2);
static_assert(ln10 == std::numbers::ln10);
static_assert(pi == std::numbers::pi);
static_assert(twopi == 2 * std::numbers::pi);
static_assert(pi_over_2 == std::numbers::pi / 2);
static_assert(pi_over_4 == std::numbers::pi / 4);
static_assert(pi_over_180 == std::numbers::pi / 180);
static_assert(one_over_pi == std::numbers::inv_pi);
static_assert(two_over_pi == 2 * std::numbers::inv_pi);
static_assert(deg_per_rad == 180 * std::numbers::inv_pi);
static_assert(two_over_sqrtpi == 2 * std::numbers::inv_sqrtpi);
static_assert(sqrt2 == std::numbers::sqrt2);
static_assert(sqrt1_2 == std::numbers::sqrt2 / 2);
static_assert(sqrt1_3 == std::numbers::inv_sqrt3);
static_assert(euler == std::numbers::egamma);
#endif
 
//: IEEE double machine precision
inline constexpr double eps = std::numeric_limits<double>::epsilon();
// sqrt(2^-52) = 2^-26 = 1.490116119384765625e-8 exactly, matching C++26 constexpr std::sqrt(eps);
// vnl_math::sqrteps used 1.490116119384766e-08, which rounds 1 ULP above (2^-26 + 2^-78).
inline constexpr double sqrteps = 0x1p-26;
//: IEEE single machine precision
inline constexpr float float_eps = std::numeric_limits<float>::epsilon();
inline constexpr float float_sqrteps = 3.4526698300e-4F;

#define itkTemplateFloatingToIntegerMacro(name)                             \
  template <typename TReturn, typename TInput>                              \
  inline constexpr auto name(TInput x)                                      \
  {                                                                         \
    if constexpr (sizeof(TReturn) <= 4)                                     \
    {                                                                       \
      return static_cast<TReturn>(Detail::name##_32(x));                    \
    }                                                                       \
    else if constexpr (sizeof(TReturn) <= 8)                                \
    {                                                                       \
      return static_cast<TReturn>(Detail::name##_64(x));                    \
    }                                                                       \
    else                                                                    \
    {                                                                       \
      return static_cast<TReturn>(Detail::name##_base<TReturn, TInput>(x)); \
    }                                                                       \
  }

itkTemplateFloatingToIntegerMacro(RoundHalfIntegerToEven)
 

itkTemplateFloatingToIntegerMacro(RoundHalfIntegerUp)

template <typename TReturn, typename TInput>
inline constexpr auto
Round(TInput x)
{
  return RoundHalfIntegerUp<TReturn, TInput>(x);
}

itkTemplateFloatingToIntegerMacro(Floor)
 

itkTemplateFloatingToIntegerMacro(Ceil)
 
#undef itkTemplateFloatingToIntegerMacro
 
template <typename TReturn, typename TInput>
inline auto
CastWithRangeCheck(TInput x)
{
  itkConceptMacro(OnlyDefinedForIntegerTypes1, (itk::Concept::IsInteger<TReturn>));
  itkConceptMacro(OnlyDefinedForIntegerTypes2, (itk::Concept::IsInteger<TInput>));
 
  auto ret = static_cast<TReturn>(x);
  if constexpr (sizeof(TReturn) > sizeof(TInput) &&
                !(!itk::NumericTraits<TReturn>::is_signed && itk::NumericTraits<TInput>::is_signed))
  {
    // if the output type is bigger and we are not converting a signed
    // integer to an unsigned integer then we have no problems
    return ret;
  }
  else if constexpr (sizeof(TReturn) >= sizeof(TInput))
  {
    if (itk::NumericTraits<TInput>::IsPositive(x) != itk::NumericTraits<TReturn>::IsPositive(ret))
    {
      itk::RangeError _e(__FILE__, __LINE__);
      throw _e;
    }
  }
  else if (static_cast<TInput>(ret) != x ||
           (itk::NumericTraits<TInput>::IsPositive(x) != itk::NumericTraits<TReturn>::IsPositive(ret)))
  {
    itk::RangeError _e(__FILE__, __LINE__);
    throw _e;
  }
  return ret;
}

template <typename T>
inline constexpr typename Detail::FloatIEEE<T>::IntType
FloatDifferenceULP(T x1, T x2)
{
  const Detail::FloatIEEE<T> x1f(x1);
  const Detail::FloatIEEE<T> x2f(x2);
  return x1f.AsULP() - x2f.AsULP();
}

template <typename T>
inline constexpr T
FloatAddULP(T x, typename Detail::FloatIEEE<T>::IntType ulps)
{
  const Detail::FloatIEEE<T> representInput(x);
  const Detail::FloatIEEE<T> representOutput(representInput.asInt + ulps);
  return representOutput.asFloat;
}

template <typename T>
inline bool
FloatAlmostEqual(T                                        x1,
                 T                                        x2,
                 typename Detail::FloatIEEE<T>::IntType   maxUlps = 4,
                 typename Detail::FloatIEEE<T>::FloatType maxAbsoluteDifference = 0.1 *
                                                                                  itk::NumericTraits<T>::epsilon())
{
  // Check if the numbers are really close -- needed
  // when comparing numbers near zero.
  const T absDifference = std::abs(x1 - x2);
  if (absDifference <= maxAbsoluteDifference)
  {
    return true;
  }
 
  // This check for different signs is necessary for several reasons, see the blog link above.
  // Subtracting the signed-magnitude representation of floats using twos-complement
  // math isn't particularly meaningful, and the subtraction would produce a 33-bit
  // result and overflow an int.
  if (std::signbit(x1) != std::signbit(x2))
  {
    return false;
  }
 
  const typename Detail::FloatIEEE<T>::IntType ulps = std::abs(FloatDifferenceULP(x1, x2));
  return ulps <= maxUlps;
}
 
// The following code cannot be moved to the itkMathDetail.h file without introducing circular dependencies
namespace Detail // The Detail namespace holds the templates used by AlmostEquals
{
// The following structs and templates are used to choose
// which version of the AlmostEquals function
// should be implemented base on input parameter types
 
// Structs for choosing AlmostEquals function
 
struct AlmostEqualsFloatVsFloat
{
  template <typename TFloatType1, typename TFloatType2>
  static bool
  AlmostEqualsFunction(TFloatType1 x1, TFloatType2 x2)
  {
    return FloatAlmostEqual<double>(x1, x2);
  }
 
  template <typename TFloatType1, typename TFloatType2>
  static bool
  AlmostEqualsFunction(double x1, double x2)
  {
    return FloatAlmostEqual<double>(x1, x2);
  }
 
  template <typename TFloatType1, typename TFloatType2>
  static bool
  AlmostEqualsFunction(double x1, float x2)
  {
    return FloatAlmostEqual<float>(x1, x2);
  }
 
  template <typename TFloatType1, typename TFloatType2>
  static bool
  AlmostEqualsFunction(float x1, double x2)
  {
    return FloatAlmostEqual<float>(x1, x2);
  }
 
  template <typename TFloatType1, typename TFloatType2>
  static bool
  AlmostEqualsFunction(float x1, float x2)
  {
    return FloatAlmostEqual<float>(x1, x2);
  }
};
 
struct AlmostEqualsFloatVsInteger
{
  template <typename TFloatType, typename TIntType>
  static bool
  AlmostEqualsFunction(TFloatType floatingVariable, TIntType integerVariable)
  {
    return FloatAlmostEqual<TFloatType>(floatingVariable, integerVariable);
  }
};
 
struct AlmostEqualsIntegerVsFloat
{
  template <typename TIntType, typename TFloatType>
  static bool
  AlmostEqualsFunction(TIntType integerVariable, TFloatType floatingVariable)
  {
    return AlmostEqualsFloatVsInteger::AlmostEqualsFunction(floatingVariable, integerVariable);
  }
};
 
struct AlmostEqualsSignedVsUnsigned
{
  template <typename TSignedInt, typename TUnsignedInt>
  static bool
  AlmostEqualsFunction(TSignedInt signedVariable, TUnsignedInt unsignedVariable)
  {
    if (signedVariable < 0)
    {
      return false;
    }
    if (unsignedVariable > static_cast<size_t>(itk::NumericTraits<TSignedInt>::max()))
    {
      return false;
    }
    return signedVariable == static_cast<TSignedInt>(unsignedVariable);
  }
};
 
struct AlmostEqualsUnsignedVsSigned
{
  template <typename TUnsignedInt, typename TSignedInt>
  static bool
  AlmostEqualsFunction(TUnsignedInt unsignedVariable, TSignedInt signedVariable)
  {
    return AlmostEqualsSignedVsUnsigned::AlmostEqualsFunction(signedVariable, unsignedVariable);
  }
};
 
struct AlmostEqualsPlainOldEquals
{
  template <typename TIntegerType1, typename TIntegerType2>
  static bool
  AlmostEqualsFunction(TIntegerType1 x1, TIntegerType2 x2)
  {
    return x1 == x2;
  }
};
// end of structs that choose the specific AlmostEquals function
 
// Selector structs, these select the correct case based on its types
//        input1 is int?  input 1 is signed? input2 is int?  input 2 is signed?
template <bool TInput1IsIntger, bool TInput1IsSigned, bool TInput2IsInteger, bool TInput2IsSigned>
struct AlmostEqualsFunctionSelector
{ // default case
  using SelectedVersion = AlmostEqualsPlainOldEquals;
};

template <>
struct AlmostEqualsFunctionSelector<false, true, false, true>
// floating type v floating type
{
  using SelectedVersion = AlmostEqualsFloatVsFloat;
};
 
template <>
struct AlmostEqualsFunctionSelector<false, true, true, true>
// float vs int
{
  using SelectedVersion = AlmostEqualsFloatVsInteger;
};
 
template <>
struct AlmostEqualsFunctionSelector<false, true, true, false>
// float vs unsigned int
{
  using SelectedVersion = AlmostEqualsFloatVsInteger;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, false, false, true>
// unsigned int vs float
{
  using SelectedVersion = AlmostEqualsIntegerVsFloat;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, true, false, true>
// int vs float
{
  using SelectedVersion = AlmostEqualsIntegerVsFloat;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, true, true, false>
// signed vs unsigned
{
  using SelectedVersion = AlmostEqualsSignedVsUnsigned;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, false, true, true>
// unsigned vs signed
{
  using SelectedVersion = AlmostEqualsUnsignedVsSigned;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, true, true, true>
//   signed vs signed
{
  using SelectedVersion = AlmostEqualsPlainOldEquals;
};
 
template <>
struct AlmostEqualsFunctionSelector<true, false, true, false>
// unsigned vs unsigned
{
  using SelectedVersion = AlmostEqualsPlainOldEquals;
};
// end of AlmostEqualsFunctionSelector structs
 
// The implementor tells the selector what to do
template <typename TInputType1, typename TInputType2>
struct AlmostEqualsScalarImplementer
{
  static constexpr bool TInputType1IsInteger = std::is_integral_v<TInputType1>;
  static constexpr bool TInputType1IsSigned = std::is_signed_v<TInputType1>;
  static constexpr bool TInputType2IsInteger = std::is_integral_v<TInputType2>;
  static constexpr bool TInputType2IsSigned = std::is_signed_v<TInputType2>;
 
  using SelectedVersion = typename AlmostEqualsFunctionSelector<TInputType1IsInteger,
                                                                TInputType1IsSigned,
                                                                TInputType2IsInteger,
                                                                TInputType2IsSigned>::SelectedVersion;
};
 
// The AlmostEqualsScalarComparer returns the result of an
// approximate comparison between two scalar values of
// potentially different data types.
template <typename TScalarType1, typename TScalarType2>
inline bool
AlmostEqualsScalarComparer(TScalarType1 x1, TScalarType2 x2)
{
  return AlmostEqualsScalarImplementer<TScalarType1, TScalarType2>::SelectedVersion::
    template AlmostEqualsFunction<TScalarType1, TScalarType2>(x1, x2);
}
 
// The following structs are used to evaluate approximate comparisons between
// complex and scalar values of potentially different types.
 
// Comparisons between scalar types use the AlmostEqualsScalarComparer function.
struct AlmostEqualsScalarVsScalar
{
  template <typename TScalarType1, typename TScalarType2>
  static bool
  AlmostEqualsFunction(TScalarType1 x1, TScalarType2 x2)
  {
    return AlmostEqualsScalarComparer(x1, x2);
  }
};
 
// Comparisons between two complex values compare the real and imaginary components
// separately with the AlmostEqualsScalarComparer function.
struct AlmostEqualsComplexVsComplex
{
  template <typename TComplexType1, typename TComplexType2>
  static bool
  AlmostEqualsFunction(TComplexType1 x1, TComplexType2 x2)
  {
    return AlmostEqualsScalarComparer(x1.real(), x2.real()) && AlmostEqualsScalarComparer(x1.imag(), x2.imag());
  }
};
 
// Comparisons between complex and scalar values first check to see if the imaginary component
// of the complex value is approximately 0. Then a ScalarComparison is done between the real
// part of the complex number and the scalar value.
struct AlmostEqualsScalarVsComplex
{
  template <typename TScalarType, typename TComplexType>
  static bool
  AlmostEqualsFunction(TScalarType scalarVariable, TComplexType complexVariable)
  {
    if (!AlmostEqualsScalarComparer(complexVariable.imag(), typename itk::NumericTraits<TComplexType>::ValueType{}))
    {
      return false;
    }
    return AlmostEqualsScalarComparer(scalarVariable, complexVariable.real());
  }
};
 
struct AlmostEqualsComplexVsScalar
{
  template <typename TComplexType, typename TScalarType>
  static bool
  AlmostEqualsFunction(TComplexType complexVariable, TScalarType scalarVariable)
  {
    return AlmostEqualsScalarVsComplex::AlmostEqualsFunction(scalarVariable, complexVariable);
  }
};
 
// The AlmostEqualsComplexChooser structs choose the correct case
// from the input parameter types' IsComplex property
// The default case is scalar vs scalar
template <bool T1IsComplex, bool T2IsComplex> // Default is false, false
struct AlmostEqualsComplexChooser
{
  using ChosenVersion = AlmostEqualsScalarVsScalar;
};
 
template <>
struct AlmostEqualsComplexChooser<true, true>
{
  using ChosenVersion = AlmostEqualsComplexVsComplex;
};
 
template <>
struct AlmostEqualsComplexChooser<false, true>
{
  using ChosenVersion = AlmostEqualsScalarVsComplex;
};
 
template <>
struct AlmostEqualsComplexChooser<true, false>
{
  using ChosenVersion = AlmostEqualsComplexVsScalar;
};
// End of AlmostEqualsComplexChooser structs.
 
// The AlmostEqualsComplexImplementer determines which of the input
// parameters are complex and which are real, and sends that information
// to the AlmostEqualsComplexChooser structs to determine the proper
// type of evaluation.
template <typename T1, typename T2>
struct AlmostEqualsComplexImplementer
{
  static constexpr bool T1IsComplex = NumericTraits<T1>::IsComplex;
  static constexpr bool T2IsComplex = NumericTraits<T2>::IsComplex;
 
  using ChosenVersion = typename AlmostEqualsComplexChooser<T1IsComplex, T2IsComplex>::ChosenVersion;
};
 
} // end namespace Detail

 
// The AlmostEquals function
template <typename T1, typename T2>
inline bool
AlmostEquals(T1 x1, T2 x2)
{
  return Detail::AlmostEqualsComplexImplementer<T1, T2>::ChosenVersion::AlmostEqualsFunction(x1, x2);
}
 
// The NotAlmostEquals function
template <typename T1, typename T2>
inline bool
NotAlmostEquals(T1 x1, T2 x2)
{
  return !AlmostEquals(x1, x2);
}
 

 
// The ExactlyEquals function
template <typename TInput1, typename TInput2>
inline constexpr bool
ExactlyEquals(const TInput1 & x1, const TInput2 & x2)
{
  ITK_GCC_PRAGMA_PUSH
  ITK_GCC_SUPPRESS_Wfloat_equal
  return x1 == x2;
  ITK_GCC_PRAGMA_POP
}
 
// The NotExactlyEquals function
template <typename TInput1, typename TInput2>
inline constexpr bool
NotExactlyEquals(const TInput1 & x1, const TInput2 & x2)
{
  return !ExactlyEquals(x1, x2);
}
 

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;
}


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;
}


template <typename TReturnType = uintmax_t>
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");
 
  // Note that numeric overflow is not "undefined behavior", for unsigned numbers.
  // This function checks if the result of a*b is mathematically correct.
  return (a == 0) || (b == 0) ||
             (((static_cast<TReturnType>(a * b) / a) == b) && ((static_cast<TReturnType>(a * b) / b) == a))
           ? static_cast<TReturnType>(a * b)
           : (assert(!"UnsignedProduct overflow!"), 0);
}
 

template <typename TReturnType = uintmax_t>
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");
 
  // Uses recursive function calls because C++11 does not support other ways of
  // iterations for a constexpr function.
  return (exponent == 0)   ? (assert(base > 0), 1)
         : (exponent == 1) ? base
                           : UnsignedProduct<TReturnType>(UnsignedPower<TReturnType>(base, exponent / 2),
                                                          UnsignedPower<TReturnType>(base, (exponent + 1) / 2));
}
 
 
/*==========================================
 * itk::Math utility functions. std:: equivalents are used where they exist;
 * the remainder are native definitions providing the historical interface.
 */
using std::isnan;
using std::isinf;
using std::isfinite;
using std::isnormal;
using std::cbrt;
using std::hypot;

template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
constexpr int
sgn(T x)
{
  return (x != T{}) ? ((x > T{}) ? 1 : -1) : 0;
}

template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
constexpr int
sgn0(T x)
{
  return (x >= T{}) ? 1 : -1;
}

constexpr bool
sqr(bool x)
{
  return x;
}
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
constexpr T
sqr(T x)
{
  return x * x;
}

constexpr bool
cube(bool x)
{
  return x;
}
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
constexpr T
cube(T x)
{
  return x * x * x;
}

constexpr unsigned int
squared_magnitude(char x)
{
  return static_cast<unsigned int>(static_cast<int>(x) * static_cast<int>(x));
}
constexpr unsigned int
squared_magnitude(unsigned char x)
{
  return static_cast<unsigned int>(static_cast<int>(x) * static_cast<int>(x));
}
constexpr unsigned int
squared_magnitude(int x)
{
  // Multiply in unsigned: avoids the signed-overflow UB of vnl_math's x * x.
  const auto ux = static_cast<unsigned int>(x);
  return ux * ux;
}
constexpr unsigned int
squared_magnitude(unsigned int x)
{
  return x * x;
}
constexpr unsigned long
squared_magnitude(long x)
{
  const auto ux = static_cast<unsigned long>(x);
  return ux * ux;
}
constexpr unsigned long
squared_magnitude(unsigned long x)
{
  return x * x;
}
constexpr unsigned long long
squared_magnitude(long long x)
{
  const auto ux = static_cast<unsigned long long>(x);
  return ux * ux;
}
constexpr unsigned long long
squared_magnitude(unsigned long long x)
{
  return x * x;
}
constexpr float
squared_magnitude(float x)
{
  return x * x;
}
constexpr double
squared_magnitude(double x)
{
  return x * x;
}
constexpr long double
squared_magnitude(long double x)
{
  return x * x;
}

template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
constexpr T
remainder_truncated(T x, T y)
{
  return x % y;
}
inline float
remainder_truncated(float x, float y)
{
  return std::fmod(x, y);
}
inline double
remainder_truncated(double x, double y)
{
  return std::fmod(x, y);
}
inline long double
remainder_truncated(long double x, long double y)
{
  return std::fmod(x, y);
}

template <typename T, std::enable_if_t<std::is_integral_v<T> && std::is_signed_v<T>, int> = 0>
constexpr T
remainder_floored(T x, T y)
{
  // Conditional add: vnl_math's ((x % y) + y) % y overflows for |y| > max/2.
  const T r = x % y;
  return (r != 0 && ((r < 0) != (y < 0))) ? r + y : r;
}
template <typename T, std::enable_if_t<std::is_integral_v<T> && std::is_unsigned_v<T>, int> = 0>
constexpr T
remainder_floored(T x, T y)
{
  return x % y;
}
inline float
remainder_floored(float x, float y)
{
  return std::fmod(std::fmod(x, y) + y, y);
}
inline double
remainder_floored(double x, double y)
{
  return std::fmod(std::fmod(x, y) + y, y);
}
inline long double
remainder_floored(long double x, long double y)
{
  return std::fmod(std::fmod(x, y) + y, y);
}

inline int
rnd_halfinttoeven(float x)
{
  return static_cast<int>(std::lrint(x));
}
inline int
rnd_halfinttoeven(double x)
{
  return static_cast<int>(std::lrint(x));
}

inline int
rnd_halfintup(float x)
{
  return rnd_halfinttoeven(2 * x + 0.5f) >> 1;
}
inline int
rnd_halfintup(double x)
{
  return rnd_halfinttoeven(2 * x + 0.5) >> 1;
}

inline int
rnd(float x)
{
  return static_cast<int>(std::lrint(x));
}
inline int
rnd(double x)
{
  return static_cast<int>(std::lrint(x));
}

inline int
floor(float x)
{
  return static_cast<int>(std::floor(x));
}
inline int
floor(double x)
{
  return static_cast<int>(std::floor(x));
}

inline int
ceil(float x)
{
  return static_cast<int>(std::ceil(x));
}
inline int
ceil(double x)
{
  return static_cast<int>(std::ceil(x));
}

inline double
angle_0_to_2pi(double angle)
{
  angle = std::fmod(angle, twopi);
  if (angle >= 0)
    return angle;
  const double a = angle + twopi;
  // Guard the boundary: tiny negative inputs round a up to exactly twopi; clamp 1 ULP below.
  return (a < twopi) ? a : 6.28318530717958575;
}

inline double
angle_minuspi_to_pi(double angle)
{
  angle = std::fmod(angle, twopi);
  if (angle > pi)
    angle -= twopi;
  if (angle < -pi)
    angle += twopi;
  return angle;
}
 

template <typename T>
constexpr auto
Absolute(T x) noexcept
{
  if constexpr (std::is_same_v<T, bool>)
  {
    // std::abs does not provide abs for bool, but ITK depends on bool support for abs.
    return x;
  }
  else if constexpr (std::is_integral_v<T>)
  {
    if constexpr (std::is_signed_v<T>)
    {
      using UnsignedType = std::make_unsigned_t<T>;
      using Limits = std::numeric_limits<T>;
 
      // The absolute value of the minimum integer value, having a two's complement signed integer representation.
      constexpr UnsignedType absoluteValueOfMin{ UnsignedType{ Limits::max() } + UnsignedType{ 1 } };
 
      return static_cast<UnsignedType>((x < 0) ? (x == Limits::min()) ? absoluteValueOfMin : -x : x);
    }
    else
    { // In C++17, the std::abs() integer overloads are only for : int, long, and long long.
      // unsigned type values std::abs() are supported through implicit type conversions to
      // the next larger sized integer type.  The 'unsigned long long' failed due to an
      // lack of larger sized signed integer type to be promoted to.
      // type of x  | default std::abs() called
      // uint8_t    | std::abs( (int32_t) x)
      // uint16_t   | std::abs( (int32_t) x)
      // uint32_t   | std::abs( (int64_t) x)
      // uint64_t   | // compiler error
      // This explicit override provides direct support for all unsigned integer types with type promotion.
      return x;
    }
  }
  else if constexpr (std::is_floating_point_v<T>) // floating point or complex<> types
  {                                               // floating point std::abs is constexpr in c++23 and later
    // Note: +0.0 and -0.0 are considered equal, according to ExactlyEquals. They are both considered equal to T{}.
    // And they both have the same absolute value: +0.0.
    return ExactlyEquals(x, T{}) ? T{} : (x < 0) ? -x : x;
  }
}
template <typename T>
auto
Absolute(const std::complex<T> & x) noexcept
{
  // std::abs<T>(std::complex<T>) is not constexpr in in any proposed c++ standard
  return std::abs<T>(x);
}
 
#if !defined(ITK_FUTURE_LEGACY_REMOVE)
// itk::Math::Absolute has different behavior than std::abs.  The use of
// abs() in the ITK context without namespace resolution can be confusing
// The abs() version is provided for backwards compatibility.
template <typename T>
auto
abs(T x) noexcept
{
  return Absolute(x);
}
#endif
 
} // namespace itk::Math
 
#endif // end of itkMath.h