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