math.hh

Go to the documentation of this file.

#pragma once
 
// DiFfRG
#include <DiFfRG/common/complex_math.hh>
#include <DiFfRG/common/cuda_prefix.hh>
#include <DiFfRG/common/utils.hh>
 
// standard library
#include <cmath>
#include <type_traits>
 
// external libraries
#include <Eigen/Dense>
#include <autodiff/forward/real.hpp>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/vector.h>
 
namespace DiFfRG
{
  template <size_t N, typename T> bool isfinite(const autodiff::Real<N, T> &x)
  {
    return std::isfinite(autodiff::val(x)) && std::isfinite(autodiff::derivative(x));
  }
  using std::isfinite;
 
  template <int n, typename NumberType>
    requires requires(NumberType x) {
      x *x;
      NumberType(1.) / x;
    }
  constexpr __forceinline__ __host__ __device__ NumberType powr(const NumberType x)
  {
    if constexpr (n == 0)
      return NumberType(1.);
    else if constexpr (n < 0)
      return NumberType(1.) / powr<-n, NumberType>(x);
    else if constexpr (n > 1)
      return x * powr<n - 1, NumberType>(x);
    else
      return x;
  }
 
  template <typename NT> constexpr __forceinline__ __host__ __device__ double V_d(NT d)
  {
    using std::pow;
    using std::tgamma;
    return pow(M_PI, d / 2.) / tgamma(d / 2. + 1.);
  }
 
  template <typename NT1, typename NT2> constexpr __forceinline__ __host__ __device__ double V_d(NT1 d, NT2 extent)
  {
    using std::pow;
    using std::tgamma;
    return pow(M_PI, d / 2.) / tgamma(d / 2. + 1.) * pow(extent, d);
  }
 
  template <typename NT> constexpr __forceinline__ __host__ __device__ double S_d(NT d)
  {
    using std::pow;
    using std::tgamma;
    return 2. * pow(M_PI, d / 2.) / tgamma(d / 2.);
  }
 
  template <typename NT> consteval NT S_d_prec(uint d)
  {
    if (d == 1)
      return 2;
    else if (d == 2)
      return static_cast<NT>(2) * static_cast<NT>(M_PI);
    else if (d == 3)
      return static_cast<NT>(4) * static_cast<NT>(M_PI);
    else if (d == 4)
      return static_cast<NT>(2) * powr<2>(static_cast<NT>(M_PI));
    else if (d == 5)
      return static_cast<NT>(8) * powr<2>(static_cast<NT>(M_PI)) / static_cast<NT>(3);
    else if (d == 6)
      return powr<3>(static_cast<NT>(M_PI));
    else if (d == 7)
      return static_cast<NT>(16) * powr<3>(static_cast<NT>(M_PI)) / static_cast<NT>(15);
    return std::numeric_limits<NT>::quiet_NaN();
  }
 
  template <typename NumberType>
    requires requires(NumberType x) { x >= 0; }
  constexpr __forceinline__ __host__ __device__ auto heaviside_theta(const NumberType x)
  {
    if constexpr (std::is_same_v<NumberType, autodiff::real>)
      return x >= 0. ? 1. : 0.;
    else
      return x >= static_cast<NumberType>(0) ? static_cast<NumberType>(1) : static_cast<NumberType>(0);
  }
 
  template <typename NumberType>
    requires requires(NumberType x) { x >= 0; }
  constexpr __forceinline__ __host__ __device__ auto sign(const NumberType x)
  {
    if constexpr (std::is_same_v<NumberType, autodiff::real>)
      return x >= 0. ? 1. : -1.;
    else
      return x >= static_cast<NumberType>(0) ? static_cast<NumberType>(1) : static_cast<NumberType>(-1);
  }
 
  template <typename T1, typename T2, typename T3>
    requires(std::is_floating_point<T1>::value || std::is_same_v<T1, autodiff::real> || is_complex<T1>::value) &&
            (std::is_floating_point<T2>::value || std::is_same_v<T2, autodiff::real> || is_complex<T2>::value) &&
            std::is_floating_point<T3>::value
  bool __forceinline__ __host__ __device__ is_close(T1 a, T2 b, T3 eps_)
  {
    if constexpr (is_complex<T1>::value || is_complex<T2>::value) {
      return is_close(real(a), real(b), eps_) && is_close(imag(a), imag(b), eps_);
    } else if constexpr (std::is_same_v<T1, autodiff::real> || std::is_same_v<T2, autodiff::real>)
      return is_close((double)a, (double)b, (double)eps_);
    else {
      T1 diff = std::fabs(a - b);
      if (diff <= eps_) return true;
      if (diff <= std::fmax(std::fabs(a), std::fabs(b)) * eps_) return true;
    }
    return false;
  }
 
  template <typename T1, typename T2>
    requires(std::is_floating_point<T1>::value || std::is_same_v<T1, autodiff::real> || is_complex<T1>::value) &&
            (std::is_floating_point<T2>::value || std::is_same_v<T2, autodiff::real> || is_complex<T1>::value)
  bool __forceinline__ __host__ __device__ is_close(T1 a, T2 b)
  {
    if constexpr (is_complex<T1>::value || is_complex<T2>::value) {
      return is_close(real(a), real(b)) && is_close(imag(a), imag(b));
    } else if constexpr (std::is_same_v<T1, autodiff::real> || std::is_same_v<T2, autodiff::real>) {
      constexpr auto eps_ = std::numeric_limits<double>::epsilon() * 10.;
      return is_close((double)a, (double)b, eps_);
    } else {
      constexpr auto eps_ = std::max(std::numeric_limits<T1>::epsilon(), std::numeric_limits<T2>::epsilon());
      return is_close(a, b, eps_);
    }
    return false;
  }
 
  template <uint n, typename NT, typename A1, typename A2>
    requires requires(A1 a1, A2 a2) { a1[0] * a2[0]; }
  NT dot(const A1 &a1, const A2 &a2)
  {
    NT ret = a1[0] * a2[0];
    for (uint i = 1; i < n; ++i)
      ret += a1[i] * a2[i];
    return ret;
  }
 
  void dealii_to_eigen(const dealii::Vector<double> &dealii, Eigen::VectorXd &eigen);
 
  void dealii_to_eigen(const dealii::BlockVector<double> &dealii, Eigen::VectorXd &eigen);
 
  void eigen_to_dealii(const Eigen::VectorXd &eigen, dealii::Vector<double> &dealii);
 
  void eigen_to_dealii(const Eigen::VectorXd &eigen, dealii::BlockVector<double> &dealii);
} // namespace DiFfRG