DUNE PDELab (git)

polynomial.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
8#define DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
9
10#include <cmath>
11#include <initializer_list>
12#include <vector>
13
14
15#include <dune/common/hybridutilities.hh>
16
17namespace Dune {
18namespace Functions {
19
20namespace Impl {
21
22 // Compute coefficients of derivative of polynomial.
23 // Overload for std::vector
24 template<class K, class Allocator>
25 auto polynomialDerivativeCoefficients(const std::vector<K, Allocator>& coefficients) {
26 if (coefficients.size()==0)
27 return std::vector<K, Allocator>();
28 std::vector<K, Allocator> dpCoefficients(coefficients.size()-1);
29 for (size_t i=1; i<coefficients.size(); ++i)
30 dpCoefficients[i-1] = coefficients[i]*K(i);
31 return dpCoefficients;
32 }
33
34 // Compute coefficients of derivative of polynomial.
35 // Overload for std::array
36 template<class K, unsigned long n>
37 auto polynomialDerivativeCoefficients(const std::array<K, n>& coefficients) {
38 if constexpr (n==0)
39 return coefficients;
40 else
41 {
42 std::array<K, n-1> dpCoefficients;
43 for (size_t i=1; i<coefficients.size(); ++i)
44 dpCoefficients[i-1] = coefficients[i]*K(i);
45 return dpCoefficients;
46 }
47 }
48
49 // Compute coefficients of derivative of polynomial.
50 // Helper function for the std::integer_sequence overload.
51 // With C++20 this can be avoided, because lambda function
52 // can partially specify template arguments which allows
53 // to do the same inline.
54 template<class I, I i0, I... i, class J, J j0, J... j>
55 auto polynomialDerivativeCoefficientsHelper(std::integer_sequence<I, i0, i...>, std::integer_sequence<J, j0, j...>) {
56 return std::integer_sequence<I, i*I(j)...>();
57 }
58
59 // Compute coefficients of derivative of polynomial.
60 // Overload for std::integer_sequence
61 template<class I, I... i>
62 auto polynomialDerivativeCoefficients(std::integer_sequence<I, i...> coefficients) {
63 if constexpr (sizeof...(i)==0)
64 return coefficients;
65 else
66 return polynomialDerivativeCoefficientsHelper(coefficients, std::make_index_sequence<sizeof...(i)>());
67 }
68
69 // Compute coefficients of derivative of polynomial.
70 // Overload for std::tuple
71 template<class...T>
72 auto polynomialDerivativeCoefficients(const std::tuple<T...>& coefficients) {
73 if constexpr (sizeof...(T)==0)
74 return coefficients;
75 else
76 {
77 // Notice that std::multiplies<void> has issues with signed types.
78 // E.g., `decltype(-2,2ul)` is `long unsigned int`.
79 // Hence the same is deduced as return type in std::multiplies.
80 // To avoid this, we explicitly pass the exponent `i+1` as signed type.
81 // If the coefficient is signed, both types are now signed and
82 // so is the deduced result type of std::multiplies.
83 auto mult = Dune::Hybrid::hybridFunctor(std::multiplies());
84 return Dune::unpackIntegerSequence([&](auto... i) {
85 return std::tuple(mult(std::get<i+1>(coefficients), std::integral_constant<long signed int, i+1>()) ...);
86 }, std::make_index_sequence<sizeof...(T)-1>());
87 }
88 }
89
90} // namespace Impl in Dune::Functions::
91
92
93
120template<class K, class C=std::vector<K>>
122{
123 template<class CC>
124 struct IsIntegerSequence : public std::false_type {};
125
126 template<class I, I... i>
127 struct IsIntegerSequence<std::integer_sequence<I, i...>> : public std::true_type {};
128
129public:
130
132 using Coefficients = C;
133
135 Polynomial() = default;
136
146 coefficients_(std::move(coefficients))
147 {}
148
150 K operator() (const K& x) const
151 {
152 auto y = K(0);
153 auto n = Dune::Hybrid::size(coefficients_);
154 Dune::Hybrid::forEach(Dune::range(n), [&](auto i) {
155 y += Dune::Hybrid::elementAt(coefficients_, i) * std::pow(x, int(i));
156 });
157 return y;
158 }
159
161 bool operator==(const Polynomial& other) const
162 {
163 if constexpr (IsIntegerSequence<Coefficients>::value)
164 return true;
165 else
166 return coefficients()==other.coefficients();
167 }
168
178 friend auto derivative(const Polynomial& p)
179 {
180 auto derivativeCoefficients = Impl::polynomialDerivativeCoefficients(p.coefficients());
181 using DerivativeCoefficients = decltype(derivativeCoefficients);
182 return Polynomial<K, DerivativeCoefficients>(std::move(derivativeCoefficients));
183 }
184
187 {
188 return coefficients_;
189 }
190
191private:
192 Coefficients coefficients_;
193};
194
195
196
197template<class K>
198Polynomial(std::vector<K>) -> Polynomial<K, std::vector<K>>;
199
200template<class K, unsigned long n>
201Polynomial(std::array<K,n>) -> Polynomial<K, std::array<K,n>>;
202
203template<class K, K... ci>
204Polynomial(std::integer_sequence<K, ci...>) -> Polynomial<K, std::integer_sequence<K,ci...>>;
205
206template<class K>
207Polynomial(std::initializer_list<K>) -> Polynomial<K, std::vector<K>>;
208
209
210
223template<class K, class Coefficients>
224auto makePolynomial(Coefficients coefficients)
225{
226 return Polynomial<K, Coefficients>(std::move(coefficients));
227}
228
238template<class K, class C>
239auto makePolynomial(std::initializer_list<C> coefficients)
240{
241 return Polynomial<K>(std::move(coefficients));
242}
243
244
245
246
247
248}} // namespace Dune::Functions
249
250
251
252#endif // DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
A univariate polynomial implementation.
Definition: polynomial.hh:122
Polynomial()=default
Default constructor.
K operator()(const K &x) const
Evaluate polynomial.
Definition: polynomial.hh:150
C Coefficients
The type of the stored coefficient container.
Definition: polynomial.hh:132
const Coefficients & coefficients() const
Obtain reference to coefficient vector.
Definition: polynomial.hh:186
Polynomial(Coefficients coefficients)
Create from container of coefficients.
Definition: polynomial.hh:145
bool operator==(const Polynomial &other) const
Comparison of coefficients.
Definition: polynomial.hh:161
decltype(auto) constexpr unpackIntegerSequence(F &&f, std::integer_sequence< I, i... > sequence)
Unpack an std::integer_sequence<I,i...> to std::integral_constant<I,i>...
Definition: indices.hh:124
friend auto derivative(const Polynomial &p)
Obtain derivative of Polynomial function.
Definition: polynomial.hh:178
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr decltype(auto) elementAt(Container &&c, Index &&i)
Get element at given position from container.
Definition: hybridutilities.hh:126
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)