DUNE PDELab (git)

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// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_MATH_HH
6#define DUNE_MATH_HH
7
12#include <cmath>
13#include <complex>
14#include <limits>
15#include <type_traits>
16
18
19namespace Dune
20{
21
32 template< class T >
34 {
38 static const T e ()
39 {
40 using std::exp;
41 static const T e = exp( T( 1 ) );
42 return e;
43 }
44
48 static const T pi ()
49 {
50 using std::acos;
51 static const T pi = acos( T( -1 ) );
52 return pi;
53 }
54 };
55
56
64 template< class Field >
66 : public StandardMathematicalConstants<Field>
67 {};
68
69
74 template <class Base, class Exponent>
75 constexpr Base power(Base m, Exponent p)
76 {
77 static_assert(std::numeric_limits<Exponent>::is_integer, "Exponent must be an integer type!");
78
79 auto result = Base(1);
80 auto absp = (p<0) ? -p : p; // This is simply abs, but std::abs is not constexpr
81 for (Exponent i = Exponent(0); i<absp; i++)
82 result *= m;
83
84 if (p<0)
85 result = Base(1)/result;
86
87 return result;
88 }
89
92 template <int m>
93 struct Factorial
94 {
96 static constexpr int factorial = m * Factorial<m-1>::factorial;
97 };
98
100 template <>
101 struct [[deprecated("Use function factorial instead! Will be removed after Dune 2.9")]] Factorial<0>
102 {
103 // 0! = 1
104 static constexpr int factorial = 1;
105 };
106
107
109 // T has to be an integral type
110 template<class T>
111 constexpr inline static T factorial(const T& n) noexcept
112 {
113 static_assert(std::numeric_limits<T>::is_integer, "`factorial(n)` has to be called with an integer type.");
114 T fac = 1;
115 for(T k = 0; k < n; ++k)
116 fac *= k+1;
117 return fac;
118 }
119
121 template<class T, T n>
122 constexpr inline static auto factorial (std::integral_constant<T, n>) noexcept
123 {
124 return std::integral_constant<T, factorial(n)>{};
125 }
126
127
129 // T has to be an integral type
130 template<class T>
131 constexpr inline static T binomial (const T& n, const T& k) noexcept
132 {
133 static_assert(std::numeric_limits<T>::is_integer, "`binomial(n, k)` has to be called with an integer type.");
134
135 if( k < 0 || k > n )
136 return 0;
137
138 if (2*k > n)
139 return binomial(n, n-k);
140
141 T bin = 1;
142 for(auto i = n-k; i < n; ++i)
143 bin *= i+1;
144 return bin / factorial(k);
145 }
146
148 template<class T, T n, T k>
149 constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, k>) noexcept
150 {
151 return std::integral_constant<T, binomial(n, k)>{};
152 }
153
154 template<class T, T n>
155 constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, n>) noexcept
156 {
157 return std::integral_constant<T, (n >= 0 ? 1 : 0)>{};
158 }
159
160
162 // conjugate complex does nothing for non-complex types
163 template<class K>
164 inline K conjugateComplex (const K& x)
165 {
166 return x;
167 }
168
169#ifndef DOXYGEN
170 // specialization for complex
171 template<class K>
172 inline std::complex<K> conjugateComplex (const std::complex<K>& c)
173 {
174 return std::complex<K>(c.real(),-c.imag());
175 }
176#endif
177
179 template <class T>
180 int sign(const T& val)
181 {
182 return (val < 0 ? -1 : 1);
183 }
184
185
186 namespace Impl {
187 // Returns whether a given type behaves like std::complex<>, i.e. whether
188 // real() and imag() are defined
189 template<class T>
190 struct isComplexLike {
191 private:
192 template<class U>
193 static auto test(U* u) -> decltype(u->real(), u->imag(), std::true_type());
194
195 template<class U>
196 static auto test(...) -> decltype(std::false_type());
197
198 public:
199 static const bool value = decltype(test<T>(0))::value;
200 };
201 } // namespace Impl
202
204
227 namespace MathOverloads {
228
230 struct ADLTag {};
231
232#define DUNE_COMMON_MATH_ISFUNCTION(function, stdfunction) \
233 template<class T> \
234 auto function(const T &t, PriorityTag<1>, ADLTag) \
235 -> decltype(function(t)) { \
236 return function(t); \
237 } \
238 template<class T> \
239 auto function(const T &t, PriorityTag<0>, ADLTag) { \
240 using std::stdfunction; \
241 return stdfunction(t); \
242 } \
243 static_assert(true, "Require semicolon to unconfuse editors")
244
245 DUNE_COMMON_MATH_ISFUNCTION(isNaN,isnan);
246 DUNE_COMMON_MATH_ISFUNCTION(isInf,isinf);
247 DUNE_COMMON_MATH_ISFUNCTION(isFinite,isfinite);
248#undef DUNE_COMMON_MATH_ISFUNCTION
249
250 template<class T>
251 auto isUnordered(const T &t1, const T &t2, PriorityTag<1>, ADLTag)
252 -> decltype(isUnordered(t1, t2)) {
253 return isUnordered(t1, t2);
254 }
255
256 template<class T>
257 auto isUnordered(const T &t1, const T &t2, PriorityTag<0>, ADLTag) {
258 using std::isunordered;
259 return isunordered(t1, t2);
260 }
261 }
262
263 namespace MathImpl {
264
265 // NOTE: it is important that these functors have names different from the
266 // names of the functions they are forwarding to. Otherwise the
267 // unqualified call would find the functor type, not a function, and ADL
268 // would never be attempted.
269#define DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(function) \
270 struct function##Impl { \
271 template<class T> \
272 constexpr auto operator()(const T &t) const { \
273 return function(t, PriorityTag<10>{}, MathOverloads::ADLTag{}); \
274 } \
275 }; \
276 static_assert(true, "Require semicolon to unconfuse editors")
277
278 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isNaN);
279 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isInf);
280 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isFinite);
281#undef DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR
282
283 struct isUnorderedImpl {
284 template<class T>
285 constexpr auto operator()(const T &t1, const T &t2) const {
286 return isUnordered(t1, t2, PriorityTag<10>{}, MathOverloads::ADLTag{});
287 }
288 };
289
290 } //MathImpl
291
292
293 namespace Impl {
294 /* This helper has a math functor as a static constexpr member. Doing
295 this as a static member of a template struct means we can do this
296 without violating the ODR or putting the definition into a separate
297 compilation unit, while still still ensuring the functor is the same
298 lvalue across all compilation units.
299 */
300 template<class T>
301 struct MathDummy
302 {
303 static constexpr T value{};
304 };
305
306 template<class T>
307 constexpr T MathDummy<T>::value;
308
309 } //namespace Impl
310
311 namespace {
312 /* Provide the math functors directly in the `Dune` namespace.
313
314 This actually declares a different name in each translation unit, but
315 they all resolve to the same lvalue.
316 */
317
319
323 constexpr auto const &isNaN = Impl::MathDummy<MathImpl::isNaNImpl>::value;
324
326
330 constexpr auto const &isInf = Impl::MathDummy<MathImpl::isInfImpl>::value;
331
333
337 constexpr auto const &isFinite = Impl::MathDummy<MathImpl::isFiniteImpl>::value;
338
340
345 constexpr auto const &isUnordered = Impl::MathDummy<MathImpl::isUnorderedImpl>::value;
346 }
347
348 namespace MathOverloads {
349 /*Overloads for complex types*/
350 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
351 auto isNaN(const T &t, PriorityTag<2>, ADLTag) {
352 return Dune::isNaN(real(t)) || Dune::isNaN(imag(t));
353 }
354
355 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
356 auto isInf(const T &t, PriorityTag<2>, ADLTag) {
357 return Dune::isInf(real(t)) || Dune::isInf(imag(t));
358 }
359
360 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
361 auto isFinite(const T &t, PriorityTag<2>, ADLTag) {
362 return Dune::isFinite(real(t)) && Dune::isFinite(imag(t));
363 }
364 } //MathOverloads
365}
366
367#endif // #ifndef DUNE_MATH_HH
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:627
bool isInf(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is infinite.
Definition: fvector.hh:615
auto isFinite(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether all entries are finite.
Definition: fvector.hh:604
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:639
Dune namespace.
Definition: alignedallocator.hh:13
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:131
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
static constexpr T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition: math.hh:111
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:164
Definition: math.hh:94
static constexpr int factorial
factorial stores m!
Definition: math.hh:96
Tag to make sure the functions in this namespace can be found by ADL.
Definition: math.hh:230
Provides commonly used mathematical constants.
Definition: math.hh:67
Helper class for tagging priorities.
Definition: typeutilities.hh:87
Helper class for tagging priorities.
Definition: typeutilities.hh:73
Standard implementation of MathematicalConstants.
Definition: math.hh:34
static const T e()
Euler's number.
Definition: math.hh:38
static const T pi()
Archimedes' constant.
Definition: math.hh:48
Utilities for type computations, constraining overloads, ...
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)