boost/math/tools/polynomial.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  (C) Copyright Jeremy William Murphy 2015.
   3
   4
   5 //  Use, modification and distribution are subject to the
   6 //  Boost Software License, Version 1.0. (See accompanying file
   7 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 #ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP
  10 #define BOOST_MATH_TOOLS_POLYNOMIAL_HPP
  11
  12 #ifdef _MSC_VER
  13 #pragma once
  14 #endif
  15
  16 #include <boost/assert.hpp>
  17 #include <boost/config.hpp>
  18 #ifdef BOOST_NO_CXX11_LAMBDAS
  19 #include <boost/lambda/lambda.hpp>
  20 #endif
  21 #include <boost/math/tools/rational.hpp>
  22 #include <boost/math/tools/real_cast.hpp>
  23 #include <boost/math/policies/error_handling.hpp>
  24 #include <boost/math/special_functions/binomial.hpp>
  25 #include <boost/core/enable_if.hpp>
  26 #include <boost/type_traits/is_convertible.hpp>
  27 #include <boost/math/tools/detail/is_const_iterable.hpp>
  28
  29 #include <vector>
  30 #include <ostream>
  31 #include <algorithm>
  32 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
  33 #include <initializer_list>
  34 #endif
  35
  36 namespace boost{ namespace math{ namespace tools{
  37
  38 template <class T>
  39 T chebyshev_coefficient(unsigned n, unsigned m)
  40 {
  41    BOOST_MATH_STD_USING
  42    if(m > n)
  43       return 0;
  44    if((n & 1) != (m & 1))
  45       return 0;
  46    if(n == 0)
  47       return 1;
  48    T result = T(n) / 2;
  49    unsigned r = n - m;
  50    r /= 2;
  51
  52    BOOST_ASSERT(n - 2 * r == m);
  53
  54    if(r & 1)
  55       result = -result;
  56    result /= n - r;
  57    result *= boost::math::binomial_coefficient<T>(n - r, r);
  58    result *= ldexp(1.0f, m);
  59    return result;
  60 }
  61
  62 template <class Seq>
  63 Seq polynomial_to_chebyshev(const Seq& s)
  64 {
  65    // Converts a Polynomial into Chebyshev form:
  66    typedef typename Seq::value_type value_type;
  67    typedef typename Seq::difference_type difference_type;
  68    Seq result(s);
  69    difference_type order = s.size() - 1;
  70    difference_type even_order = order & 1 ? order - 1 : order;
  71    difference_type odd_order = order & 1 ? order : order - 1;
  72
  73    for(difference_type i = even_order; i >= 0; i -= 2)
  74    {
  75       value_type val = s[i];
  76       for(difference_type k = even_order; k > i; k -= 2)
  77       {
  78          val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
  79       }
  80       val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
  81       result[i] = val;
  82    }
  83    result[0] *= 2;
  84
  85    for(difference_type i = odd_order; i >= 0; i -= 2)
  86    {
  87       value_type val = s[i];
  88       for(difference_type k = odd_order; k > i; k -= 2)
  89       {
  90          val -= result[k] * chebyshev_coefficient<value_type>(static_cast<unsigned>(k), static_cast<unsigned>(i));
  91       }
  92       val /= chebyshev_coefficient<value_type>(static_cast<unsigned>(i), static_cast<unsigned>(i));
  93       result[i] = val;
  94    }
  95    return result;
  96 }
  97
  98 template <class Seq, class T>
  99 T evaluate_chebyshev(const Seq& a, const T& x)
 100 {
 101    // Clenshaw's formula:
 102    typedef typename Seq::difference_type difference_type;
 103    T yk2 = 0;
 104    T yk1 = 0;
 105    T yk = 0;
 106    for(difference_type i = a.size() - 1; i >= 1; --i)
 107    {
 108       yk2 = yk1;
 109       yk1 = yk;
 110       yk = 2 * x * yk1 - yk2 + a[i];
 111    }
 112    return a[0] / 2 + yk * x - yk1;
 113 }
 114
 115
 116 template <typename T>
 117 class polynomial;
 118
 119 namespace detail {
 120
 121 /**
 122 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
 123 * Chapter 4.6.1, Algorithm D: Division of polynomials over a field.
 124 *
 125 * @tparam  T   Coefficient type, must be not be an integer.
 126 *
 127 * Template-parameter T actually must be a field but we don't currently have that
 128 * subtlety of distinction.
 129 */
 130 template <typename T, typename N>
 131 BOOST_DEDUCED_TYPENAME disable_if_c<std::numeric_limits<T>::is_integer, void >::type
 132 division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
 133 {
 134     q[k] = u[n + k] / v[n];
 135     for (N j = n + k; j > k;)
 136     {
 137         j--;
 138         u[j] -= q[k] * v[j - k];
 139     }
 140 }
 141
 142 template <class T, class N>
 143 T integer_power(T t, N n)
 144 {
 145    switch(n)
 146    {
 147    case 0:
 148       return static_cast<T>(1u);
 149    case 1:
 150       return t;
 151    case 2:
 152       return t * t;
 153    case 3:
 154       return t * t * t;
 155    }
 156    T result = integer_power(t, n / 2);
 157    result *= result;
 158    if(n & 1)
 159       result *= t;
 160    return result;
 161 }
 162
 163
 164 /**
 165 * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
 166 * Chapter 4.6.1, Algorithm R: Pseudo-division of polynomials.
 167 *
 168 * @tparam  T   Coefficient type, must be an integer.
 169 *
 170 * Template-parameter T actually must be a unique factorization domain but we
 171 * don't currently have that subtlety of distinction.
 172 */
 173 template <typename T, typename N>
 174 BOOST_DEDUCED_TYPENAME enable_if_c<std::numeric_limits<T>::is_integer, void >::type
 175 division_impl(polynomial<T> &q, polynomial<T> &u, const polynomial<T>& v, N n, N k)
 176 {
 177     q[k] = u[n + k] * integer_power(v[n], k);
 178     for (N j = n + k; j > 0;)
 179     {
 180         j--;
 181         u[j] = v[n] * u[j] - (j < k ? T(0) : u[n + k] * v[j - k]);
 182     }
 183 }
 184
 185
 186 /**
 187  * Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998
 188  * Chapter 4.6.1, Algorithm D and R: Main loop.
 189  *
 190  * @param   u   Dividend.
 191  * @param   v   Divisor.
 192  */
 193 template <typename T>
 194 std::pair< polynomial<T>, polynomial<T> >
 195 division(polynomial<T> u, const polynomial<T>& v)
 196 {
 197     BOOST_ASSERT(v.size() <= u.size());
 198     BOOST_ASSERT(v);
 199     BOOST_ASSERT(u);
 200
 201     typedef typename polynomial<T>::size_type N;
 202
 203     N const m = u.size() - 1, n = v.size() - 1;
 204     N k = m - n;
 205     polynomial<T> q;
 206     q.data().resize(m - n + 1);
 207
 208     do
 209     {
 210         division_impl(q, u, v, n, k);
 211     }
 212     while (k-- != 0);
 213     u.data().resize(n);
 214     u.normalize(); // Occasionally, the remainder is zeroes.
 215     return std::make_pair(q, u);
 216 }
 217
 218 //
 219 // These structures are the same as the void specializations of the functors of the same name
 220 // in the std lib from C++14 onwards:
 221 //
 222 struct negate
 223 {
 224    template <class T>
 225    T operator()(T const &x) const
 226    {
 227       return -x;
 228    }
 229 };
 230
 231 struct plus
 232 {
 233    template <class T, class U>
 234    T operator()(T const &x, U const& y) const
 235    {
 236       return x + y;
 237    }
 238 };
 239
 240 struct minus
 241 {
 242    template <class T, class U>
 243    T operator()(T const &x, U const& y) const
 244    {
 245       return x - y;
 246    }
 247 };
 248
 249 } // namespace detail
 250
 251 /**
 252  * Returns the zero element for multiplication of polynomials.
 253  */
 254 template <class T>
 255 polynomial<T> zero_element(std::multiplies< polynomial<T> >)
 256 {
 257     return polynomial<T>();
 258 }
 259
 260 template <class T>
 261 polynomial<T> identity_element(std::multiplies< polynomial<T> >)
 262 {
 263     return polynomial<T>(T(1));
 264 }
 265
 266 /* Calculates a / b and a % b, returning the pair (quotient, remainder) together
 267  * because the same amount of computation yields both.
 268  * This function is not defined for division by zero: user beware.
 269  */
 270 template <typename T>
 271 std::pair< polynomial<T>, polynomial<T> >
 272 quotient_remainder(const polynomial<T>& dividend, const polynomial<T>& divisor)
 273 {
 274     BOOST_ASSERT(divisor);
 275     if (dividend.size() < divisor.size())
 276         return std::make_pair(polynomial<T>(), dividend);
 277     return detail::division(dividend, divisor);
 278 }
 279
 280
 281 template <class T>
 282 class polynomial
 283 {
 284 public:
 285    // typedefs:
 286    typedef typename std::vector<T>::value_type value_type;
 287    typedef typename std::vector<T>::size_type size_type;
 288
 289    // construct:
 290    polynomial(){}
 291
 292    template <class U>
 293    polynomial(const U* data, unsigned order)
 294       : m_data(data, data + order + 1)
 295    {
 296        normalize();
 297    }
 298
 299    template <class I>
 300    polynomial(I first, I last)
 301    : m_data(first, last)
 302    {
 303        normalize();
 304    }
 305
 306 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 307    polynomial(std::vector<T>&& p) : m_data(std::move(p))
 308    {
 309       normalize();
 310    }
 311 #endif
 312
 313    template <class U>
 314    explicit polynomial(const U& point, typename boost::enable_if<boost::is_convertible<U, T> >::type* = 0)
 315    {
 316        if (point != U(0))
 317           m_data.push_back(point);
 318    }
 319
 320 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 321    // move:
 322    polynomial(polynomial&& p) BOOST_NOEXCEPT
 323       : m_data(std::move(p.m_data)) { }
 324 #endif
 325
 326    // copy:
 327    polynomial(const polynomial& p)
 328       : m_data(p.m_data) { }
 329
 330    template <class U>
 331    polynomial(const polynomial<U>& p)
 332    {
 333       m_data.resize(p.size());
 334       for(unsigned i = 0; i < p.size(); ++i)
 335       {
 336          m_data[i] = boost::math::tools::real_cast<T>(p[i]);
 337       }
 338    }
 339 #ifdef BOOST_MATH_HAS_IS_CONST_ITERABLE
 340     template <class Range>
 341     explicit polynomial(const Range& r, typename boost::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = 0)
 342        : polynomial(r.begin(), r.end())
 343     {
 344     }
 345 #endif
 346 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !BOOST_WORKAROUND(BOOST_GCC_VERSION, < 40500)
 347     polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l))
 348     {
 349     }
 350
 351     polynomial&
 352     operator=(std::initializer_list<T> l)
 353     {
 354         m_data.assign(std::begin(l), std::end(l));
 355         normalize();
 356         return *this;
 357     }
 358 #endif
 359
 360
 361    // access:
 362    size_type size() const { return m_data.size(); }
 363    size_type degree() const
 364    {
 365        if (size() == 0)
 366            throw std::logic_error("degree() is undefined for the zero polynomial.");
 367        return m_data.size() - 1;
 368    }
 369    value_type& operator[](size_type i)
 370    {
 371       return m_data[i];
 372    }
 373    const value_type& operator[](size_type i) const
 374    {
 375       return m_data[i];
 376    }
 377
 378    T evaluate(T z) const
 379    {
 380       return this->operator()(z);
 381    }
 382
 383    T operator()(T z) const
 384    {
 385       return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial(&m_data[0], z, m_data.size()) : T(0);
 386    }
 387    std::vector<T> chebyshev() const
 388    {
 389       return polynomial_to_chebyshev(m_data);
 390    }
 391
 392    std::vector<T> const& data() const
 393    {
 394        return m_data;
 395    }
 396
 397    std::vector<T> & data()
 398    {
 399        return m_data;
 400    }
 401
 402 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 403    polynomial<T> prime() const
 404    {
 405 #ifdef BOOST_MSVC
 406       // Disable int->float conversion warning:
 407 #pragma warning(push)
 408 #pragma warning(disable:4244)
 409 #endif
 410       if (m_data.size() == 0)
 411       {
 412         return polynomial<T>({});
 413       }
 414
 415       std::vector<T> p_data(m_data.size() - 1);
 416       for (size_t i = 0; i < p_data.size(); ++i) {
 417           p_data[i] = m_data[i+1]*static_cast<T>(i+1);
 418       }
 419       return polynomial<T>(std::move(p_data));
 420 #ifdef BOOST_MSVC
 421 #pragma warning(pop)
 422 #endif
 423    }
 424
 425    polynomial<T> integrate() const
 426    {
 427       std::vector<T> i_data(m_data.size() + 1);
 428       // Choose integration constant such that P(0) = 0.
 429       i_data[0] = T(0);
 430       for (size_t i = 1; i < i_data.size(); ++i)
 431       {
 432           i_data[i] = m_data[i-1]/static_cast<T>(i);
 433       }
 434       return polynomial<T>(std::move(i_data));
 435    }
 436
 437    // operators:
 438    polynomial& operator =(polynomial&& p) BOOST_NOEXCEPT
 439    {
 440        m_data = std::move(p.m_data);
 441        return *this;
 442    }
 443 #endif
 444    polynomial& operator =(const polynomial& p)
 445    {
 446        m_data = p.m_data;
 447        return *this;
 448    }
 449
 450    template <class U>
 451    typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator +=(const U& value)
 452    {
 453        addition(value);
 454        normalize();
 455        return *this;
 456    }
 457
 458    template <class U>
 459    typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator -=(const U& value)
 460    {
 461        subtraction(value);
 462        normalize();
 463        return *this;
 464    }
 465
 466    template <class U>
 467    typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator *=(const U& value)
 468    {
 469       multiplication(value);
 470       normalize();
 471       return *this;
 472    }
 473
 474    template <class U>
 475    typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator /=(const U& value)
 476    {
 477        division(value);
 478        normalize();
 479        return *this;
 480    }
 481
 482    template <class U>
 483    typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial&>::type operator %=(const U& /*value*/)
 484    {
 485        // We can always divide by a scalar, so there is no remainder:
 486        this->set_zero();
 487        return *this;
 488    }
 489
 490    template <class U>
 491    polynomial& operator +=(const polynomial<U>& value)
 492    {
 493       addition(value);
 494       normalize();
 495       return *this;
 496    }
 497
 498    template <class U>
 499    polynomial& operator -=(const polynomial<U>& value)
 500    {
 501        subtraction(value);
 502        normalize();
 503        return *this;
 504    }
 505
 506    template <typename U, typename V>
 507    void multiply(const polynomial<U>& a, const polynomial<V>& b) {
 508        if (!a || !b)
 509        {
 510            this->set_zero();
 511            return;
 512        }
 513        std::vector<T> prod(a.size() + b.size() - 1, T(0));
 514        for (unsigned i = 0; i < a.size(); ++i)
 515            for (unsigned j = 0; j < b.size(); ++j)
 516                prod[i+j] += a.m_data[i] * b.m_data[j];
 517        m_data.swap(prod);
 518    }
 519
 520    template <class U>
 521    polynomial& operator *=(const polynomial<U>& value)
 522    {
 523       this->multiply(*this, value);
 524       return *this;
 525    }
 526
 527    template <typename U>
 528    polynomial& operator /=(const polynomial<U>& value)
 529    {
 530        *this = quotient_remainder(*this, value).first;
 531        return *this;
 532    }
 533
 534    template <typename U>
 535    polynomial& operator %=(const polynomial<U>& value)
 536    {
 537        *this = quotient_remainder(*this, value).second;
 538        return *this;
 539    }
 540
 541    template <typename U>
 542    polynomial& operator >>=(U const &n)
 543    {
 544        BOOST_ASSERT(n <= m_data.size());
 545        m_data.erase(m_data.begin(), m_data.begin() + n);
 546        return *this;
 547    }
 548
 549    template <typename U>
 550    polynomial& operator <<=(U const &n)
 551    {
 552        m_data.insert(m_data.begin(), n, static_cast<T>(0));
 553        normalize();
 554        return *this;
 555    }
 556
 557    // Convenient and efficient query for zero.
 558    bool is_zero() const
 559    {
 560        return m_data.empty();
 561    }
 562
 563    // Conversion to bool.
 564 #ifdef BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS
 565    typedef bool (polynomial::*unmentionable_type)() const;
 566
 567    BOOST_FORCEINLINE operator unmentionable_type() const
 568    {
 569        return is_zero() ? false : &polynomial::is_zero;
 570    }
 571 #else
 572    BOOST_FORCEINLINE explicit operator bool() const
 573    {
 574        return !m_data.empty();
 575    }
 576 #endif
 577
 578    // Fast way to set a polynomial to zero.
 579    void set_zero()
 580    {
 581        m_data.clear();
 582    }
 583
 584     /** Remove zero coefficients 'from the top', that is for which there are no
 585     *        non-zero coefficients of higher degree. */
 586    void normalize()
 587    {
 588 #ifndef BOOST_NO_CXX11_LAMBDAS
 589       m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), [](const T& x)->bool { return x != T(0); }).base(), m_data.end());
 590 #else
 591        using namespace boost::lambda;
 592        m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end());
 593 #endif
 594    }
 595
 596 private:
 597     template <class U, class R>
 598     polynomial& addition(const U& value, R op)
 599     {
 600         if(m_data.size() == 0)
 601             m_data.resize(1, 0);
 602         m_data[0] = op(m_data[0], value);
 603         return *this;
 604     }
 605
 606     template <class U>
 607     polynomial& addition(const U& value)
 608     {
 609         return addition(value, detail::plus());
 610     }
 611
 612     template <class U>
 613     polynomial& subtraction(const U& value)
 614     {
 615         return addition(value, detail::minus());
 616     }
 617
 618     template <class U, class R>
 619     polynomial& addition(const polynomial<U>& value, R op)
 620     {
 621         if (m_data.size() < value.size())
 622             m_data.resize(value.size(), 0);
 623         for(size_type i = 0; i < value.size(); ++i)
 624             m_data[i] = op(m_data[i], value[i]);
 625         return *this;
 626     }
 627
 628     template <class U>
 629     polynomial& addition(const polynomial<U>& value)
 630     {
 631         return addition(value, detail::plus());
 632     }
 633
 634     template <class U>
 635     polynomial& subtraction(const polynomial<U>& value)
 636     {
 637         return addition(value, detail::minus());
 638     }
 639
 640     template <class U>
 641     polynomial& multiplication(const U& value)
 642     {
 643 #ifndef BOOST_NO_CXX11_LAMBDAS
 644        std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x * value; });
 645 #else
 646         using namespace boost::lambda;
 647         std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 * value));
 648 #endif
 649         return *this;
 650     }
 651
 652     template <class U>
 653     polynomial& division(const U& value)
 654     {
 655 #ifndef BOOST_NO_CXX11_LAMBDAS
 656        std::transform(m_data.begin(), m_data.end(), m_data.begin(), [&](const T& x)->T { return x / value; });
 657 #else
 658         using namespace boost::lambda;
 659         std::transform(m_data.begin(), m_data.end(), m_data.begin(), ret<T>(_1 / value));
 660 #endif
 661         return *this;
 662     }
 663
 664     std::vector<T> m_data;
 665 };
 666
 667
 668 template <class T>
 669 inline polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b)
 670 {
 671    polynomial<T> result(a);
 672    result += b;
 673    return result;
 674 }
 675 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 676 template <class T>
 677 inline polynomial<T> operator + (polynomial<T>&& a, const polynomial<T>& b)
 678 {
 679    a += b;
 680    return a;
 681 }
 682 template <class T>
 683 inline polynomial<T> operator + (const polynomial<T>& a, polynomial<T>&& b)
 684 {
 685    b += a;
 686    return b;
 687 }
 688 template <class T>
 689 inline polynomial<T> operator + (polynomial<T>&& a, polynomial<T>&& b)
 690 {
 691    a += b;
 692    return a;
 693 }
 694 #endif
 695
 696 template <class T>
 697 inline polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b)
 698 {
 699    polynomial<T> result(a);
 700    result -= b;
 701    return result;
 702 }
 703 #ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
 704 template <class T>
 705 inline polynomial<T> operator - (polynomial<T>&& a, const polynomial<T>& b)
 706 {
 707    a -= b;
 708    return a;
 709 }
 710 template <class T>
 711 inline polynomial<T> operator - (const polynomial<T>& a, polynomial<T>&& b)
 712 {
 713    b -= a;
 714    return -b;
 715 }
 716 template <class T>
 717 inline polynomial<T> operator - (polynomial<T>&& a, polynomial<T>&& b)
 718 {
 719    a -= b;
 720    return a;
 721 }
 722 #endif
 723
 724 template <class T>
 725 inline polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b)
 726 {
 727    polynomial<T> result;
 728    result.multiply(a, b);
 729    return result;
 730 }
 731
 732 template <class T>
 733 inline polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b)
 734 {
 735    return quotient_remainder(a, b).first;
 736 }
 737
 738 template <class T>
 739 inline polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b)
 740 {
 741    return quotient_remainder(a, b).second;
 742 }
 743
 744 template <class T, class U>
 745 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (polynomial<T> a, const U& b)
 746 {
 747    a += b;
 748    return a;
 749 }
 750
 751 template <class T, class U>
 752 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (polynomial<T> a, const U& b)
 753 {
 754    a -= b;
 755    return a;
 756 }
 757
 758 template <class T, class U>
 759 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (polynomial<T> a, const U& b)
 760 {
 761    a *= b;
 762    return a;
 763 }
 764
 765 template <class T, class U>
 766 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator / (polynomial<T> a, const U& b)
 767 {
 768    a /= b;
 769    return a;
 770 }
 771
 772 template <class T, class U>
 773 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator % (const polynomial<T>&, const U&)
 774 {
 775    // Since we can always divide by a scalar, result is always an empty polynomial:
 776    return polynomial<T>();
 777 }
 778
 779 template <class U, class T>
 780 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator + (const U& a, polynomial<T> b)
 781 {
 782    b += a;
 783    return b;
 784 }
 785
 786 template <class U, class T>
 787 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator - (const U& a, polynomial<T> b)
 788 {
 789    b -= a;
 790    return -b;
 791 }
 792
 793 template <class U, class T>
 794 inline typename boost::enable_if_c<boost::is_constructible<T, U>::value, polynomial<T> >::type operator * (const U& a, polynomial<T> b)
 795 {
 796    b *= a;
 797    return b;
 798 }
 799
 800 template <class T>
 801 bool operator == (const polynomial<T> &a, const polynomial<T> &b)
 802 {
 803     return a.data() == b.data();
 804 }
 805
 806 template <class T>
 807 bool operator != (const polynomial<T> &a, const polynomial<T> &b)
 808 {
 809     return a.data() != b.data();
 810 }
 811
 812 template <typename T, typename U>
 813 polynomial<T> operator >> (polynomial<T> a, const U& b)
 814 {
 815     a >>= b;
 816     return a;
 817 }
 818
 819 template <typename T, typename U>
 820 polynomial<T> operator << (polynomial<T> a, const U& b)
 821 {
 822     a <<= b;
 823     return a;
 824 }
 825
 826 // Unary minus (negate).
 827 template <class T>
 828 polynomial<T> operator - (polynomial<T> a)
 829 {
 830     std::transform(a.data().begin(), a.data().end(), a.data().begin(), detail::negate());
 831     return a;
 832 }
 833
 834 template <class T>
 835 bool odd(polynomial<T> const &a)
 836 {
 837     return a.size() > 0 && a[0] != static_cast<T>(0);
 838 }
 839
 840 template <class T>
 841 bool even(polynomial<T> const &a)
 842 {
 843     return !odd(a);
 844 }
 845
 846 template <class T>
 847 polynomial<T> pow(polynomial<T> base, int exp)
 848 {
 849     if (exp < 0)
 850         return policies::raise_domain_error(
 851                 "boost::math::tools::pow<%1%>",
 852                 "Negative powers are not supported for polynomials.",
 853                 base, policies::policy<>());
 854         // if the policy is ignore_error or errno_on_error, raise_domain_error
 855         // will return std::numeric_limits<polynomial<T>>::quiet_NaN(), which
 856         // defaults to polynomial<T>(), which is the zero polynomial
 857     polynomial<T> result(T(1));
 858     if (exp & 1)
 859         result = base;
 860     /* "Exponentiation by squaring" */
 861     while (exp >>= 1)
 862     {
 863         base *= base;
 864         if (exp & 1)
 865             result *= base;
 866     }
 867     return result;
 868 }
 869
 870 template <class charT, class traits, class T>
 871 inline std::basic_ostream<charT, traits>& operator << (std::basic_ostream<charT, traits>& os, const polynomial<T>& poly)
 872 {
 873    os << "{ ";
 874    for(unsigned i = 0; i < poly.size(); ++i)
 875    {
 876       if(i) os << ", ";
 877       os << poly[i];
 878    }
 879    os << " }";
 880    return os;
 881 }
 882
 883 } // namespace tools
 884 } // namespace math
 885 } // namespace boost
 886
 887 //
 888 // Polynomial specific overload of gcd algorithm:
 889 //
 890 #include <boost/math/tools/polynomial_gcd.hpp>
 891
 892 #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP