Dune Core Modules (unstable)

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
91 // T has to be an integral type
92 template<class T>
93 constexpr inline static T factorial(const T& n) noexcept
94 {
95 static_assert(std::numeric_limits<T>::is_integer, "`factorial(n)` has to be called with an integer type.");
96 T fac = 1;
97 for(T k = 0; k < n; ++k)
98 fac *= k+1;
99 return fac;
100 }
101
103 template<class T, T n>
104 constexpr inline static auto factorial (std::integral_constant<T, n>) noexcept
105 {
106 return std::integral_constant<T, factorial(n)>{};
107 }
108
109
111 // T has to be an integral type
112 template<class T>
113 constexpr inline static T binomial (const T& n, const T& k) noexcept
114 {
115 static_assert(std::numeric_limits<T>::is_integer, "`binomial(n, k)` has to be called with an integer type.");
116
117 if( k < 0 || k > n )
118 return 0;
119
120 if (2*k > n)
121 return binomial(n, n-k);
122
123 T bin = 1;
124 for(auto i = n-k; i < n; ++i)
125 bin *= i+1;
126 return bin / factorial(k);
127 }
128
130 template<class T, T n, T k>
131 constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, k>) noexcept
132 {
133 return std::integral_constant<T, binomial(n, k)>{};
134 }
135
136 template<class T, T n>
137 constexpr inline static auto binomial (std::integral_constant<T, n>, std::integral_constant<T, n>) noexcept
138 {
139 return std::integral_constant<T, (n >= 0 ? 1 : 0)>{};
140 }
141
142
144 // conjugate complex does nothing for non-complex types
145 template<class K>
146 constexpr K conjugateComplex (const K& x)
147 {
148 return x;
149 }
150
151#ifndef DOXYGEN
152 // specialization for complex
153 template<class K>
154 constexpr std::complex<K> conjugateComplex (const std::complex<K>& c)
155 {
156 return std::complex<K>(c.real(),-c.imag());
157 }
158#endif
159
161 template <class T>
162 constexpr int sign(const T& val)
163 {
164 return (val < 0 ? -1 : 1);
165 }
166
167
168 namespace Impl {
169 // Returns whether a given type behaves like std::complex<>, i.e. whether
170 // real() and imag() are defined
171 template<class T>
172 struct isComplexLike {
173 private:
174 template<class U>
175 static auto test(U* u) -> decltype(u->real(), u->imag(), std::true_type());
176
177 template<class U>
178 static auto test(...) -> decltype(std::false_type());
179
180 public:
181 static const bool value = decltype(test<T>(0))::value;
182 };
183 } // namespace Impl
184
186
209 namespace MathOverloads {
210
212 struct ADLTag {};
213
214#define DUNE_COMMON_MATH_ISFUNCTION(function, stdfunction) \
215 template<class T> \
216 auto function(const T &t, PriorityTag<1>, ADLTag) \
217 -> decltype(function(t)) { \
218 return function(t); \
219 } \
220 template<class T> \
221 auto function(const T &t, PriorityTag<0>, ADLTag) { \
222 using std::stdfunction; \
223 return stdfunction(t); \
224 } \
225 static_assert(true, "Require semicolon to unconfuse editors")
226
227 DUNE_COMMON_MATH_ISFUNCTION(isNaN,isnan);
228 DUNE_COMMON_MATH_ISFUNCTION(isInf,isinf);
229 DUNE_COMMON_MATH_ISFUNCTION(isFinite,isfinite);
230#undef DUNE_COMMON_MATH_ISFUNCTION
231
232 template<class T>
233 auto isUnordered(const T &t1, const T &t2, PriorityTag<1>, ADLTag)
234 -> decltype(isUnordered(t1, t2)) {
235 return isUnordered(t1, t2);
236 }
237
238 template<class T>
239 auto isUnordered(const T &t1, const T &t2, PriorityTag<0>, ADLTag) {
240 using std::isunordered;
241 return isunordered(t1, t2);
242 }
243 }
244
245 namespace MathImpl {
246
247 // NOTE: it is important that these functors have names different from the
248 // names of the functions they are forwarding to. Otherwise the
249 // unqualified call would find the functor type, not a function, and ADL
250 // would never be attempted.
251#define DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(function) \
252 struct function##Impl { \
253 template<class T> \
254 constexpr auto operator()(const T &t) const { \
255 return function(t, PriorityTag<10>{}, MathOverloads::ADLTag{}); \
256 } \
257 }; \
258 static_assert(true, "Require semicolon to unconfuse editors")
259
260 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isNaN);
261 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isInf);
262 DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR(isFinite);
263#undef DUNE_COMMON_MATH_ISFUNCTION_FUNCTOR
264
265 struct isUnorderedImpl {
266 template<class T>
267 constexpr auto operator()(const T &t1, const T &t2) const {
268 return isUnordered(t1, t2, PriorityTag<10>{}, MathOverloads::ADLTag{});
269 }
270 };
271
272 } //MathImpl
273
274
275 namespace Impl {
276 /* This helper has a math functor as a static constexpr member. Doing
277 this as a static member of a template struct means we can do this
278 without violating the ODR or putting the definition into a separate
279 compilation unit, while still still ensuring the functor is the same
280 lvalue across all compilation units.
281 */
282 template<class T>
283 struct MathDummy
284 {
285 static constexpr T value{};
286 };
287
288 template<class T>
289 constexpr T MathDummy<T>::value;
290
291 } //namespace Impl
292
293 namespace {
294 /* Provide the math functors directly in the `Dune` namespace.
295
296 This actually declares a different name in each translation unit, but
297 they all resolve to the same lvalue.
298 */
299
301
305 constexpr auto const &isNaN = Impl::MathDummy<MathImpl::isNaNImpl>::value;
306
308
312 constexpr auto const &isInf = Impl::MathDummy<MathImpl::isInfImpl>::value;
313
315
319 constexpr auto const &isFinite = Impl::MathDummy<MathImpl::isFiniteImpl>::value;
320
322
327 constexpr auto const &isUnordered = Impl::MathDummy<MathImpl::isUnorderedImpl>::value;
328 }
329
330 namespace MathOverloads {
331 /*Overloads for complex types*/
332 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
333 auto isNaN(const T &t, PriorityTag<2>, ADLTag) {
334 return Dune::isNaN(real(t)) || Dune::isNaN(imag(t));
335 }
336
337 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
338 auto isInf(const T &t, PriorityTag<2>, ADLTag) {
339 return Dune::isInf(real(t)) || Dune::isInf(imag(t));
340 }
341
342 template<class T, class = std::enable_if_t<Impl::isComplexLike<T>::value> >
343 auto isFinite(const T &t, PriorityTag<2>, ADLTag) {
344 return Dune::isFinite(real(t)) && Dune::isFinite(imag(t));
345 }
346 } //MathOverloads
347}
348
349#endif // #ifndef DUNE_MATH_HH
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:458
bool isInf(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is infinite.
Definition: fvector.hh:446
auto isFinite(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether all entries are finite.
Definition: fvector.hh:435
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 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:113
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
constexpr int sign(const T &val)
Return the sign of the value.
Definition: math.hh:162
static constexpr T factorial(const T &n) noexcept
calculate the factorial of n as a constexpr
Definition: math.hh:93
constexpr K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:146
Tag to make sure the functions in this namespace can be found by ADL.
Definition: math.hh:212
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 & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 21, 22:40, 2025)