-
Notifications
You must be signed in to change notification settings - Fork 264
Expand file tree
/
Copy pathpolynomial.qbk
More file actions
205 lines (161 loc) · 7.18 KB
/
Copy pathpolynomial.qbk
File metadata and controls
205 lines (161 loc) · 7.18 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
[section:polynomials Polynomials]
[h4 Synopsis]
``
#include <boost/math/tools/polynomial.hpp>
``
namespace boost { namespace math {
namespace tools {
template <class T>
class polynomial
{
public:
// typedefs:
typedef typename std::vector<T>::value_type value_type;
typedef typename std::vector<T>::size_type size_type;
// construct:
polynomial(){}
template <class U>
polynomial(const U* data, unsigned order);
template <class Iterator>
polynomial(Iterator first, Iterator last);
template <class U>
explicit polynomial(const U& point,
typename std::enable_if<std::is_convertible<U, T> >::type* = nullptr);
template <class Range>
explicit polynomial(const Range& r,
typename std::enable_if<boost::math::tools::detail::is_const_iterable<Range> >::type* = nullptr); // C++14
polynomial(std::initializer_list<T> l); // C++11
polynomial(std::vector<T>&& p);
// access:
size_type size()const;
size_type degree()const;
value_type& operator[](size_type i);
const value_type& operator[](size_type i)const;
std::vector<T> const& data() const;
std::vector<T>& data();
explicit operator bool() const;
polynomial<T> prime() const;
polynomial<T> integrate() const;
T operator()(T z) const;
// Only available if Eigen library is available:
std::vector<std::complex<T>> roots() const;
// modify:
void set_zero();
// operators:
template <class U>
polynomial& operator +=(const U& value);
template <class U>
polynomial& operator -=(const U& value);
template <class U>
polynomial& operator *=(const U& value);
template <class U>
polynomial& operator /=(const U& value);
template <class U>
polynomial& operator %=(const U& value);
template <class U>
polynomial& operator +=(const polynomial<U>& value);
template <class U>
polynomial& operator -=(const polynomial<U>& value);
template <class U>
polynomial& operator *=(const polynomial<U>& value);
template <class U>
polynomial& operator /=(const polynomial<U>& value);
template <class U>
polynomial& operator %=(const polynomial<U>& value);
};
template <class T>
polynomial<T> operator + (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator - (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator * (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator / (const polynomial<T>& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator % (const polynomial<T>& a, const polynomial<T>& b);
template <class T, class U>
polynomial<T> operator + (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator - (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator * (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator / (const polynomial<T>& a, const U& b);
template <class T, class U>
polynomial<T> operator % (const polynomial<T>& a, const U& b);
template <class U, class T>
polynomial<T> operator + (const U& a, const polynomial<T>& b);
template <class U, class T>
polynomial<T> operator - (const U& a, const polynomial<T>& b);
template <class U, class T>
polynomial<T> operator * (const U& a, const polynomial<T>& b);
template <class T>
polynomial<T> operator - (const polynomial<T>& a);
template <class T>
polynomial<T> operator >>= (const U& a);
template <class T>
polynomial<T> operator >> (polynomial<T> const &a, const U& b);
template <class T>
polynomial<T> operator <<= (const U& a);
template <class T>
polynomial<T> operator << (polynomial<T> const &a, const U& b);
template <class T>
bool operator == (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
bool operator != (const polynomial<T> &a, const polynomial<T> &b);
template <class T>
polynomial<T> pow(polynomial<T> base, int exp);
template <class charT, class traits, class T>
std::basic_ostream<charT, traits>& operator <<
(std::basic_ostream<charT, traits>& os, const polynomial<T>& poly);
template <typename T>
std::pair< polynomial<T>, polynomial<T> >
quotient_remainder(const polynomial<T>& a, const polynomial<T>& b);
} // namespace tools
}} // namespace boost { namespace math
[h4 Description]
This is a somewhat trivial class for polynomial manipulation.
See
* [@https://en.wikipedia.org/wiki/Polynomial Polynomial] and
* Donald E. Knuth, The Art of Computer Programming: Volume 2, Third edition, (1998)
Chapter 4.6.1, Algorithm D: Division of polynomials over a field.
Implementation is currently of the "naive" variety, with [bigo](N[super 2])
multiplication, for example. This class should not be used in
high-performance computing environments: it is intended for the
simple manipulation of small polynomials, typically generated
for special function approximation.
It does has division for polynomials over a [@https://en.wikipedia.org/wiki/Field_%28mathematics%29 field]
(here floating point, complex, etc)
and over a unique factorization domain (integers).
Division of polynomials over a field is compatible with
[@https://en.wikipedia.org/wiki/Euclidean_algorithm Euclidean GCD].
Division of polynomials over a UFD is compatible with the subresultant algorithm for GCD (implemented as subresultant_gcd), but a serious word of warning is required: the intermediate value swell of that algorithm will cause single-precision integral types to overflow very easily. So although the algorithm will work on single-precision integral types, an overload of the gcd function is only provided for polynomials with multi-precision integral types, to prevent nasty surprises. This is done somewhat crudely by disabling the overload for non-POD integral types.
Advanced manipulations: the FFT, factorisation etc are
not currently provided. Submissions for these are of course welcome :-)
[h4:polynomial_examples Polynomial Arithmetic Examples]
[import ../../example/polynomial_arithmetic.cpp]
[polynomial_arithmetic_0]
[polynomial_arithmetic_1]
[polynomial_arithmetic_2]
for output:
[polynomial_output_1]
[polynomial_arithmetic_3]
for output
[polynomial_output_2]
[h5 Division, Quotient and Remainder]
[polynomial_arithmetic_4]
The source code is at [@../../example/polynomial_arithmetic.cpp polynomial_arithmetic.cpp]
[h4 Roots]
Polynomial roots are recovered by the eigenvalues of the companion matrix.
This requires that the Eigen C++ library to be available; otherwise calling `.roots()` is a compile-time error.
In addition, the polynomial coefficients must be of floating point type, or a static assertion will fail.
[endsect]
[/section:polynomials Polynomials]
[/
Copyright 2006 John Maddock and Paul A. Bristow.
Copyright 2015 Jeremy William Murphy.
Copyright 2024 Nick Thompson.
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or copy at
http://www.boost.org/LICENSE_1_0.txt).
]