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