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