/usr/include/boost/math/tools
Edit: /usr/include/boost/math/tools/polynomial.hpp (21946B)
// (C) Copyright John Maddock 2006.
// (C) Copyright Jeremy William Murphy 2015.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP
#define BOOST_MATH_TOOLS_POLYNOMIAL_HPP
#ifdef _MSC_VER
#pragma once
#endif
#include
#include
#include
#ifdef BOOST_NO_CXX11_LAMBDAS
#include
#endif
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
#include
#endif
namespace boost{ namespace math{ namespace tools{
template
T chebyshev_coefficient(unsigned n, unsigned m)
{
BOOST_MATH_STD_USING
if(m > n)
return 0;
if((n & 1) != (m & 1))
return 0;
if(n == 0)
return 1;
T result = T(n) / 2;
unsigned r = n - m;
r /= 2;
BOOST_ASSERT(n - 2 * r == m);
if(r & 1)
result = -result;
result /= n - r;
result *= boost::math::binomial_coefficient(n - r, r);
result *= ldexp(1.0f, m);
return result;
}
template
Seq polynomial_to_chebyshev(const Seq& s)
{
// Converts a Polynomial into Chebyshev form:
typedef typename Seq::value_type value_type;
typedef typename Seq::difference_type difference_type;
Seq result(s);
difference_type order = s.size() - 1;
difference_type even_order = order & 1 ? order - 1 : order;
difference_type odd_order = order & 1 ? order : order - 1;
for(difference_type i = even_order; i >= 0; i -= 2)
{
value_type val = s[i];
for(difference_type k = even_order; k > i; k -= 2)
{
val -= result[k] * chebyshev_coefficient(static_cast(k), static_cast(i));
}
val /= chebyshev_coefficient(static_cast(i), static_cast(i));
result[i] = val;
}
result[0] *= 2;
for(difference_type i = odd_order; i >= 0; i -= 2)
{
value_type val = s[i];
for(difference_type k = odd_order; k > i; k -= 2)
{
val -= result[k] * chebyshev_coefficient(static_cast(k), static_cast(i));
}
val /= chebyshev_coefficient(static_cast(i), static_cast(i));
result[i] = val;
}
return result;
}
template
T evaluate_chebyshev(const Seq& a, const T& x)
{
// Clenshaw's formula:
typedef typename Seq::difference_type difference_type;
T yk2 = 0;
T yk1 = 0;
T yk = 0;
for(difference_type i = a.size() - 1; i >= 1; --i)
{
yk2 = yk1;
yk1 = yk;
yk = 2 * x * yk1 - yk2 + a[i];
}
return a[0] / 2 + yk * x - yk1;
}
template
class polynomial;
namespace detail {
/**
* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
* Chapter 4.6.1, Algorithm D: Division of polynomials over a field.
*
* @tparam T Coefficient type, must be not be an integer.
*
* Template-parameter T actually must be a field but we don't currently have that
* subtlety of distinction.
*/
template
BOOST_DEDUCED_TYPENAME disable_if_c::is_integer, void >::type
division_impl(polynomial &q, polynomial &u, const polynomial& v, N n, N k)
{
q[k] = u[n + k] / v[n];
for (N j = n + k; j > k;)
{
j--;
u[j] -= q[k] * v[j - k];
}
}
template
T integer_power(T t, N n)
{
switch(n)
{
case 0:
return static_cast(1u);
case 1:
return t;
case 2:
return t * t;
case 3:
return t * t * t;
}
T result = integer_power(t, n / 2);
result *= result;
if(n & 1)
result *= t;
return result;
}
/**
* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
* Chapter 4.6.1, Algorithm R: Pseudo-division of polynomials.
*
* @tparam T Coefficient type, must be an integer.
*
* Template-parameter T actually must be a unique factorization domain but we
* don't currently have that subtlety of distinction.
*/
template
BOOST_DEDUCED_TYPENAME enable_if_c::is_integer, void >::type
division_impl(polynomial &q, polynomial &u, const polynomial& v, N n, N k)
{
q[k] = u[n + k] * integer_power(v[n], k);
for (N j = n + k; j > 0;)
{
j--;
u[j] = v[n] * u[j] - (j < k ? T(0) : u[n + k] * v[j - k]);
}
}
/**
* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
* Chapter 4.6.1, Algorithm D and R: Main loop.
*
* @param u Dividend.
* @param v Divisor.
*/
template
std::pair< polynomial, polynomial >
division(polynomial u, const polynomial& v)
{
BOOST_ASSERT(v.size() <= u.size());
BOOST_ASSERT(v);
BOOST_ASSERT(u);
typedef typename polynomial::size_type N;
N const m = u.size() - 1, n = v.size() - 1;
N k = m - n;
polynomial q;
q.data().resize(m - n + 1);
do
{
division_impl(q, u, v, n, k);
}
while (k-- != 0);
u.data().resize(n);
u.normalize(); // Occasionally, the remainder is zeroes.
return std::make_pair(q, u);
}
//
// These structures are the same as the void specializations of the functors of the same name
// in the std lib from C++14 onwards:
//
struct negate
{
template
T operator()(T const &x) const
{
return -x;
}
};
struct plus
{
template
T operator()(T const &x, U const& y) const
{
return x + y;
}
};
struct minus
{
template
T operator()(T const &x, U const& y) const
{
return x - y;
}
};
} // namespace detail
/**
* Returns the zero element for multiplication of polynomials.
*/
template
polynomial zero_element(std::multiplies< polynomial >)
{
return polynomial();
}
template
polynomial identity_element(std::multiplies< polynomial >)
{
return polynomial(T(1));
}
/* Calculates a / b and a % b, returning the pair (quotient, remainder) together
* because the same amount of computation yields both.
* This function is not defined for division by zero: user beware.
*/
template
std::pair< polynomial, polynomial >
quotient_remainder(const polynomial& dividend, const polynomial& divisor)
{
BOOST_ASSERT(divisor);
if (dividend.size() < divisor.size())
return std::make_pair(polynomial(), dividend);
return detail::division(dividend, divisor);
}
template
class polynomial
{
public:
// typedefs:
typedef typename std::vector::value_type value_type;
typedef typename std::vector::size_type size_type;
// construct:
polynomial(){}
template
polynomial(const U* data, unsigned order)
: m_data(data, data + order + 1)
{
normalize();
}
template
polynomial(I first, I last)
: m_data(first, last)
{
normalize();
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
polynomial(std::vector&& p) : m_data(std::move(p))
{
normalize();
}
#endif
template
explicit polynomial(const U& point, typename boost::enable_if >::type* = 0)
{
if (point != U(0))
m_data.push_back(point);
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
// move:
polynomial(polynomial&& p) BOOST_NOEXCEPT
: m_data(std::move(p.m_data)) { }
#endif
// copy:
polynomial(const polynomial& p)
: m_data(p.m_data) { }
template
polynomial(const polynomial& p)
{
m_data.resize(p.size());
for(unsigned i = 0; i < p.size(); ++i)
{
m_data[i] = boost::math::tools::real_cast(p[i]);
}
}
#ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
template
explicit polynomial(const Range& r, typename boost::enable_if >::type* = 0)
: polynomial(r.begin(), r.end())
{
}
#endif
#if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
polynomial(std::initializer_list l) : polynomial(std::begin(l), std::end(l))
{
}
polynomial&
operator=(std::initializer_list l)
{
m_data.assign(std::begin(l), std::end(l));
normalize();
return *this;
}
#endif
// access:
size_type size() const { return m_data.size(); }
size_type degree() const
{
if (size() == 0)
throw std::logic_error("degree() is undefined for the zero polynomial.");
return m_data.size() - 1;
}
value_type& operator[](size_type i)
{
return m_data[i];
}
const value_type& operator[](size_type i) const
{
return m_data[i];
}
T evaluate(T z) const
{
return this->operator()(z);
}
T operator()(T z) const
{
return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : T(0);
}
std::vector chebyshev() const
{
return polynomial_to_chebyshev(m_data);
}
std::vector const& data() const
{
return m_data;
}
std::vector & data()
{
return m_data;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
polynomial prime() const
{
#ifdef BOOST_MSVC
// Disable int->float conversion warning:
#pragma warning(push)
#pragma warning(disable:4244)
#endif
if (m_data.size() == 0)
{
return polynomial({});
}
std::vector p_data(m_data.size() - 1);
for (size_t i = 0; i < p_data.size(); ++i) {
p_data[i] = m_data[i+1]*static_cast(i+1);
}
return polynomial(std::move(p_data));
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
}
polynomial integrate() const
{
std::vector i_data(m_data.size() + 1);
// Choose integration constant such that P(0) = 0.
i_data[0] = T(0);
for (size_t i = 1; i < i_data.size(); ++i)
{
i_data[i] = m_data[i-1]/static_cast(i);
}
return polynomial(std::move(i_data));
}
// operators:
polynomial& operator =(polynomial&& p) BOOST_NOEXCEPT
{
m_data = std::move(p.m_data);
return *this;
}
#endif
polynomial& operator =(const polynomial& p)
{
m_data = p.m_data;
return *this;
}
template
typename boost::enable_if_c::value, polynomial&>::type operator +=(const U& value)
{
addition(value);
normalize();
return *this;
}
template
typename boost::enable_if_c::value, polynomial&>::type operator -=(const U& value)
{
subtraction(value);
normalize();
return *this;
}
template
typename boost::enable_if_c::value, polynomial&>::type operator *=(const U& value)
{
multiplication(value);
normalize();
return *this;
}
template
typename boost::enable_if_c::value, polynomial&>::type operator /=(const U& value)
{
division(value);
normalize();
return *this;
}
template
typename boost::enable_if_c::value, polynomial&>::type operator %=(const U& /*value*/)
{
// We can always divide by a scalar, so there is no remainder:
this->set_zero();
return *this;
}
template
polynomial& operator +=(const polynomial& value)
{
addition(value);
normalize();
return *this;
}
template
polynomial& operator -=(const polynomial& value)
{
subtraction(value);
normalize();
return *this;
}
template
void multiply(const polynomial& a, const polynomial& b) {
if (!a || !b)
{
this->set_zero();
return;
}
std::vector prod(a.size() + b.size() - 1, T(0));
for (unsigned i = 0; i < a.size(); ++i)
for (unsigned j = 0; j < b.size(); ++j)
prod[i+j] += a.m_data[i] * b.m_data[j];
m_data.swap(prod);
}
template
polynomial& operator *=(const polynomial& value)
{
this->multiply(*this, value);
return *this;
}
template
polynomial& operator /=(const polynomial& value)
{
*this = quotient_remainder(*this, value).first;
return *this;
}
template
polynomial& operator %=(const polynomial& value)
{
*this = quotient_remainder(*this, value).second;
return *this;
}
template
polynomial& operator >>=(U const &n)
{
BOOST_ASSERT(n <= m_data.size());
m_data.erase(m_data.begin(), m_data.begin() + n);
return *this;
}
template
polynomial& operator <<=(U const &n)
{
m_data.insert(m_data.begin(), n, static_cast(0));
normalize();
return *this;
}
// Convenient and efficient query for zero.
bool is_zero() const
{
return m_data.empty();
}
// Conversion to bool.
#ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
typedef bool (polynomial::*unmentionable_type)() const;
BOOST_FORCEINLINE operator unmentionable_type() const
{
return is_zero() ? false : &polynomial::is_zero;
}
#else
BOOST_FORCEINLINE explicit operator bool() const
{
return !m_data.empty();
}
#endif
// Fast way to set a polynomial to zero.
void set_zero()
{
m_data.clear();
}
/** Remove zero coefficients 'from the top', that is for which there are no
* non-zero coefficients of higher degree. */
void normalize()
{
#ifndef BOOST_NO_CXX11_LAMBDAS
m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end());
#else
using namespace boost::lambda;
m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
#endif
}
private:
template
polynomial& addition(const U& value, R op)
{
if(m_data.size() == 0)
m_data.resize(1, 0);
m_data[0] = op(m_data[0], value);
return *this;
}
template
polynomial& addition(const U& value)
{
return addition(value, detail::plus());
}
template
polynomial& subtraction(const U& value)
{
return addition(value, detail::minus());
}
template
polynomial& addition(const polynomial& value, R op)
{
if (m_data.size() < value.size())
m_data.resize(value.size(), 0);
for(size_type i = 0; i < value.size(); ++i)
m_data[i] = op(m_data[i], value[i]);
return *this;
}
template
polynomial& addition(const polynomial& value)
{
return addition(value, detail::plus());
}
template
polynomial& subtraction(const polynomial& value)
{
return addition(value, detail::minus());
}
template
polynomial& multiplication(const U& value)
{
#ifndef BOOST_NO_CXX11_LAMBDAS
std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x * value; });
#else
using namespace boost::lambda;
std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret(_1 * value));
#endif
return *this;
}
template
polynomial& division(const U& value)
{
#ifndef BOOST_NO_CXX11_LAMBDAS
std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x / value; });
#else
using namespace boost::lambda;
std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret(_1 / value));
#endif
return *this;
}
std::vector m_data;
};
template
inline polynomial operator + (const polynomial& a, const polynomial& b)
{
polynomial result(a);
result += b;
return result;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template
inline polynomial operator + (polynomial&& a, const polynomial& b)
{
a += b;
return a;
}
template
inline polynomial operator + (const polynomial& a, polynomial&& b)
{
b += a;
return b;
}
template
inline polynomial operator + (polynomial&& a, polynomial&& b)
{
a += b;
return a;
}
#endif
template
inline polynomial operator - (const polynomial& a, const polynomial& b)
{
polynomial result(a);
result -= b;
return result;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template
inline polynomial operator - (polynomial&& a, const polynomial& b)
{
a -= b;
return a;
}
template
inline polynomial operator - (const polynomial& a, polynomial&& b)
{
b -= a;
return -b;
}
template
inline polynomial operator - (polynomial&& a, polynomial&& b)
{
a -= b;
return a;
}
#endif
template
inline polynomial operator * (const polynomial& a, const polynomial& b)
{
polynomial result;
result.multiply(a, b);
return result;
}
template
inline polynomial operator / (const polynomial& a, const polynomial& b)
{
return quotient_remainder(a, b).first;
}
template
inline polynomial operator % (const polynomial& a, const polynomial& b)
{
return quotient_remainder(a, b).second;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator + (polynomial a, const U& b)
{
a += b;
return a;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator - (polynomial a, const U& b)
{
a -= b;
return a;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator * (polynomial a, const U& b)
{
a *= b;
return a;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator / (polynomial a, const U& b)
{
a /= b;
return a;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator % (const polynomial&, const U&)
{
// Since we can always divide by a scalar, result is always an empty polynomial:
return polynomial();
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator + (const U& a, polynomial b)
{
b += a;
return b;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator - (const U& a, polynomial b)
{
b -= a;
return -b;
}
template
inline typename boost::enable_if_c::value, polynomial >::type operator * (const U& a, polynomial b)
{
b *= a;
return b;
}
template
bool operator == (const polynomial &a, const polynomial &b)
{
return a.data() == b.data();
}
template
bool operator != (const polynomial &a, const polynomial &b)
{
return a.data() != b.data();
}
template
polynomial operator >> (polynomial a, const U& b)
{
a >>= b;
return a;
}
template
polynomial operator << (polynomial a, const U& b)
{
a <<= b;
return a;
}
// Unary minus (negate).
template
polynomial operator - (polynomial a)
{
std::transform(a.data().begin(), a.data().end(), a.data().begin(), detail::negate());
return a;
}
template
bool odd(polynomial const &a)
{
return a.size() > 0 && a[0] != static_cast(0);
}
template
bool even(polynomial const &a)
{
return !odd(a);
}
template
polynomial pow(polynomial base, int exp)
{
if (exp < 0)
return policies::raise_domain_error(
"boost::math::tools::pow<%1%>",
"Negative powers are not supported for polynomials.",
base, policies::policy<>());
// if the policy is ignore_error or errno_on_error, raise_domain_error
// will return std::numeric_limits>::quiet_NaN(), which
// defaults to polynomial(), which is the zero polynomial
polynomial result(T(1));
if (exp & 1)
result = base;
/* "Exponentiation by squaring" */
while (exp >>= 1)
{
base *= base;
if (exp & 1)
result *= base;
}
return result;
}
template
inline std::basic_ostream& operator << (std::basic_ostream& os, const polynomial& poly)
{
os << "{ ";
for(unsigned i = 0; i < poly.size(); ++i)
{
if(i) os << ", ";
os << poly[i];
}
os << " }";
return os;
}
} // namespace tools
} // namespace math
} // namespace boost
//
// Polynomial specific overload of gcd algorithm:
//
#include
#endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP