863 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Jeremy William Murphy 2015.3 4 5// Use, modification and distribution are subject to the6// Boost Software License, Version 1.0. (See accompanying file7// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)8 9#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_HPP10#define BOOST_MATH_TOOLS_POLYNOMIAL_HPP11 12#ifdef _MSC_VER13#pragma once14#endif15 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 32namespace boost{ namespace math{ namespace tools{33 34template <class T>35BOOST_MATH_GPU_ENABLED T chebyshev_coefficient(unsigned n, unsigned m)36{37 BOOST_MATH_STD_USING38 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 58template <class Seq>59BOOST_MATH_GPU_ENABLED 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 94template <class Seq, class T>95BOOST_MATH_GPU_ENABLED 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 112template <typename T>113class polynomial;114 115namespace detail {116 117/**118* Knuth, The Art of Computer Programming: Volume 2, Third edition, 1998119* 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 that124* subtlety of distinction.125*/126template <typename T, typename N>127BOOST_MATH_GPU_ENABLED typename std::enable_if<!std::numeric_limits<T>::is_integer, void >::type128division_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 138template <class T, class N>139BOOST_MATH_GPU_ENABLED 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, 1998162* 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 we167* don't currently have that subtlety of distinction.168*/169template <typename T, typename N>170BOOST_MATH_GPU_ENABLED typename std::enable_if<std::numeric_limits<T>::is_integer, void >::type171division_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, 1998184 * Chapter 4.6.1, Algorithm D and R: Main loop.185 *186 * @param u Dividend.187 * @param v Divisor.188 */189template <typename T>190BOOST_MATH_GPU_ENABLED std::pair< polynomial<T>, polynomial<T> >191division(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 do205 {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 name216// in the std lib from C++14 onwards:217//218struct negate219{220 template <class T>221 BOOST_MATH_GPU_ENABLED T operator()(T const &x) const222 {223 return -x;224 }225};226 227struct plus228{229 template <class T, class U>230 BOOST_MATH_GPU_ENABLED T operator()(T const &x, U const& y) const231 {232 return x + y;233 }234};235 236struct minus237{238 template <class T, class U>239 BOOST_MATH_GPU_ENABLED T operator()(T const &x, U const& y) const240 {241 return x - y;242 }243};244 245} // namespace detail246 247/**248 * Returns the zero element for multiplication of polynomials.249 */250template <class T>251BOOST_MATH_GPU_ENABLED polynomial<T> zero_element(std::multiplies< polynomial<T> >)252{253 return polynomial<T>();254}255 256template <class T>257BOOST_MATH_GPU_ENABLED 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) together263 * because the same amount of computation yields both.264 * This function is not defined for division by zero: user beware.265 */266template <typename T>267BOOST_MATH_GPU_ENABLED std::pair< polynomial<T>, polynomial<T> >268quotient_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 277template <class T>278class polynomial279{280public: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 BOOST_MATH_GPU_ENABLED polynomial()= default;287 288 template <class U>289 BOOST_MATH_GPU_ENABLED polynomial(const U* data, unsigned order)290 : m_data(data, data + order + 1)291 {292 normalize();293 }294 295 template <class Iterator>296 BOOST_MATH_GPU_ENABLED polynomial(Iterator first, Iterator last)297 : m_data(first, last)298 {299 normalize();300 }301 302 template <class Iterator>303 BOOST_MATH_GPU_ENABLED polynomial(Iterator first, unsigned length)304 : m_data(first, std::next(first, length + 1))305 {306 normalize();307 }308 309 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED explicit polynomial(const U& point)316 {317 if (point != U(0))318 m_data.push_back(point);319 }320 321 // move:322 BOOST_MATH_GPU_ENABLED polynomial(polynomial&& p) noexcept323 : m_data(std::move(p.m_data)) { }324 325 // copy:326 BOOST_MATH_GPU_ENABLED polynomial(const polynomial& p)327 : m_data(p.m_data) { }328 329 template <class U>330 BOOST_MATH_GPU_ENABLED 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_ITERABLE339 template <class Range, typename std::enable_if<boost::math::tools::detail::is_const_iterable<Range>::value, bool>::type = true>340 BOOST_MATH_GPU_ENABLED explicit polynomial(const Range& r)341 : polynomial(r.begin(), r.end()) 342 {343 }344#endif345 BOOST_MATH_GPU_ENABLED polynomial(std::initializer_list<T> l) : polynomial(std::begin(l), std::end(l))346 {347 }348 349 polynomial&350 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED size_type size() const { return m_data.size(); }360 BOOST_MATH_GPU_ENABLED size_type degree() const361 {362 if (size() == 0)363 BOOST_MATH_THROW_EXCEPTION(std::logic_error("degree() is undefined for the zero polynomial."));364 return m_data.size() - 1;365 }366 BOOST_MATH_GPU_ENABLED value_type& operator[](size_type i)367 {368 return m_data[i];369 }370 BOOST_MATH_GPU_ENABLED const value_type& operator[](size_type i) const371 {372 return m_data[i];373 }374 375 BOOST_MATH_GPU_ENABLED T evaluate(T z) const376 {377 return this->operator()(z);378 }379 380 BOOST_MATH_GPU_ENABLED T operator()(T z) const381 {382 return m_data.size() > 0 ? boost::math::tools::evaluate_polynomial((m_data).data(), z, m_data.size()) : T(0);383 }384 BOOST_MATH_GPU_ENABLED std::vector<T> chebyshev() const385 {386 return polynomial_to_chebyshev(m_data);387 }388 389 BOOST_MATH_GPU_ENABLED std::vector<T> const& data() const390 {391 return m_data;392 }393 394 BOOST_MATH_GPU_ENABLED std::vector<T> & data()395 {396 return m_data;397 }398 399 BOOST_MATH_GPU_ENABLED polynomial<T> prime() const400 {401#ifdef _MSC_VER402 // Disable int->float conversion warning:403#pragma warning(push)404#pragma warning(disable:4244)405#endif406 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_VER417#pragma warning(pop)418#endif419 }420 421 BOOST_MATH_GPU_ENABLED polynomial<T> integrate() const422 {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 BOOST_MATH_GPU_ENABLED polynomial& operator =(polynomial&& p) noexcept435 {436 m_data = std::move(p.m_data);437 return *this;438 }439 440 BOOST_MATH_GPU_ENABLED polynomial& operator =(const polynomial& p)441 {442 m_data = p.m_data;443 return *this;444 }445 446 template <class U>447 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED polynomial& operator +=(const polynomial<U>& value)488 {489 addition(value);490 normalize();491 return *this;492 }493 494 template <class U>495 BOOST_MATH_GPU_ENABLED polynomial& operator -=(const polynomial<U>& value)496 {497 subtraction(value);498 normalize();499 return *this;500 }501 502 template <typename U, typename V>503 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED polynomial& operator *=(const polynomial<U>& value)518 {519 this->multiply(*this, value);520 return *this;521 }522 523 template <typename U>524 BOOST_MATH_GPU_ENABLED polynomial& operator /=(const polynomial<U>& value)525 {526 *this = quotient_remainder(*this, value).first;527 return *this;528 }529 530 template <typename U>531 BOOST_MATH_GPU_ENABLED polynomial& operator %=(const polynomial<U>& value)532 {533 *this = quotient_remainder(*this, value).second;534 return *this;535 }536 537 template <typename U>538 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED bool is_zero() const555 {556 return m_data.empty();557 }558 559 // Conversion to bool.560 BOOST_MATH_GPU_ENABLED inline explicit operator bool() const561 {562 return !m_data.empty();563 }564 565 // Fast way to set a polynomial to zero.566 BOOST_MATH_GPU_ENABLED void set_zero()567 {568 m_data.clear();569 }570 571 /** Remove zero coefficients 'from the top', that is for which there are no572 * non-zero coefficients of higher degree. */573 BOOST_MATH_GPU_ENABLED 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 578private:579 template <class U, class R>580 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED polynomial& addition(const U& value)590 {591 return addition(value, detail::plus());592 }593 594 template <class U>595 BOOST_MATH_GPU_ENABLED polynomial& subtraction(const U& value)596 {597 return addition(value, detail::minus());598 }599 600 template <class U, class R>601 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED polynomial& addition(const polynomial<U>& value)612 {613 return addition(value, detail::plus());614 }615 616 template <class U>617 BOOST_MATH_GPU_ENABLED polynomial& subtraction(const polynomial<U>& value)618 {619 return addition(value, detail::minus());620 }621 622 template <class U>623 BOOST_MATH_GPU_ENABLED 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 BOOST_MATH_GPU_ENABLED 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 640template <class T>641BOOST_MATH_GPU_ENABLED 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 648template <class T>649BOOST_MATH_GPU_ENABLED inline polynomial<T> operator + (polynomial<T>&& a, const polynomial<T>& b)650{651 a += b;652 return std::move(a);653}654template <class T>655BOOST_MATH_GPU_ENABLED inline polynomial<T> operator + (const polynomial<T>& a, polynomial<T>&& b)656{657 b += a;658 return b;659}660template <class T>661BOOST_MATH_GPU_ENABLED inline polynomial<T> operator + (polynomial<T>&& a, polynomial<T>&& b)662{663 a += b;664 return a;665}666 667template <class T>668BOOST_MATH_GPU_ENABLED 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 675template <class T>676BOOST_MATH_GPU_ENABLED inline polynomial<T> operator - (polynomial<T>&& a, const polynomial<T>& b)677{678 a -= b;679 return a;680}681template <class T>682BOOST_MATH_GPU_ENABLED inline polynomial<T> operator - (const polynomial<T>& a, polynomial<T>&& b)683{684 b -= a;685 return -b;686}687template <class T>688BOOST_MATH_GPU_ENABLED inline polynomial<T> operator - (polynomial<T>&& a, polynomial<T>&& b)689{690 a -= b;691 return a;692}693 694template <class T>695BOOST_MATH_GPU_ENABLED 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 702template <class T>703BOOST_MATH_GPU_ENABLED inline polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b)704{705 return quotient_remainder(a, b).first;706}707 708template <class T>709BOOST_MATH_GPU_ENABLED inline polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b)710{711 return quotient_remainder(a, b).second;712}713 714template <class T, class U>715BOOST_MATH_GPU_ENABLED 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 721template <class T, class U>722BOOST_MATH_GPU_ENABLED 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 728template <class T, class U>729BOOST_MATH_GPU_ENABLED 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 735template <class T, class U>736BOOST_MATH_GPU_ENABLED 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 742template <class T, class U>743BOOST_MATH_GPU_ENABLED 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 749template <class U, class T>750BOOST_MATH_GPU_ENABLED 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 756template <class U, class T>757BOOST_MATH_GPU_ENABLED 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 763template <class U, class T>764BOOST_MATH_GPU_ENABLED 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 770template <class T>771BOOST_MATH_GPU_ENABLED bool operator == (const polynomial<T> &a, const polynomial<T> &b)772{773 return a.data() == b.data();774}775 776template <class T>777BOOST_MATH_GPU_ENABLED bool operator != (const polynomial<T> &a, const polynomial<T> &b)778{779 return a.data() != b.data();780}781 782template <typename T, typename U>783BOOST_MATH_GPU_ENABLED polynomial<T> operator >> (polynomial<T> a, const U& b)784{785 a >>= b;786 return a;787}788 789template <typename T, typename U>790BOOST_MATH_GPU_ENABLED polynomial<T> operator << (polynomial<T> a, const U& b)791{792 a <<= b;793 return a;794}795 796// Unary minus (negate).797template <class T>798BOOST_MATH_GPU_ENABLED 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 804template <class T>805BOOST_MATH_GPU_ENABLED bool odd(polynomial<T> const &a)806{807 return a.size() > 0 && a[0] != static_cast<T>(0);808}809 810template <class T>811BOOST_MATH_GPU_ENABLED bool even(polynomial<T> const &a)812{813 return !odd(a);814}815 816template <class T>817BOOST_MATH_GPU_ENABLED 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_error825 // will return std::numeric_limits<polynomial<T>>::quiet_NaN(), which826 // defaults to polynomial<T>(), which is the zero polynomial827 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 840template <class charT, class traits, class T>841BOOST_MATH_GPU_ENABLED 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 tools854} // namespace math855} // namespace boost856 857//858// Polynomial specific overload of gcd algorithm:859//860#include <boost/math/tools/polynomial_gcd.hpp>861 862#endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP863