Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/third-party/expolit/include/expolit/polynomial.hh
e0a01ad22e7emaxwell_viscoelastic…

Last Change
Raw File

polynomial.hh
View Options

/*
* Copyright (©) 2018-2022 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Expolit is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Expolit is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
* more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Expolit. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef POLYNOMIALS_HH
#define POLYNOMIALS_HH
#include "operations.hh"
#include "types.hh"
#include <algorithm>
#include <array>
#include <numeric>
#include <ostream>
namespace expolit {
template <typename From, typename To>
using compatible = std::enable_if_t<std::is_convertible<From, To>::value>;
/**
* @brief Polynomial type
*
* Type for polynomial expressions of arbitrary order, variable type and
* coefficient type
*/
template <typename CoeffType, UInt order>
struct Polynomial : Expression<Polynomial<CoeffType, order>> {
/// Number of terms
constexpr static std::size_t terms = order + 1;
constexpr Polynomial() = default;
constexpr explicit Polynomial(const std::array<CoeffType, terms>& coeffs)
: coeffs(coeffs) {}
/// Evaluating polynomial
template <typename Val, typename = compatible<Val, CoeffType>>
constexpr auto operator()(Val z) const {
using Res = decltype(std::declval<CoeffType>() * z);
Res res(coeffs.front());
Val power(z);
for (UInt i = 1; i < terms; ++i) {
res += coeffs[i] * power;
power *= z;
}
return res;
}
/// Storing coefficients
std::array<CoeffType, terms> coeffs = {0};
};
template <typename CoeffType>
using Constant = Polynomial<CoeffType, 0>;
namespace detail {
template <typename Poly>
struct get_terms
: std::integral_constant<
std::size_t, std::remove_reference_t<std::remove_cv_t<Poly>>::terms> {
};
} // namespace detail
/* -------------------------------------------------------------------------- */
/* Algebraic rules */
/* -------------------------------------------------------------------------- */
/// Sum with scalar
template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
constexpr auto operator+(DT x, const Polynomial<CT, order>& p) {
Polynomial<decltype(std::declval<CT>() + x), order> res(p);
std::get<0>(res.coeffs) += x;
return res;
}
/// Sum with scalar
template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
constexpr auto operator+(const Polynomial<CT, order>& p, DT x) {
return x + p;
}
namespace detail {
template <std::size_t i>
struct for_scalar {
template <typename CT, std::size_t n, typename DT>
constexpr static void prod(std::array<CT, n>& res, DT x) {
std::get<i>(res) *= x;
for_scalar<i - 1>::prod(res, x);
}
template <typename CT, std::size_t n, typename DT>
constexpr static void div(std::array<CT, n>& res, DT x) {
std::get<i>(res) /= x;
for_scalar<i - 1>::prod(res, x);
}
};
template <>
struct for_scalar<0> {
template <typename CT, std::size_t n, typename DT>
constexpr static void prod(std::array<CT, n>& res, DT x) {
std::get<0>(res) *= x;
}
template <typename CT, std::size_t n, typename DT>
constexpr static void div(std::array<CT, n>& res, DT x) {
std::get<0>(res) /= x;
}
};
} // namespace detail
/// Product with scalar
template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
constexpr auto operator*(DT x, const Polynomial<CT, order>& p) {
Polynomial<decltype(std::declval<CT>() * x), order> res(p);
detail::for_scalar<res.terms - 1>::prod(res.coeffs, x);
return res;
}
/// Product with scalar
template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
constexpr auto operator*(const Polynomial<CT, order>& p, DT x) {
return x * p;
}
/// Division with scalar
template <typename CT, UInt order, typename DT, typename = compatible<DT, CT>>
constexpr auto operator/(const Polynomial<CT, order>& p, DT x) {
Polynomial<decltype(std::declval<CT>() / x), order> res(p);
detail::for_scalar<res.terms - 1>::div(res.coeffs, x);
return res;
}
/// Division with constant
template <typename CT, UInt order, typename CT2, typename = compatible<CT2, CT>>
constexpr auto operator/(const Polynomial<CT, order>& p,
const Constant<CT2>& x) {
Polynomial<decltype(std::declval<CT>() / x.coeffs.front()), order> res(p);
detail::for_scalar<res.terms - 1>::div(res.coeffs, std::get<0>(x.coeffs));
return res;
}
namespace detail {
template <std::size_t i, std::size_t j>
struct coeff_product {
template <typename CT1, typename CT2, typename CT3, std::size_t n,
std::size_t m, std::size_t l>
static constexpr void prod(std::array<CT1, n>& res,
const std::array<CT2, m>& a,
const std::array<CT3, l>& b) {
std::get<i + j>(res) += a[i] * b[j];
coeff_product<i, j - 1>::prod(res, a, b);
}
};
template <std::size_t i>
struct coeff_product<i, 0> {
template <typename CT1, typename CT2, typename CT3, std::size_t n,
std::size_t m, std::size_t l>
static constexpr void prod(std::array<CT1, n>& res,
const std::array<CT2, m>& a,
const std::array<CT3, l>& b) {
std::get<i>(res) += a[i] * b[0];
coeff_product<i - 1, l - 1>::prod(res, a, b);
}
};
template <>
struct coeff_product<0, 0> {
template <typename CT1, typename CT2, typename CT3, std::size_t n,
std::size_t m, std::size_t l>
static constexpr void prod(std::array<CT1, n>& res,
const std::array<CT2, m>& a,
const std::array<CT3, l>& b) {
std::get<0>(res) += a[0] * b[0];
}
};
} // namespace detail
/// Product of polynomials
template <typename CT1, typename CT2, UInt o1, UInt o2>
constexpr auto operator*(const Polynomial<CT1, o1>& p1,
const Polynomial<CT2, o2>& p2) {
Polynomial<decltype(std::declval<CT1>() * std::declval<CT2>()), o1 + o2> res;
detail::coeff_product<detail::get_terms<decltype(p1)>::value - 1,
detail::get_terms<decltype(p2)>::value -
1>::prod(res.coeffs, p1.coeffs, p2.coeffs);
return res;
}
namespace detail {
template <std::size_t i>
struct coeff_sum {
template <typename CT1, typename CT2, typename CT3, std::size_t n,
std::size_t m, std::size_t l>
static constexpr void sum(std::array<CT1, n>& res,
const std::array<CT2, m>& a,
const std::array<CT3, l>& b) {
std::get<i>(res) += ((i < m) ? a[i] : 0) + ((i < l) ? b[i] : 0);
coeff_sum<i - 1>::sum(res, a, b);
}
};
template <>
struct coeff_sum<0> {
template <typename CT1, typename CT2, typename CT3, std::size_t n,
std::size_t m, std::size_t l>
static constexpr void sum(std::array<CT1, n>& res,
const std::array<CT2, m>& a,
const std::array<CT3, l>& b) {
std::get<0>(res) += a[0] + b[0];
}
};
} // namespace detail
/// Sum of polynomials
template <typename CT1, typename CT2, UInt o1, UInt o2>
constexpr auto operator+(const Polynomial<CT1, o1>& p1,
const Polynomial<CT2, o2>& p2) {
Polynomial<decltype(std::declval<CT1>() + std::declval<CT2>()),
std::max(o1, o2)>
res;
detail::coeff_sum<res.terms - 1>::sum(res.coeffs, p1.coeffs, p2.coeffs);
return res;
}
/* -------------------------------------------------------------------------- */
/* Simplification rules */
/* -------------------------------------------------------------------------- */
#define POLY_SIMPLIFICATION(opsign, std_op) \
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<CT1, o1>& p1, \
const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) { \
return p2.op(p2.op(p1, p2.operands.first), p2.operands.second); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const BinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2, \
const Polynomial<CT1, o1>& p1) { \
return p2.op(p1, p2); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<CT1, o1>& p1, \
const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2) { \
return p2.op(p1, p2.commute()); \
} \
\
template <typename CT1, typename CT2, UInt o1, UInt o2, typename Derived> \
constexpr auto operator opsign( \
const SBinaryOperation<Polynomial<CT2, o2>, Derived, std_op>& p2, \
const Polynomial<CT1, o1>& p1) { \
return p2.op(p1, p2.commute()); \
}
POLY_SIMPLIFICATION(+, std::plus<>)
POLY_SIMPLIFICATION(*, std::multiplies<>)
#undef POLY_SIMPLIFICATION
template <typename CT1, UInt o, typename CT2, typename Derived>
constexpr auto operator/(const Product<Polynomial<CT1, o>, Derived>& p,
const Constant<CT2>& c) {
return p.operands.second * (p.operands.first / c);
}
template <typename CT1, UInt o, typename CT2, typename Derived>
constexpr auto operator/(const Product<Derived, Polynomial<CT1, o>>& p,
const Constant<CT2>& c) {
return p.operands.second * (p.operands.first / c);
}
template <typename CT1, UInt o, typename CT2, typename Derived>
constexpr auto operator*(const Polynomial<CT1, o>& p,
const Division<Derived, Constant<CT2>>& d) {
return (p / d.operands.second) * d.operands.first;
}
template <typename CT1, UInt o, typename CT2, typename Derived>
constexpr auto operator*(const Division<Derived, Constant<CT2>>& d,
const Polynomial<CT1, o>& p) {
return (p / d.operands.second) * d.operands.first;
}
/* -------------------------------------------------------------------------- */
/* Output */
/* -------------------------------------------------------------------------- */
template <typename CT, UInt order>
std::ostream& operator<<(std::ostream& os, const Polynomial<CT, order>& p) {
os << p.coeffs[0];
if (p.terms >= 2) {
os << " + " << p.coeffs[1] << "*z";
for (UInt i = 2; i < p.terms; ++i)
os << " + " << p.coeffs[i] << "*z^" << i;
}
return os;
}
} // namespace expolit
#endif
c4science · Help