Page MenuHomec4science
Files F74454181

polynomial.hh
No OneTemporary
Actions

File Metadata

Created
Sat, Jul 27, 22:14

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.
*
* Tamaas 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 <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;
Polynomial() = default;
explicit Polynomial(const std::array<CoeffType, terms>& coeffs)
: coeffs(coeffs) {}
/// Evaluating polynomial
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 T1, typename T2, typename CT1, typename CT2, UInt o1,
UInt o2>
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, typename Derived, typename Op>
auto
operator*(const Polynomial<T1, o1, CT1>& p1,
const BinaryOperation<Polynomial<T2, o2, CT2>, Derived, 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, typename Op>
auto operator*(const BinaryOperation<Polynomial<T2, o2, CT2>, Derived, 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, typename Op>
auto
operator*(const Polynomial<T1, o1, CT1>& p1,
const SBinaryOperation<Polynomial<T2, o2, CT2>, Derived, Op>& p2) {
return p2.op(p1, p2.commute());
}
template <typename T1, typename T2, typename CT1, typename CT2, UInt o1,
UInt o2, typename Derived, typename Op>
auto operator*(const SBinaryOperation<Polynomial<T2, o2, CT2>, Derived, Op>& p2,
const Polynomial<T1, o1, CT1>& p1) {
return p2.op(p1, p2.commute());
}
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

Event Timeline

c4science · Help