Loading...
Searching...
No Matches
polynomial.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-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
#define DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
#include <cmath>
#include <initializer_list>
#include <vector>
#include <dune/common/hybridutilities.hh>
namespace Dune {
namespace Functions {
namespace Impl {
// Compute coefficients of derivative of polynomial.
// Overload for std::vector
template<class K, class Allocator>
auto polynomialDerivativeCoefficients(const std::vector<K, Allocator>& coefficients) {
if (coefficients.size()==0)
return std::vector<K, Allocator>();
std::vector<K, Allocator> dpCoefficients(coefficients.size()-1);
for (size_t i=1; i<coefficients.size(); ++i)
dpCoefficients[i-1] = coefficients[i]*K(i);
return dpCoefficients;
}
// Compute coefficients of derivative of polynomial.
// Overload for std::array
template<class K, std::size_t n>
auto polynomialDerivativeCoefficients(const std::array<K, n>& coefficients) {
if constexpr (n==0)
return coefficients;
else
{
std::array<K, n-1> dpCoefficients;
for (size_t i=1; i<coefficients.size(); ++i)
dpCoefficients[i-1] = coefficients[i]*K(i);
return dpCoefficients;
}
}
// Compute coefficients of derivative of polynomial.
// Helper function for the std::integer_sequence overload.
// With C++20 this can be avoided, because lambda function
// can partially specify template arguments which allows
// to do the same inline.
template<class I, I i0, I... i, class J, J j0, J... j>
auto polynomialDerivativeCoefficientsHelper(std::integer_sequence<I, i0, i...>, std::integer_sequence<J, j0, j...>) {
return std::integer_sequence<I, i*I(j)...>();
}
// Compute coefficients of derivative of polynomial.
// Overload for std::integer_sequence
template<class I, I... i>
auto polynomialDerivativeCoefficients(std::integer_sequence<I, i...> coefficients) {
if constexpr (sizeof...(i)==0)
return coefficients;
else
return polynomialDerivativeCoefficientsHelper(coefficients, std::make_index_sequence<sizeof...(i)>());
}
// Compute coefficients of derivative of polynomial.
// Overload for std::tuple
template<class...T>
auto polynomialDerivativeCoefficients(const std::tuple<T...>& coefficients) {
if constexpr (sizeof...(T)==0)
return coefficients;
else
{
// Notice that std::multiplies<void> has issues with signed types.
// E.g., `decltype(-2,2ul)` is `long unsigned int`.
// Hence the same is deduced as return type in std::multiplies.
// To avoid this, we explicitly pass the exponent `i+1` as signed type.
// If the coefficient is signed, both types are now signed and
// so is the deduced result type of std::multiplies.
auto mult = Dune::Hybrid::hybridFunctor(std::multiplies());
return Dune::unpackIntegerSequence([&](auto... i) {
using signed_type = std::conditional_t<std::is_same_v<std::size_t, long unsigned int>,
long signed int, signed int>;
return std::tuple(mult(std::get<i+1>(coefficients),
std::integral_constant<signed_type, (signed_type)(i+1)>()) ...);
}, std::make_index_sequence<sizeof...(T)-1>());
}
}
} // namespace Impl in Dune::Functions::
template<class K, class C=std::vector<K>>
class Polynomial
{
template<class CC>
struct IsIntegerSequence : public std::false_type {};
template<class I, I... i>
struct IsIntegerSequence<std::integer_sequence<I, i...>> : public std::true_type {};
template <typename Coefficient>
static void add(K& y, const Coefficient &c)
{
if constexpr (!IsIntegralConstant<Coefficient>::value)
{
if (c!=0)
y += c;
}
else
y += c;
}
public:
using Coefficients = C;
Polynomial() = default;
Polynomial(Coefficients coefficients) :
coefficients_(std::move(coefficients))
{}
K operator() (const K& x) const
{
using namespace Dune::Indices;
auto n = Dune::Hybrid::size(coefficients_);
// Explicitly handling the corner case of an empty coefficient set
// allows to save one multiplication.
return Hybrid::ifElse(Hybrid::equal_to(n, _0),
[&](auto id) { /* then */
// No coefficients at all
return K(0);
},
[&](auto id) { /* else */
// Do the Horner scheme knowing that there is at least one coefficient
K y = Hybrid::elementAt(coefficients_, Hybrid::minus(id(n), _1) );
Dune::Hybrid::forEach(Dune::range(Hybrid::minus(id(n), _1) ), [&](auto i)
{
y *= x;
// Do y+= _next coefficient_, unless that coefficient is std::integral_constant<0>.
add(y,Hybrid::elementAt(coefficients_, Hybrid::minus(Hybrid::minus(id(n),_2), i)));
});
return y;
});
}
bool operator==(const Polynomial& other) const
{
if constexpr (IsIntegerSequence<Coefficients>::value)
return true;
else
return coefficients()==other.coefficients();
}
friend auto derivative(const Polynomial& p)
{
auto derivativeCoefficients = Impl::polynomialDerivativeCoefficients(p.coefficients());
using DerivativeCoefficients = decltype(derivativeCoefficients);
return Polynomial<K, DerivativeCoefficients>(std::move(derivativeCoefficients));
}
const Coefficients& coefficients() const
{
return coefficients_;
}
private:
Coefficients coefficients_;
};
template<class K>
Polynomial(std::vector<K>) -> Polynomial<K, std::vector<K>>;
template<class K, std::size_t n>
Polynomial(std::array<K,n>) -> Polynomial<K, std::array<K,n>>;
template<class K, K... ci>
Polynomial(std::integer_sequence<K, ci...>) -> Polynomial<K, std::integer_sequence<K,ci...>>;
template<class K>
Polynomial(std::initializer_list<K>) -> Polynomial<K, std::vector<K>>;
template<class K, class Coefficients>
auto makePolynomial(Coefficients coefficients)
{
return Polynomial<K, Coefficients>(std::move(coefficients));
}
template<class K, class C>
auto makePolynomial(std::initializer_list<C> coefficients)
{
return Polynomial<K>(std::move(coefficients));
}
}} // namespace Dune::Functions
#endif // DUNE_FUNCTIONS_ANALYTICFUNCTIONS_POLYNOMIAL_HH
Dune::Functions::Polynomial::derivative
friend auto derivative(const Polynomial &p)
Obtain derivative of Polynomial function.
Definition polynomial.hh:222
Dune
Definition monomialset.hh:19
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::Polynomial
Polynomial(std::vector< K >) -> Polynomial< K, std::vector< K > >
Dune::Functions::makePolynomial
auto makePolynomial(Coefficients coefficients)
Create Polynomial.
Definition polynomial.hh:268
Dune::Functions::Polynomial
A univariate polynomial implementation.
Definition polynomial.hh:125
Dune::Functions::Polynomial::Polynomial
Polynomial()=default
Default constructor.
Dune::Functions::Polynomial::operator()
K operator()(const K &x) const
Evaluate polynomial using the Horner scheme.
Definition polynomial.hh:178
Dune::Functions::Polynomial::Coefficients
C Coefficients
The type of the stored coefficient container.
Definition polynomial.hh:156
Dune::Functions::Polynomial::coefficients
const Coefficients & coefficients() const
Obtain reference to coefficient vector.
Definition polynomial.hh:230
Dune::Functions::Polynomial::Polynomial
Polynomial(Coefficients coefficients)
Create from container of coefficients.
Definition polynomial.hh:169
Dune::Functions::Polynomial::operator==
bool operator==(const Polynomial &other) const
Comparison of coefficients.
Definition polynomial.hh:205