Theoretica
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 const Type& at(unsigned int i) const {
47 return coeff.at(i);
48 }
49
51 inline Type& at(unsigned int i) {
52 return coeff.at(i);
53 }
54
55
57 inline const Type& operator[](unsigned int i) const {
58 return coeff[i];
59 }
60
61
63 inline Type& operator[](unsigned int i) {
64 return coeff[i];
65 }
66
67
69 template<typename EvalType = Type>
70 inline EvalType eval(EvalType x) const {
71
72 EvalType sum = EvalType(0.0);
73
74 // Evaluate using Horner's method
75 for (unsigned int i = 0; i < coeff.size(); ++i)
76 sum = coeff[coeff.size() - i - 1] + x * sum;
77
78 return sum;
79 }
80
81
83 template<typename EvalType = Type>
84 inline EvalType operator()(EvalType x) const {
85 return eval(x);
86 }
87
88
91 inline unsigned int find_order() const {
92
93 for (int i = coeff.size() - 1; i >= 0; --i) {
94 if(abs(coeff[i]) < MACH_EPSILON)
95 return i;
96 }
97
98 return 0;
99 }
100
101
103 inline void trim() {
104
105 for (unsigned int i = coeff.size() - 1; i >= 0; --i) {
106 if(abs(coeff[i]) > MACH_EPSILON)
107 break;
108
109 coeff.pop_back();
110 }
111 }
112
113
115 inline size_t size() const {
116 return coeff.size();
117 }
118
119
121 inline polynomial operator+(const polynomial& p) const {
122
123 polynomial r;
124 r.coeff.resize(max(coeff.size(), p.size()));
125
126 unsigned int i = 0;
127
128 for (; i < min(coeff.size(), p.size()); ++i)
129 r[i] = coeff[i] + p[i];
130
131 if (coeff.size() > p.size())
132 for (; i < coeff.size(); ++i)
133 r[i] = coeff[i];
134 else
135 for (; i < p.size(); ++i)
136 r[i] = p[i];
137
138 return r;
139 }
140
141
143 inline polynomial operator-(const polynomial& p) const {
144
145 polynomial r;
146 r.coeff.resize(max(coeff.size(), p.size()));
147
148 unsigned int i = 0;
149
150 for (; i < min(r.size(), p.size()); ++i)
151 r[i] = coeff[i] - p[i];
152
153 if (coeff.size() > p.size())
154 for (; i < coeff.size(); ++i)
155 r[i] = coeff[i];
156 else
157 for (; i < p.size(); ++i)
158 r[i] = -p[i];
159
160 return r;
161 }
162
163
165 inline polynomial operator*(const polynomial& p) const {
166
168 r.coeff.resize(this->size() + p.size() - 1);
169
170 for (unsigned int i = 0; i < size(); ++i) {
171 for (unsigned int j = 0; j < p.size(); ++j) {
172 r[i + j] += coeff[i] * p[j];
173 }
174 }
175
176 return r;
177 }
178
179
181 inline polynomial operator/(const polynomial& d) const {
182
183 const unsigned int d_order = d.find_order();
184 const unsigned int this_order = find_order();
185
186 if(d_order == 0 && d[0] == 0) {
187 TH_MATH_ERROR("polynomial::operator/", d[0], DIV_BY_ZERO);
188 return polynomial(nan());
189 }
190
191 // Remainder
192 polynomial r = *this;
193
194 // Quotient
196 unsigned int i = 0;
197
198 while(i < this_order) {
199
200 // Compute only once the degree of the polynomial
201 const unsigned int r_order = r.find_order();
202
203 // Stop execution if the division is complete
204 // (when the remainder is 0 or has lower degree)
205 if((r_order == 0 && (abs(r[0]) < MACH_EPSILON)) || r_order < d_order)
206 break;
207
208 // Simple division between highest degree terms
210 r[r_order] / d[d_order],
211 r_order - d_order);
212
213 // Add monomial to quotient and subtract the
214 // monomial times the dividend from the remainder
215 q += t;
216 r -= t * d;
217
218 i++;
219 }
220
221 // The algorithm has stopped iterating after a number
222 // of the dividend's degree counts
223 if(i == this_order) {
224 TH_MATH_ERROR("polynomial::operator/", i, NO_ALGO_CONVERGENCE);
225 return polynomial(nan());
226 }
227
228 return q;
229 }
230
231
233 inline polynomial operator*(Type a) const {
234
235 polynomial r = polynomial(*this);
236
237 for (unsigned int i = 0; i < size(); ++i)
238 r.coeff[i] *= a;
239
240 return r;
241 }
242
243
245 inline polynomial operator/(Type a) const {
246
247 if(abs(a) < MACH_EPSILON) {
248 TH_MATH_ERROR("polynomial::operator/", a, DIV_BY_ZERO);
249 }
250
251 polynomial r = polynomial(*this);
252
253 for (unsigned int i = 0; i < size(); ++i)
254 r.coeff[i] /= a;
255
256 return r;
257 }
258
259
261 inline polynomial& operator+=(const polynomial& p) {
262
263 // Make room for the new coefficients
264 if(coeff.size() < p.size())
265 coeff.resize(p.size(), Type(0));
266
267 for (unsigned int i = 0; i < min(size(), p.size()); ++i)
268 coeff[i] += p[i];
269
270 return *this;
271 }
272
273
275 inline polynomial& operator-=(const polynomial& p) {
276
277 // Make room for the new coefficients
278 if(coeff.size() < p.size())
279 coeff.resize(p.size(), Type(0));
280
281 for (unsigned int i = 0; i < min(coeff.size(), p.size()); ++i)
282 coeff[i] -= p[i];
283
284 return *this;
285 }
286
287
289 inline polynomial& operator*=(const polynomial& p) {
290
292 r.coeff.resize(this->size() + p.size() - 1, Type(0));
293
294 for (unsigned int i = 0; i < coeff.size(); ++i)
295 for (unsigned int j = 0; j < p.size(); ++j)
296 r[i + j] += coeff[i] * p[j];
297
298 *this = r;
299 return *this;
300 }
301
302
304 inline polynomial& operator*=(Type a) {
305
306 for (unsigned int i = 0; i < coeff.size(); ++i)
307 coeff[i] *= a;
308
309 return *this;
310 }
311
312
314 inline polynomial& operator/=(const polynomial& a) {
315 return (*this = (*this / a));
316 }
317
318
320 inline polynomial& operator/=(Type a) {
321
322 if(abs(a) < MACH_EPSILON) {
323 TH_MATH_ERROR("polynomial::operator/=", a, DIV_BY_ZERO);
324 return *this;
325 }
326
327 for (unsigned int i = 0; i < coeff.size(); ++i)
328 coeff[i] /= a;
329
330 return *this;
331 }
332
333
335 inline bool operator==(const polynomial& other) const {
336
337 const unsigned int n = min(other.size(), this->size());
338 for (unsigned int i = 0; i < n; ++i) {
339 if(abs(other.coeff[i] - coeff[i]) >= MACH_EPSILON)
340 return false;
341 }
342
343 if(size() > other.size()) {
344 for (unsigned int i = 0; i < size(); ++i) {
345 if(abs(coeff[i]) >= MACH_EPSILON)
346 return false;
347 }
348 } else {
349 for (unsigned int i = 0; i < other.size(); ++i) {
350 if(abs(other.coeff[i]) >= MACH_EPSILON)
351 return false;
352 }
353 }
354
355 return true;
356 }
357
358
361 inline auto begin() {
362 return coeff.begin();
363 }
364
365
368 inline auto end() {
369 return coeff.end();
370 }
371
372
375 inline auto cbegin() {
376 return coeff.cbegin();
377 }
378
379
382 inline auto cend() {
383 return coeff.cend();
384 }
385
386
388 inline vec<complex<>, 2> quadratic_roots() const {
389
390 const int order = find_order();
391
392 // Check that the polynomial is quadratic
393 if(order != 2) {
394 TH_MATH_ERROR("quadratic_roots", order, IMPOSSIBLE_OPERATION);
395 return vec<complex<>, 2>({nan(), nan()});
396 }
397
398 const Type p = coeff[1] / coeff[2];
399 const Type q = coeff[0] / coeff[2];
400
401 // Case when 0 is a root
402 if(abs(q) < MACH_EPSILON)
403 return {complex<>(-p), complex<>(0)};
404
405 complex<> z1, z2;
406
407 // Use Vieta's theorem to avoid catastrophic cancellation
408 if(abs(p) > 1) {
409
410 z1 = -sgn(p) * (abs(p) / 2.0
411 + abs(p) * sqrt(complex<>(0.25 - (q / p) / p)));
412 z2 = q / z1;
413
414 } else {
415
416 const complex<> s = sqrt(complex<>(0.25 * square(p) - q));
417 z1 = -p / 2.0 + s;
418 z2 = -p / 2.0 - s;
419 }
420
421 return {z1, z2};
422 }
423
424
426 inline static polynomial<Type> from_roots(
427 const std::vector<Type>& roots) {
428
429 polynomial<Type> P = {1};
430
431 // P = product((x - r_i))
432 for (unsigned int i = 0; i < roots.size(); ++i)
433 P *= polynomial<Type>({roots[i] * -1, 1});
434
435 return P;
436 }
437
438
440 inline static polynomial<Type> monomial(Type c, unsigned int order) {
441
442 polynomial m;
443 m.coeff = std::vector<Type>(order + 1, Type(0));
444 m.coeff[order] = c;
445 return m;
446 }
447
448
449 // Friend operators to enable equations of the form
450 // (Type) op. (polynomial<Type>)
451
452 inline friend polynomial<Type> operator+(
453 Type r, const polynomial<Type>& z) {
454
455 return z + polynomial(r);
456 }
457
458 inline friend polynomial<Type> operator-(
459 Type r, const polynomial<Type>& z) {
460
461 return (z * -1) + polynomial(r);
462 }
463
464 inline friend polynomial<Type> operator*(
465 Type r, const polynomial<Type>& z) {
466
467 return z * r;
468 }
469
470
471#ifndef THEORETICA_NO_PRINT
472
474 inline std::string to_string(
475 const std::string& unknown = "x",
476 const std::string& exponentiation = "^") const {
477
478 std::stringstream res;
479 bool flag = false;
480 const int sz = coeff.size();
481
482 for (int i = sz - 1; i >= 0; --i) {
483
484 if(abs(coeff[i]) < MACH_EPSILON)
485 continue;
486
487 res << (coeff[i] >= 0 ? "+ " : "- ");
488 res << abs(coeff[i]);
489
490 if(i) {
491 res << "*" << unknown << exponentiation << i;
492 res << " ";
493 }
494
495 flag = true;
496 }
497
498 if(!flag)
499 res << "0";
500
501 return res.str();
502 }
503
504
506 inline operator std::string() {
507 return to_string();
508 }
509
510
513 friend std::ostream& operator<<(std::ostream& out, const polynomial& obj) {
514 return out << obj.to_string();
515 }
516
517#endif
518
519
520 };
521
522}
523
524#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:314
polynomial & operator+=(const polynomial &p)
Sum a polynomial to this one.
Definition polynomial.h:261
~polynomial()
Default destructor.
Definition polynomial.h:42
auto begin()
Get an iterator for the first coefficient of the polynomial.
Definition polynomial.h:361
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:375
EvalType eval(EvalType x) const
Evaluate the polynomial using x as variable.
Definition polynomial.h:70
polynomial operator/(const polynomial &d) const
Polynomial division.
Definition polynomial.h:181
void trim()
Remove trailing zero coefficients.
Definition polynomial.h:103
EvalType operator()(EvalType x) const
Evaluate the polynomial using x as variable.
Definition polynomial.h:84
polynomial & operator*=(const polynomial &p)
Multiply two polynomials.
Definition polynomial.h:289
Type & at(unsigned int i)
Access i-th coefficient by reference, with bound checking.
Definition polynomial.h:51
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:368
vec< complex<>, 2 > quadratic_roots() const
Compute the roots of a quadratic polynomial.
Definition polynomial.h:388
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:513
polynomial & operator-=(const polynomial &p)
Subtract a polynomial from this one.
Definition polynomial.h:275
const Type & operator[](unsigned int i) const
Get the n-th order coefficient by constant reference.
Definition polynomial.h:57
static polynomial< Type > monomial(Type c, unsigned int order)
Returns a monomial of the given degree and coefficient.
Definition polynomial.h:440
bool operator==(const polynomial &other) const
Check whether two polynomials are equal.
Definition polynomial.h:335
polynomial operator*(const polynomial &p) const
Multiply two polynomials.
Definition polynomial.h:165
polynomial & operator/=(Type a)
Divide a polynomial by a scalar value.
Definition polynomial.h:320
const Type & at(unsigned int i) const
Get i-th coefficient by constant reference, with bound checking.
Definition polynomial.h:46
Type & operator[](unsigned int i)
Get the n-th order coefficient by reference.
Definition polynomial.h:63
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:91
auto cend()
Get a const iterator for the one plus last coefficient of the polynomial.
Definition polynomial.h:382
size_t size() const
Get the number of coefficients.
Definition polynomial.h:115
polynomial operator+(const polynomial &p) const
Sum two polynomials.
Definition polynomial.h:121
polynomial operator-(const polynomial &p) const
Subtract a polynomial from another.
Definition polynomial.h:143
polynomial operator/(Type a) const
Divide a polynomial by a scalar.
Definition polynomial.h:245
std::string to_string(const std::string &unknown="x", const std::string &exponentiation="^") const
Convert the polynomial to string representation.
Definition polynomial.h:474
polynomial & operator*=(Type a)
Multiply a polynomial by a scalar value.
Definition polynomial.h:304
polynomial operator*(Type a) const
Multiply a polynomial by a scalar.
Definition polynomial.h:233
static polynomial< Type > from_roots(const std::vector< Type > &roots)
Construct a polynomial from its roots.
Definition polynomial.h:426
#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:225
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
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
TH_CONSTEXPR real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
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