Homec4science
Diffusion expolit (master) Browse

expolit/polynomial.hh
046a43293eecmaster

Last Change
Raw File

polynomial.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2018 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 Tamaas. 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 <iostream>
#include <numeric>
#include <ostream>
namespace expolit {
template <typename T, UInt order, typename CoeffType = T>
struct Polynomial : public Expression<Polynomial<T, order, CoeffType>> {
using value_type = T;
/// Number of terms
constexpr static UInt terms = order + 1;
constexpr Polynomial() = default;
constexpr explicit Polynomial(const std::array<CoeffType, terms>& coeffs)
: coeffs(coeffs) {}
/// Evaluating polynomial
constexpr auto operator()(const T& z) const {
T res(0);
for (UInt i = 0; i < terms; ++i) {
T prod(coeffs[i]);
for (UInt j = 0; j < i; ++j)
prod *= z;
res += prod;
}
return res;
}
/// Storing coefficients
std::array<CoeffType, terms> coeffs = {0};
};
template <typename T, typename CoeffType = T>
using Constant = Polynomial<T, 0, CoeffType>;
template <typename T, UInt order, typename CT, typename DT>
constexpr auto operator+(DT x, const Polynomial<T, order, CT>& p) {
Polynomial<T, order, CT> res(p);
res.coeffs.front() += x;
return res;
}
template <typename T, UInt order, typename CT>
constexpr auto operator+(const Polynomial<T, order, CT>& p, CT x) {
return x + p;
}
template <typename T, UInt order, typename CT>
constexpr auto operator*(CT x, const Polynomial<T, order, CT>& p) {
Polynomial<T, order, CT> res(p);
for (auto&& v : res.coeffs)
v *= x;
return res;
}
template <typename T, UInt order, typename CT>
constexpr auto operator*(const Polynomial<T, order, CT>& p, CT x) {
return x * p;
}
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1,
UInt o2>
constexpr auto operator*(const Polynomial<T1, o1, CT1>& p1,
const Polynomial<T2, o2, CT2>& p2) {
Polynomial<decltype(T1() * T2()), o1 + o2, decltype(CT1() * CT2())> res;
for (UInt i = 0; i < p1.terms; ++i)
for (UInt j = 0; j < p2.terms; ++j)
res.coeffs[i + j] += p1.coeffs[i] * p2.coeffs[j];
return res;
}
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1,
UInt o2>
constexpr auto operator+(const Polynomial<T1, o1, CT1>& p1,
const Polynomial<T2, o2, CT2>& p2) {
Polynomial<decltype(T1() + T2()), std::max(o1, o2), decltype(CT1() + CT2())>
res;
for (UInt i = 0; i < res.terms; ++i)
res.coeffs[i] += ((i < p1.terms) ? p1.coeffs[i] : 0) +
((i < p2.terms) ? p2.coeffs[i] : 0);
return res;
}
#define POLY_SIMPLIFICATION(opsign, std_op) \
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1, \
UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<T1, o1, CT1>& p1, \
const BinaryOperation<Polynomial<T2, o2, CT2>, Derived, std_op>& p2) { \
return p2.op(p2.op(p1, p2.operands.first), p2.operands.second); \
} \
\
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1, \
UInt o2, typename Derived> \
constexpr auto operator opsign( \
const BinaryOperation<Polynomial<T2, o2, CT2>, Derived, std_op>& p2, \
const Polynomial<T1, o1, CT1>& p1) { \
return p2.op(p1, p2); \
} \
\
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1, \
UInt o2, typename Derived> \
constexpr auto operator opsign( \
const Polynomial<T1, o1, CT1>& p1, \
const SBinaryOperation<Polynomial<T2, o2, CT2>, Derived, std_op>& p2) { \
return p2.op(p1, p2.commute()); \
} \
\
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1, \
UInt o2, typename Derived> \
constexpr auto operator opsign( \
const SBinaryOperation<Polynomial<T2, o2, CT2>, Derived, std_op>& p2, \
const Polynomial<T1, o1, CT1>& p1) { \
return p2.op(p1, p2.commute()); \
}
POLY_SIMPLIFICATION(+, std::plus<>);
POLY_SIMPLIFICATION(*, std::multiplies<>);
#undef POLY_SIMPLIFICATION
template <typename T, UInt order, typename CoeffType>
std::ostream& operator<<(std::ostream& os,
const Polynomial<T, order, CoeffType>& 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