Dune Core Modules (2.7.0)

math.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_MATH_HH
4 #define DUNE_MATH_HH
5 
10 #include <cmath>
11 #include <complex>
12 #include <limits>
13 #include <type_traits>
14 
16 
17 namespace Dune
18 {
19 
30  template< class T >
32  {
36  static const T e ()
37  {
38  using std::exp;
39  static const T e = exp( T( 1 ) );
40  return e;
41  }
42 
46  static const T pi ()
47  {
48  using std::acos;
49  static const T pi = acos( T( -1 ) );
50  return pi;
51  }
52  };
53 
54 
62  template< class Field >
64  : public StandardMathematicalConstants<Field>
65  {};
66 
67 
72  template <class Mantissa, class Exponent>
73  constexpr Mantissa power(Mantissa m, Exponent p)
74  {
75  static_assert(std::numeric_limits<Exponent>::is_integer, "Exponent must be an integer type!");
76 
77  auto result = Mantissa(1);
78  auto absp = (p<0) ? -p : p; // This is simply abs, but std::abs is not constexpr
79  for (Exponent i = Exponent(0); i<absp; i++)
80  result *= m;
81 
82  if (p<0)
83  result = Mantissa(1)/result;
84 
85  return result;
86  }
87 
89  template <int m>
90  struct Factorial
91  {
93  enum { factorial = m * Factorial<m-1>::factorial };
94  };
95 
97  template <>
98  struct Factorial<0>
99  {
100  // 0! = 1
101  enum { factorial = 1 };
102  };
103 
104 
106  // T has to be an integral type
107  template<class T>
108  constexpr inline static T factorial(const T& n) noexcept
109  {
110  static_assert(std::numeric_limits<T>::is_integer, "`factorial(n)` has to be called with an integer type.");
111  T fac = 1;
112  for(T k = 0; k < n; ++k)
113  fac *= k+1;
114  return fac;
115  }
116 
118  template<class T, T n>
119  constexpr inline static auto factorial (std::integral_constant<T, n>) noexcept
120  {
121  return std::integral_constant<T, factorial(n)>{};
122  }
123 
124 
126  // T has to be an integral type
127  template<class T>
128  constexpr inline static T binomial (const T& n, const T& k) noexcept
129  {
130  static_assert(std::numeric_limits<T>::is_integer, "`binomial(n, k)` has to be called with an integer type.");
131 
132  if( k < 0 || k > n )
133  return 0;
134 
135  if (2*k > n)
136  return binomial(n, n-k);
137 
138  T bin = 1;
139  for(auto i = n-k; i < n; ++i)
140  bin *= i+1;
141  return bin / factorial(k);
142  }
143 
145  template<class T, T n, T k>
146  constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, k>) noexcept
147  {
148  return std::integral_constant<T, binomial(n, k)>{};
149  }
150 
151  template<class T, T n>
152  constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, n>) noexcept
153  {
154  return std::integral_constant<T, (n >= 0 ? 1 : 0)>{};
155  }
156 
157 
159  // conjugate complex does nothing for non-complex types
160  template<class K>
161  inline K conjugateComplex (const K& x)
162  {
163  return x;
164  }
165 
166 #ifndef DOXYGEN
167  // specialization for complex
168  template<class K>
169  inline std::complex<K> conjugateComplex (const std::complex<K>& c)
170  {
171  return std::complex<K>(c.real(),-c.imag());
172  }
173 #endif
174 
176  template <class T>
177  int sign(const T& val)
178  {
179  return (val < 0 ? -1 : 1);
180  }
181 
182 
183  namespace Impl {
184  // Returns whether a given type behaves like std::complex<>, i.e. whether
185  // real() and imag() are defined
186  template<class T>
187  struct isComplexLike {
188  private:
189  template<class U>
190  static auto test(U* u) -> decltype(u->real(), u->imag(), std::true_type());
191 
192  template<class U>
193  static auto test(...) -> decltype(std::false_type());
194 
195  public:
196  static const bool value = decltype(test<T>(0))::value;
197  };
198  } // namespace Impl
199 
201 
224  namespace MathOverloads {
225 
227  struct ADLTag {};
228 
229 #define DUNE_COMMON_MATH_ISFUNCTION(function, stdfunction) \
230  template<class T> \
231  auto function(const T &t, PriorityTag<1>, ADLTag) \
232  -> decltype(function(t)) { \
233  return function(t); \
234  } \
235  template<class T> \
236  auto function(const T &t, PriorityTag<0>, ADLTag) { \
237  using std::stdfunction; \
238  return stdfunction(t); \
239  } \
240  static_assert(true, "Require semicolon to unconfuse editors")
241 
242  DUNE_COMMON_MATH_ISFUNCTION(isNaN,isnan);
243  DUNE_COMMON_MATH_ISFUNCTION(isInf,isinf);
244  DUNE_COMMON_MATH_ISFUNCTION(isFinite,isfinite);
245 #undef DUNE_COMMON_MATH_ISFUNCTION
246 
247  template<class T>
248  auto isUnordered(const T &t1, const T &t2, PriorityTag<1>, ADLTag)
249  -> decltype(isUnordered(t1, t2)) {
250  return isUnordered(t1, t2);
251  }
252 
253  template<class T>
254  auto isUnordered(const T &t1, const T &t2, PriorityTag<0>, ADLTag) {
255  using std::isunordered;
256  return isunordered(t1, t2);
257  }
258  }
259 
260  namespace MathImpl {
261 
262  // NOTE: it is important that these functors have names different from the
263  // names of the functions they are forwarding to. Otherwise the
264  // unqualified call would find the functor type, not a function, and ADL
265  // would never be attempted.
266 #define DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(function) \
267  struct function##Impl { \
268  template<class T> \
269  constexpr auto operator()(const T &t) const { \
270  return function(t, PriorityTag<10>{}, MathOverloads::ADLTag{}); \
271  } \
272  }; \
273  static_assert(true, "Require semicolon to unconfuse editors")
274 
275  DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isNaN);
276  DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isInf);
277  DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isFinite);
278 #undef DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR
279 
280  struct isUnorderedImpl {
281  template<class T>
282  constexpr auto operator()(const T &t1, const T &t2) const {
283  return isUnordered(t1, t2, PriorityTag<10>{}, MathOverloads::ADLTag{});
284  }
285  };
286 
287  } //MathImpl
288 
289 
290  namespace Impl {
291  /* This helper has a math functor as a static constexpr member. Doing
292  this as a static member of a template struct means we can do this
293  without violating the ODR or putting the definition into a seperate
294  compilation unit, while still still ensuring the functor is the same
295  lvalue across all compilation units.
296  */
297  template<class T>
298  struct MathDummy
299  {
300  static constexpr T value{};
301  };
302 
303  template<class T>
304  constexpr T MathDummy<T>::value;
305 
306  } //namespace Impl
307 
308  namespace {
309  /* Provide the math functors directly in the `Dune` namespace.
310 
311  This actually declares a different name in each translation unit, but
312  they all resolve to the same lvalue.
313  */
314 
316 
320  constexpr auto const &isNaN = Impl::MathDummy<MathImpl::isNaNImpl>::value;
321 
323 
327  constexpr auto const &isInf = Impl::MathDummy<MathImpl::isInfImpl>::value;
328 
330 
334  constexpr auto const &isFinite = Impl::MathDummy<MathImpl::isFiniteImpl>::value;
335 
337 
342  constexpr auto const &isUnordered = Impl::MathDummy<MathImpl::isUnorderedImpl>::value;
343  }
344 
345  namespace MathOverloads {
346  /*Overloads for complex types*/
347  template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
348  auto isNaN(const T &t, PriorityTag<2>, ADLTag) {
349  return Dune::isNaN(real(t)) || Dune::isNaN(imag(t));
350  }
351 
352  template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
353  auto isInf(const T &t, PriorityTag<2>, ADLTag) {
354  return Dune::isInf(real(t)) || Dune::isInf(imag(t));
355  }
356 
357  template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
358  auto isFinite(const T &t, PriorityTag<2>, ADLTag) {
359  return Dune::isFinite(real(t)) && Dune::isFinite(imag(t));
360  }
361  } //MathOverloads
362 }
363 
364 #endif // #ifndef DUNE_MATH_HH
Dune namespace.
Definition: alignedallocator.hh:14
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
constexpr static T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition: math.hh:108
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:177
constexpr static T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:128
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:161
Calculates the factorial of m at compile time.
Definition: math.hh:91
Tag to make sure the functions in this namespace can be found by ADL.
Definition: math.hh:227
Provides commonly used mathematical constants.
Definition: math.hh:65
Helper class for tagging priorities.
Definition: typeutilities.hh:85
Helper class for tagging priorities.
Definition: typeutilities.hh:71
Standard implementation of MathematicalConstants.
Definition: math.hh:32
static const T e()
Euler's number.
Definition: math.hh:36
static const T pi()
Archimedes' constant.
Definition: math.hh:46
Utilities for type computations, constraining overloads, ...
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (May 13, 22:30, 2024)