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