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