Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F59596104
polynomial.hh
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Apr 23, 18:52
Size
11 KB
Mime Type
text/x-c++
Expires
Thu, Apr 25, 18:52 (2 d)
Engine
blob
Format
Raw Data
Handle
17217208
Attached To
rEXPOLIT expolit
polynomial.hh
View Options
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2018-2021 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
Event Timeline
Log In to Comment