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