Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
polynomial.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_POLYNOMIAL_H
7#define THEORETICA_POLYNOMIAL_H
8
9#ifndef THEORETICA_NO_PRINT
10#include <sstream>
11#include <ostream>
12#endif
13
14#include "../core/real_analysis.h"
15#include "../algebra/vec.h"
16#include "../complex/complex.h"
17#include "../complex/complex_analysis.h"
18
19
20namespace theoretica {
21
24 template<typename Type = real>
25 class polynomial {
26 public:
27 std::vector<Type> coeff;
28
30 polynomial() : coeff() {}
31
33 polynomial(Type a) : coeff({a}) {}
34
36 polynomial(const std::vector<Type>& c) : coeff(c) {}
37
39 polynomial(std::initializer_list<Type> l) : coeff(l) {}
40
43
44
46 inline Type& at(int i) {
47 return coeff[i];
48 }
49
50
52 inline Type get(int i) const {
53 return coeff[i];
54 }
55
56
58 inline const Type& operator[](unsigned int i) const {
59 return coeff[i];
60 }
61
62
64 inline Type& operator[](unsigned int i) {
65 return coeff[i];
66 }
67
68
70 template<typename EvalType = Type>
71 inline EvalType eval(EvalType x) const {
72
73 if(!coeff.size())
74 return 0;
75
76 EvalType sum = 0;
77
78 // Evaluate using Horner's method
79 for (unsigned int i = 0; i < coeff.size(); ++i)
80 sum = coeff[coeff.size() - i - 1] + x * sum;
81
82 return sum;
83 }
84
85
87 template<typename EvalType = Type>
88 inline EvalType operator()(EvalType x) const {
89 return eval(x);
90 }
91
92
95 inline unsigned int find_order() const {
96
97 for (int i = coeff.size() - 1; i >= 0; --i) {
98 if(abs(coeff[i]) < MACH_EPSILON)
99 return i;
100 }
101
102 return 0;
103 }
104
105
107 inline void trim() {
108
109 for (unsigned int i = coeff.size() - 1; i >= 0; --i) {
110 if(abs(coeff[i]) > MACH_EPSILON)
111 break;
112
113 coeff.pop_back();
114 }
115 }
116
117
119 inline size_t size() const {
120 return coeff.size();
121 }
122
123
125 inline polynomial operator+(const polynomial& p) const {
126
128 r.coeff.resize(max(coeff.size(), p.size()));
129
130 unsigned int i = 0;
131
132 for (; i < min(coeff.size(), p.size()); ++i)
133 r[i] = coeff[i] + p[i];
134
135 if (coeff.size() > p.size())
136 for (; i < coeff.size(); ++i)
137 r[i] = coeff[i];
138 else
139 for (; i < p.size(); ++i)
140 r[i] = p[i];
141
142 return r;
143 }
144
145
147 inline polynomial operator-(const polynomial& p) const {
148
150 r.coeff.resize(max(coeff.size(), p.size()));
151
152 unsigned int i = 0;
153
154 for (; i < min(r.size(), p.size()); ++i)
155 r[i] = coeff[i] - p.get(i);
156
157 if (coeff.size() > p.size())
158 for (; i < coeff.size(); ++i)
159 r[i] = coeff[i];
160 else
161 for (; i < p.size(); ++i)
162 r[i] = -p[i];
163
164 return r;
165 }
166
167
169 inline polynomial operator*(const polynomial& p) const {
170
172 r.coeff.resize(this->size() + p.size() - 1);
173
174 for (unsigned int i = 0; i < size(); ++i) {
175 for (unsigned int j = 0; j < p.size(); ++j) {
176 r[i + j] += coeff[i] * p.get(j);
177 }
178 }
179
180 return r;
181 }
182
183
185 inline polynomial operator/(const polynomial& d) const {
186
187 const unsigned int d_order = d.find_order();
188 const unsigned int this_order = find_order();
189
190 if(d_order == 0 && d.get(0) == 0) {
191 TH_MATH_ERROR("polynomial::operator/", d.get(0), DIV_BY_ZERO);
192 return polynomial(nan());
193 }
194
195 // Remainder
196 polynomial r = *this;
197
198 // Quotient
200 unsigned int i = 0;
201
202 while(i < this_order) {
203
204 // Compute only once the degree of the polynomial
205 const unsigned int r_order = r.find_order();
206
207 // Stop execution if the division is complete
208 // (when the remainder is 0 or has lower degree)
209 if((r_order == 0 && (abs(r.get(0)) < MACH_EPSILON)) || r_order < d_order)
210 break;
211
212 // Simple division between highest degree terms
214 r.get(r_order) / d.get(d_order),
215 r_order - d_order);
216
217 // Add monomial to quotient and subtract the
218 // monomial times the dividend from the remainder
219 q += t;
220 r -= t * d;
221
222 i++;
223 }
224
225 // The algorithm has stopped iterating after a number
226 // of the dividend's degree counts
227 if(i == this_order) {
228 TH_MATH_ERROR("polynomial::operator/", i, NO_ALGO_CONVERGENCE);
229 return polynomial(nan());
230 }
231
232 return q;
233 }
234
235
237 inline polynomial operator*(Type a) const {
238
239 polynomial r = polynomial(*this);
240
241 for (unsigned int i = 0; i < size(); ++i)
242 r.coeff[i] *= a;
243
244 return r;
245 }
246
247
249 inline polynomial operator/(Type a) const {
250
251 if(abs(a) < MACH_EPSILON) {
252 TH_MATH_ERROR("polynomial::operator/", a, DIV_BY_ZERO);
253 }
254
255 polynomial r = polynomial(*this);
256
257 for (unsigned int i = 0; i < size(); ++i)
258 r.coeff[i] /= a;
259
260 return r;
261 }
262
263
266
267 // Make room for the new coefficients
268 if(coeff.size() < p.size())
269 coeff.resize(p.size(), Type(0));
270
271 for (unsigned int i = 0; i < min(size(), p.size()); ++i)
272 coeff[i] += p.get(i);
273
274 return *this;
275 }
276
277
280
281 // Make room for the new coefficients
282 if(coeff.size() < p.size())
283 coeff.resize(p.size(), Type(0));
284
285 for (unsigned int i = 0; i < min(coeff.size(), p.size()); ++i)
286 coeff[i] -= p.get(i);
287
288 return *this;
289 }
290
291
294
296 r.coeff.resize(this->size() + p.size() - 1, Type(0));
297
298 for (unsigned int i = 0; i < coeff.size(); ++i)
299 for (unsigned int j = 0; j < p.size(); ++j)
300 r[i + j] += coeff[i] * p.get(j);
301
302 *this = r;
303 return *this;
304 }
305
306
309
310 for (unsigned int i = 0; i < coeff.size(); ++i)
311 coeff[i] *= a;
312
313 return *this;
314 }
315
316
318 inline polynomial& operator/=(const polynomial& a) {
319 return (*this = (*this / a));
320 }
321
322
325
326 if(abs(a) < MACH_EPSILON) {
327 TH_MATH_ERROR("polynomial::operator/=", a, DIV_BY_ZERO);
328 return *this;
329 }
330
331 for (unsigned int i = 0; i < coeff.size(); ++i)
332 coeff[i] /= a;
333
334 return *this;
335 }
336
337
339 inline bool operator==(const polynomial& other) const {
340
341 const unsigned int n = min(other.size(), this->size());
342 for (unsigned int i = 0; i < n; ++i) {
343 if(abs(other.coeff[i] - coeff[i]) >= MACH_EPSILON)
344 return false;
345 }
346
347 if(size() > other.size()) {
348 for (unsigned int i = 0; i < size(); ++i) {
349 if(abs(coeff[i]) >= MACH_EPSILON)
350 return false;
351 }
352 } else {
353 for (unsigned int i = 0; i < other.size(); ++i) {
354 if(abs(other.coeff[i]) >= MACH_EPSILON)
355 return false;
356 }
357 }
358
359 return true;
360 }
361
362
365 inline auto begin() {
366 return coeff.begin();
367 }
368
369
372 inline auto end() {
373 return coeff.end();
374 }
375
376
379 inline auto cbegin() {
380 return coeff.cbegin();
381 }
382
383
386 inline auto cend() {
387 return coeff.cend();
388 }
389
390
392 inline vec<complex<>, 2> quadratic_roots() const {
393
394 const int order = find_order();
395
396 // Check that the polynomial is quadratic
397 if(order != 2) {
398 TH_MATH_ERROR("quadratic_roots", order, IMPOSSIBLE_OPERATION);
399 return vec<complex<>, 2>({nan(), nan()});
400 }
401
402 const Type p = coeff[1] / coeff[2];
403 const Type q = coeff[0] / coeff[2];
404
405 // Case when 0 is a root
406 if(abs(q) < MACH_EPSILON)
407 return {complex<>(-p), complex<>(0)};
408
409 complex<> z1, z2;
410
411 // Use Vieta's theorem to avoid catastrophic cancellation
412 if(abs(p) > 1) {
413
414 z1 = -sgn(p) * (abs(p) / 2.0
415 + abs(p) * sqrt(complex<>(0.25 - (q / p) / p)));
416 z2 = q / z1;
417
418 } else {
419
420 const complex<> s = sqrt(complex<>(0.25 * square(p) - q));
421 z1 = -p / 2.0 + s;
422 z2 = -p / 2.0 - s;
423 }
424
425 return {z1, z2};
426 }
427
428
431 const std::vector<Type>& roots) {
432
433 polynomial<Type> P = {1};
434
435 // P = product((x - r_i))
436 for (unsigned int i = 0; i < roots.size(); ++i)
437 P *= polynomial<Type>({roots[i] * -1, 1});
438
439 return P;
440 }
441
442
444 inline static polynomial<Type> monomial(Type c, unsigned int order) {
445
447 m.coeff = std::vector<Type>(order + 1, Type(0));
448 m.coeff[order] = c;
449 return m;
450 }
451
452
453 // Friend operators to enable equations of the form
454 // (Type) op. (polynomial<Type>)
455
456 inline friend polynomial<Type> operator+(
457 Type r, const polynomial<Type>& z) {
458
459 return z + polynomial(r);
460 }
461
462 inline friend polynomial<Type> operator-(
463 Type r, const polynomial<Type>& z) {
464
465 return (z * -1) + polynomial(r);
466 }
467
468 inline friend polynomial<Type> operator*(
469 Type r, const polynomial<Type>& z) {
470
471 return z * r;
472 }
473
474
475#ifndef THEORETICA_NO_PRINT
476
478 inline std::string to_string(
479 const std::string& unknown = "x",
480 const std::string& exponentiation = "^") const {
481
482 std::stringstream res;
483 bool flag = false;
484 const int sz = coeff.size();
485
486 for (int i = sz - 1; i >= 0; --i) {
487
488 if(abs(coeff[i]) < MACH_EPSILON)
489 continue;
490
491 res << (coeff[i] >= 0 ? "+ " : "- ");
492 res << abs(coeff[i]);
493
494 if(i) {
495 res << "*" << unknown << exponentiation << i;
496 res << " ";
497 }
498
499 flag = true;
500 }
501
502 if(!flag)
503 res << "0";
504
505 return res.str();
506 }
507
508
510 inline operator std::string() {
511 return to_string();
512 }
513
514
517 friend std::ostream& operator<<(std::ostream& out, const polynomial& obj) {
518 return out << obj.to_string();
519 }
520
521#endif
522
523
524 };
525
526}
527
528#endif
A polynomial of arbitrary order.
Definition polynomial.h:25
polynomial & operator/=(const polynomial &a)
Multiply a polynomial by a scalar value.
Definition polynomial.h:318
polynomial & operator+=(const polynomial &p)
Sum a polynomial to this one.
Definition polynomial.h:265
~polynomial()
Default destructor.
Definition polynomial.h:42
auto begin()
Get an iterator for the first coefficient of the polynomial.
Definition polynomial.h:365
polynomial(std::initializer_list< Type > l)
Initialize from an std::initializer_list.
Definition polynomial.h:39
polynomial(const std::vector< Type > &c)
Initialize from an std::vector.
Definition polynomial.h:36
auto cbegin()
Get a const iterator for the first coefficient of the polynomial.
Definition polynomial.h:379
EvalType eval(EvalType x) const
Evaluate the polynomial using x as variable.
Definition polynomial.h:71
polynomial operator/(const polynomial &d) const
Polynomial division.
Definition polynomial.h:185
void trim()
Remove trailing zero coefficients.
Definition polynomial.h:107
Type get(int i) const
Get the i-th by value.
Definition polynomial.h:52
EvalType operator()(EvalType x) const
Evaluate the polynomial using x as variable.
Definition polynomial.h:88
polynomial & operator*=(const polynomial &p)
Multiply two polynomials.
Definition polynomial.h:293
polynomial(Type a)
Initialize as a constant.
Definition polynomial.h:33
auto end()
Get an iterator for the one plus last coefficient of the polynomial.
Definition polynomial.h:372
vec< complex<>, 2 > quadratic_roots() const
Compute the roots of a quadratic polynomial.
Definition polynomial.h:392
friend std::ostream & operator<<(std::ostream &out, const polynomial &obj)
Stream the polynomial in string representation to an output stream (std::ostream)
Definition polynomial.h:517
polynomial & operator-=(const polynomial &p)
Subtract a polynomial from this one.
Definition polynomial.h:279
const Type & operator[](unsigned int i) const
Return the nth order coefficient.
Definition polynomial.h:58
static polynomial< Type > monomial(Type c, unsigned int order)
Returns a monomial of the given degree and coefficient.
Definition polynomial.h:444
bool operator==(const polynomial &other) const
Check whether two polynomials are equal.
Definition polynomial.h:339
polynomial operator*(const polynomial &p) const
Multiply two polynomials.
Definition polynomial.h:169
polynomial & operator/=(Type a)
Divide a polynomial by a scalar value.
Definition polynomial.h:324
Type & operator[](unsigned int i)
Return the nth order coefficient.
Definition polynomial.h:64
polynomial()
Initialize as an empty polynomial.
Definition polynomial.h:30
unsigned int find_order() const
Find the true order of the polynomial (ignoring trailing null coefficients)
Definition polynomial.h:95
auto cend()
Get a const iterator for the one plus last coefficient of the polynomial.
Definition polynomial.h:386
size_t size() const
Get the number of coefficients.
Definition polynomial.h:119
polynomial operator+(const polynomial &p) const
Sum two polynomials.
Definition polynomial.h:125
polynomial operator-(const polynomial &p) const
Subtract a polynomial from another.
Definition polynomial.h:147
polynomial operator/(Type a) const
Divide a polynomial by a scalar.
Definition polynomial.h:249
std::string to_string(const std::string &unknown="x", const std::string &exponentiation="^") const
Convert the polynomial to string representation.
Definition polynomial.h:478
polynomial & operator*=(Type a)
Multiply a polynomial by a scalar value.
Definition polynomial.h:308
polynomial operator*(Type a) const
Multiply a polynomial by a scalar.
Definition polynomial.h:237
static polynomial< Type > from_roots(const std::vector< Type > &roots)
Construct a polynomial from its roots.
Definition polynomial.h:430
Type & at(int i)
Access i-th coefficient.
Definition polynomial.h:46
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:219
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
auto min(const Vector &X)
Finds the minimum value inside a dataset.
Definition dataset.h:351
dual2 sqrt(dual2 x)
Compute the square root of a second order dual number.
Definition dual2_functions.h:54
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
auto sum(const Vector &X)
Compute the sum of a vector of real values using pairwise summation to reduce round-off error.
Definition dataset.h:219
std::vector< real > roots(RealFunction f, real a, real b, real tol=OPTIMIZATION_TOL, real steps=10)
Find the roots of a univariate real function inside a given interval, by first searching for candidat...
Definition roots.h:650
auto max(const Vector &X)
Finds the maximum value inside a dataset.
Definition dataset.h:330
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
int sgn(real x)
Return the sign of x (1 if positive, -1 if negative, 0 if null)
Definition real_analysis.h:259
dual2 square(dual2 x)
Return the square of a second order dual number.
Definition dual2_functions.h:23
real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54