Synopsis

#include <boost/math/tools/polynomial.hpp>

namespace boost { namespace math {
namespace tools {

template <class T>
class polynomial
{
public:
   // typedefs:
   typedef typename std::vector<T>::value_type value_type;
   typedef typename std::vector<T>::size_type  size_type;

   // construct:
   polynomial(){}
   template <class U>
   polynomial(const U* data, unsigned order);
   template <class I>
   polynomial(I first, I last);
   template <class U>
   explicit polynomial(const U& point,
         typename boost::enable_if<boost::is_convertible<U, T> >::type* = 0);
   template <class Range>
   explicit polynomial(const Range& r,
         typename boost::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = 0); // C++14
   polynomial(std::initializer_list<T> l); // C++11

   polynomial(std::vector<T>&& p);

   // access:
   size_type size()const;
   size_type degree()const;
   value_type& operator[](size_type i);
   const value_type& operator[](size_type i)const;
   std::vector<T> const& data() const;
   std::vector<T>& data();
   explicit operator bool() const;

   polynomial<T> prime() const;

   polynomial<T> integrate() const;

   T operator()(T z) const;



   // modify:
   void set_zero();

   // operators:
   template <class U>
   polynomial& operator +=(const U& value);
   template <class U>
   polynomial& operator -=(const U& value);
   template <class U>
   polynomial& operator *=(const U& value);
   template <class U>
   polynomial& operator /=(const U& value);
   template <class U>
   polynomial& operator %=(const U& value);
   template <class U>
   polynomial& operator +=(const polynomial<U>& value);
   template <class U>
   polynomial& operator -=(const polynomial<U>& value);
   template <class U>
   polynomial& operator *=(const polynomial<U>& value);
   template <class U>
   polynomial& operator /=(const polynomial<U>& value);
   template <class U>
   polynomial& operator %=(const polynomial<U>& value);
};

template <class T>
polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b);

template <class T, class U>
polynomial<T> operator + (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator - (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator * (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator / (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator % (const polynomial<T>& a, const U& b);

template <class U, class T>
polynomial<T> operator + (const U& a, const polynomial<T>& b);
template <class U, class T>
polynomial<T> operator - (const U& a, const polynomial<T>& b);
template <class U, class T>
polynomial<T> operator * (const U& a, const polynomial<T>& b);

template <class T>
polynomial<T> operator - (const polynomial<T>& a);

template <class T>
polynomial<T> operator >>= (const U& a);
template <class T>
polynomial<T> operator >> (polynomial<T> const &a, const U& b);
template <class T>
polynomial<T> operator <<= (const U& a);
template <class T>
polynomial<T> operator << (polynomial<T> const &a, const U& b);

template <class T>
bool operator == (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator != (const polynomial<T> &a, const polynomial<T> &b);

template <class T>
polynomial<T> pow(polynomial<T> base, int exp);

template <class charT, class traits, class T>
std::basic_ostream<charT, traits>& operator <<
   (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly);

template <typename T>
std::pair< polynomial<T>, polynomial<T> >
quotient_remainder(const polynomial<T>& a, const polynomial<T>& b);

} //    namespace tools
}} //    namespace boost { namespace math

Description

This is a somewhat trivial class for polynomial manipulation.

See

Implementation is currently of the "naive" variety, with 𝑶(N²) multiplication, for example. This class should not be used in high-performance computing environments: it is intended for the simple manipulation of small polynomials, typically generated for special function approximation.

It does has division for polynomials over a field (here floating point, complex, etc) and over a unique factorization domain (integers). Division of polynomials over a field is compatible with Euclidean GCD.

Division of polynomials over a UFD is compatible with the subresultant algorithm for GCD (implemented as subresultant_gcd), but a serious word of warning is required: the intermediate value swell of that algorithm will cause single-precision integral types to overflow very easily. So although the algorithm will work on single-precision integral types, an overload of the gcd function is only provided for polynomials with multi-precision integral types, to prevent nasty surprises. This is done somewhat crudely by disabling the overload for non-POD integral types.

Advanced manipulations: the FFT, factorisation etc are not currently provided. Submissions for these are of course welcome :-)

Polynomial Arithmetic Examples

First include the essential polynomial header (and others) to make the example:

#include <boost/math/tools/polynomial.hpp>

and some using statements are convenient:

using std::string;
using std::exception;
using std::cout;
using std::abs;
using std::pair;

using namespace boost::math;
using namespace boost::math::tools; // for polynomial
using boost::lexical_cast;

Store the coefficients in a convenient way to access them, then create some polynomials using construction from an iterator range, and finally output in a 'pretty' formula format.

	Tip
	Although we might conventionally write a polynomial from left to right in descending order of degree, Boost.Math stores in ascending order of degree.

Read/write for humans:    3x^3 - 4x^2 - 6x + 10
Boost polynomial storage: [ 10, -6, -4, 3 ]

boost::array<double, 4> const d3a = {{10, -6, -4, 3}};
polynomial<double> const a(d3a.begin(), d3a.end());

// With C++11 and later, you can also use initializer_list construction.
polynomial<double> const b{{-2.0, 1.0}};

// formula_format() converts from Boost storage to human notation.
cout << "a = " << formula_format(a)
<< "\nb = " << formula_format(b) << "\n\n";

for output:

a = 3x^3 - 4x^2 - 6x + 10
b = x - 2

// Now we can do arithmetic with the usual infix operators: + - * / and %.
polynomial<double> s = a + b;
cout << "a + b = " << formula_format(s) << "\n";
polynomial<double> d = a - b;
cout << "a - b = " << formula_format(d) << "\n";
polynomial<double> p = a * b;
cout << "a * b = " << formula_format(p) << "\n";
polynomial<double> q = a / b;
cout << "a / b = " << formula_format(q) << "\n";
polynomial<double> r = a % b;
cout << "a % b = " << formula_format(r) << "\n";

for output

a + b = 3x^3 - 4x^2 - 5x + 8
a - b = 3x^3 - 4x^2 - 7x + 12
a * b = 3x^4 - 10x^3 + 2x^2 + 22x - 20
a / b = 3x^2 + 2x - 2
a % b = 6

Division, Quotient and Remainder

Division is a special case where you can calculate two for the price of one.

Actually, quotient and remainder are always calculated together due to the nature of the algorithm: the infix operators return one result and throw the other away.

If you are doing a lot of division and want both the quotient and remainder, then you don't want to do twice the work necessary.

In that case you can call the underlying function, quotient_remainder, to get both results together as a pair.

  pair< polynomial<double>, polynomial<double> > result;
  result = quotient_remainder(a, b);
// Reassure ourselves that the result is the same.
  BOOST_ASSERT(result.first == q);
  BOOST_ASSERT(result.second == r);

The source code is at polynomial_arithmetic.cpp