Loading...
Searching...
No Matches
math.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_MATH_HH
#define DUNE_MATH_HH
#include <cmath>
#include <complex>
#include <limits>
#include <type_traits>
#include <dune/common/typeutilities.hh>
namespace Dune
{
template< class T >
struct StandardMathematicalConstants
{
static const T e ()
{
using std::exp;
static const T e = exp( T( 1 ) );
return e;
}
static const T pi ()
{
using std::acos;
static const T pi = acos( T( -1 ) );
return pi;
}
};
template< class Field >
struct MathematicalConstants
: public StandardMathematicalConstants<Field>
{};
template <class Base, class Exponent>
constexpr Base power(Base m, Exponent p)
{
static_assert(std::numeric_limits<Exponent>::is_integer, "Exponent must be an integer type!");
auto result = Base(1);
auto absp = (p<0) ? -p : p; // This is simply abs, but std::abs is not constexpr
for (Exponent i = Exponent(0); i<absp; i++)
result *= m;
if (p<0)
result = Base(1)/result;
return result;
}
// T has to be an integral type
template<class T>
constexpr inline static T factorial(const T& n) noexcept
{
static_assert(std::numeric_limits<T>::is_integer, "`factorial(n)` has to be called with an integer type.");
T fac = 1;
for(T k = 0; k < n; ++k)
fac *= k+1;
return fac;
}
template<class T, T n>
constexpr inline static auto factorial (std::integral_constant<T, n>) noexcept
{
return std::integral_constant<T, factorial(n)>{};
}
// T has to be an integral type
template<class T>
constexpr inline static T binomial (const T& n, const T& k) noexcept
{
static_assert(std::numeric_limits<T>::is_integer, "`binomial(n, k)` has to be called with an integer type.");
if( k < 0 || k > n )
return 0;
if (2*k > n)
return binomial(n, n-k);
T bin = 1;
for(auto i = n-k; i < n; ++i)
bin *= i+1;
return bin / factorial(k);
}
template<class T, T n, T k>
constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, k>) noexcept
{
return std::integral_constant<T, binomial(n, k)>{};
}
template<class T, T n>
constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, n>) noexcept
{
return std::integral_constant<T, (n >= 0 ? 1 : 0)>{};
}
// conjugate complex does nothing for non-complex types
template<class K>
constexpr K conjugateComplex (const K& x)
{
return x;
}
#ifndef DOXYGEN
// specialization for complex
template<class K>
constexpr std::complex<K> conjugateComplex (const std::complex<K>& c)
{
return std::complex<K>(c.real(),-c.imag());
}
#endif
template <class T>
constexpr int sign(const T& val)
{
return (val < 0 ? -1 : 1);
}
namespace Impl {
// Returns whether a given type behaves like std::complex<>, i.e. whether
// real() and imag() are defined
template<class T>
struct isComplexLike {
private:
template<class U>
static auto test(U* u) -> decltype(u->real(), u->imag(), std::true_type());
template<class U>
static auto test(...) -> decltype(std::false_type());
public:
static const bool value = decltype(test<T>(0))::value;
};
} // namespace Impl
namespace MathOverloads {
struct ADLTag {};
#define DUNE_COMMON_MATH_ISFUNCTION(function, stdfunction) \
template<class T> \
auto function(const T &t, PriorityTag<1>, ADLTag) \
-> decltype(function(t)) { \
return function(t); \
} \
template<class T> \
auto function(const T &t, PriorityTag<0>, ADLTag) { \
using std::stdfunction; \
return stdfunction(t); \
} \
static_assert(true, "Require semicolon to unconfuse editors")
DUNE_COMMON_MATH_ISFUNCTION(isNaN,isnan);
DUNE_COMMON_MATH_ISFUNCTION(isInf,isinf);
DUNE_COMMON_MATH_ISFUNCTION(isFinite,isfinite);
#undef DUNE_COMMON_MATH_ISFUNCTION
template<class T>
auto isUnordered(const T &t1, const T &t2, PriorityTag<1>, ADLTag)
-> decltype(isUnordered(t1, t2)) {
return isUnordered(t1, t2);
}
template<class T>
auto isUnordered(const T &t1, const T &t2, PriorityTag<0>, ADLTag) {
using std::isunordered;
return isunordered(t1, t2);
}
}
namespace MathImpl {
// NOTE: it is important that these functors have names different from the
// names of the functions they are forwarding to. Otherwise the
// unqualified call would find the functor type, not a function, and ADL
// would never be attempted.
#define DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(function) \
struct function##Impl { \
template<class T> \
constexpr auto operator()(const T &t) const { \
return function(t, PriorityTag<10>{}, MathOverloads::ADLTag{}); \
} \
}; \
static_assert(true, "Require semicolon to unconfuse editors")
DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isNaN);
DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isInf);
DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isFinite);
#undef DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR
struct isUnorderedImpl {
template<class T>
constexpr auto operator()(const T &t1, const T &t2) const {
return isUnordered(t1, t2, PriorityTag<10>{}, MathOverloads::ADLTag{});
}
};
} //MathImpl
namespace Impl {
/* This helper has a math functor as a static constexpr member. Doing
this as a static member of a template struct means we can do this
without violating the ODR or putting the definition into a separate
compilation unit, while still still ensuring the functor is the same
lvalue across all compilation units.
*/
template<class T>
struct MathDummy
{
static constexpr T value{};
};
template<class T>
constexpr T MathDummy<T>::value;
} //namespace Impl
namespace {
/* Provide the math functors directly in the `Dune` namespace.
This actually declares a different name in each translation unit, but
they all resolve to the same lvalue.
*/
constexpr auto const &isNaN = Impl::MathDummy<MathImpl::isNaNImpl>::value;
constexpr auto const &isInf = Impl::MathDummy<MathImpl::isInfImpl>::value;
constexpr auto const &isFinite = Impl::MathDummy<MathImpl::isFiniteImpl>::value;
constexpr auto const &isUnordered = Impl::MathDummy<MathImpl::isUnorderedImpl>::value;
}
namespace MathOverloads {
/*Overloads for complex types*/
template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
auto isNaN(const T &t, PriorityTag<2>, ADLTag) {
return Dune::isNaN(real(t)) || Dune::isNaN(imag(t));
}
template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
auto isInf(const T &t, PriorityTag<2>, ADLTag) {
return Dune::isInf(real(t)) || Dune::isInf(imag(t));
}
template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
auto isFinite(const T &t, PriorityTag<2>, ADLTag) {
return Dune::isFinite(real(t)) && Dune::isFinite(imag(t));
}
} //MathOverloads
}
#endif // #ifndef DUNE_MATH_HH
typeutilities.hh
Utilities for type computations, constraining overloads, ...
DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR
#define DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(function)
Definition math.hh:251
DUNE_COMMON_MATH_ISFUNCTION
#define DUNE_COMMON_MATH_ISFUNCTION(function, stdfunction)
Definition math.hh:214
Dune
Dune namespace
Definition alignedallocator.hh:13
Dune::binomial
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition math.hh:113
Dune::power
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition math.hh:75
Dune::sign
constexpr int sign(const T &val)
Return the sign of the value.
Definition math.hh:162
Dune::factorial
static constexpr T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition math.hh:93
Dune::conjugateComplex
constexpr K conjugateComplex(const K &x)
compute conjugate complex of x
Definition math.hh:146
Dune::MathOverloads
namespace for customization of math functions with Dune-Semantics
Definition fvector.hh:431
Dune::MathOverloads::isNaN
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition fvector.hh:458
Dune::MathOverloads::isInf
bool isInf(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is infinite.
Definition fvector.hh:446
Dune::MathOverloads::isFinite
auto isFinite(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether all entries are finite.
Definition fvector.hh:435
Dune::MathOverloads::isUnordered
bool isUnordered(const FieldVector< K, 1 > &b, const FieldVector< K, 1 > &c, PriorityTag< 2 >, ADLTag)
Returns true if either b or c is NaN.
Definition fvector.hh:470
Dune::StandardMathematicalConstants
Standard implementation of MathematicalConstants.
Definition math.hh:34
Dune::StandardMathematicalConstants::e
static const T e()
Euler's number.
Definition math.hh:38
Dune::StandardMathematicalConstants::pi
static const T pi()
Archimedes' constant.
Definition math.hh:48
Dune::MathematicalConstants
Provides commonly used mathematical constants.
Definition math.hh:67
Dune::MathOverloads::ADLTag
Tag to make sure the functions in this namespace can be found by ADL.
Definition math.hh:212
Dune::PriorityTag
Helper class for tagging priorities.
Definition typeutilities.hh:73