GCC Code Coverage Report

Directory: ./
File: src/2geom/polynomial.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 142 190 74.7%
Functions: 10 14 71.4%
Branches: 107 214 50.0%
Line Branch Exec Source
1 /**
2 * \file
3 * \brief Polynomial in canonical (monomial) basis
4 *//*
5 * Authors:
6 * MenTaLguY <mental@rydia.net>
7 * Krzysztof Kosiński <tweenk.pl@gmail.com>
8 * Rafał Siejakowski <rs@rs-math.net>
9 *
10 * Copyright 2007-2015 Authors
11 *
12 * This library is free software; you can redistribute it and/or
13 * modify it either under the terms of the GNU Lesser General Public
14 * License version 2.1 as published by the Free Software Foundation
15 * (the "LGPL") or, at your option, under the terms of the Mozilla
16 * Public License Version 1.1 (the "MPL"). If you do not alter this
17 * notice, a recipient may use your version of this file under either
18 * the MPL or the LGPL.
19 *
20 * You should have received a copy of the LGPL along with this library
21 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 * You should have received a copy of the MPL along with this library
24 * in the file COPYING-MPL-1.1
25 *
26 * The contents of this file are subject to the Mozilla Public License
27 * Version 1.1 (the "License"); you may not use this file except in
28 * compliance with the License. You may obtain a copy of the License at
29 * http://www.mozilla.org/MPL/
30 *
31 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
32 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
33 * the specific language governing rights and limitations.
34 */
35
36 #include <algorithm>
37 #include <2geom/polynomial.h>
38 #include <2geom/math-utils.h>
39 #include <math.h>
40
41 #ifdef HAVE_GSL
42 #include <gsl/gsl_poly.h>
43 #endif
44
45 namespace Geom {
46
47 #ifndef M_PI
48 # define M_PI 3.14159265358979323846
49 #endif
50
51 122 Poly Poly::operator*(const Poly& p) const {
52 122 Poly result;
53
1/2
✓ Branch 3 taken 122 times.
✗ Branch 4 not taken.
 122 result.resize(degree() + p.degree()+1);
54
55
2/2
✓ Branch 1 taken 341 times.
✓ Branch 2 taken 122 times.
 463 for(unsigned i = 0; i < size(); i++) {
56
2/2
✓ Branch 1 taken 2624 times.
✓ Branch 2 taken 341 times.
 2965 for(unsigned j = 0; j < p.size(); j++) {
57 2624 result[i+j] += (*this)[i] * p[j];
58 }
59 }
60 122 return result;
61 ✗ }
62
63 /*double Poly::eval(double x) const {
64 return gsl_poly_eval(&coeff[0], size(), x);
65 }*/
66
67 51 void Poly::normalize() {
68
2/2
✓ Branch 1 taken 77 times.
✓ Branch 2 taken 51 times.
 128 while(back() == 0)
69 77 pop_back();
70 51 }
71
72 ✗ void Poly::monicify() {
73 ✗ normalize();
74
75 ✗ double scale = 1./back(); // unitize
76
77 ✗ for(unsigned i = 0; i < size(); i++) {
78 ✗ (*this)[i] *= scale;
79 }
80 ✗ }
81
82
83 #ifdef HAVE_GSL
84 30 std::vector<std::complex<double> > solve(Poly const & pp) {
85
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 Poly p(pp);
86
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 p.normalize();
87 gsl_poly_complex_workspace * w
88
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 = gsl_poly_complex_workspace_alloc (p.size());
89
90
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 30 times.
✗ Branch 6 not taken.
 30 gsl_complex_packed_ptr z = new double[p.degree()*2];
91
2/4
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 30 times.
✗ Branch 6 not taken.
 30 double* a = new double[p.size()];
92
2/2
✓ Branch 1 taken 150 times.
✓ Branch 2 taken 30 times.
 180 for(unsigned int i = 0; i < p.size(); i++)
93 150 a[i] = p[i];
94 30 std::vector<std::complex<double> > roots;
95 //roots.resize(p.degree());
96
97
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 gsl_poly_complex_solve (a, p.size(), w, z);
98
1/2
✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
 30 delete[]a;
99
100
1/2
✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
 30 gsl_poly_complex_workspace_free (w);
101
102
2/2
✓ Branch 1 taken 120 times.
✓ Branch 2 taken 30 times.
 150 for (unsigned int i = 0; i < p.degree(); i++) {
103
1/2
✓ Branch 1 taken 120 times.
✗ Branch 2 not taken.
 120 roots.emplace_back(z[2*i] ,z[2*i+1]);
104 //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
105 }
106
1/2
✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
 30 delete[] z;
107 60 return roots;
108 30 }
109
110 30 std::vector<double > solve_reals(Poly const & p) {
111
1/2
✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
 30 std::vector<std::complex<double> > roots = solve(p);
112 30 std::vector<double> real_roots;
113
114
2/2
✓ Branch 7 taken 120 times.
✓ Branch 8 taken 30 times.
 150 for(auto & root : roots) {
115
2/2
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 50 times.
 120 if(root.imag() == 0) // should be more lenient perhaps
116
1/2
✓ Branch 3 taken 70 times.
✗ Branch 4 not taken.
 70 real_roots.push_back(root.real());
117 }
118 30 return real_roots;
119 30 }
120 #endif
121
122 ✗ double polish_root(Poly const & p, double guess, double tol) {
123 ✗ Poly dp = derivative(p);
124
125 ✗ double fn = p(guess);
126 ✗ while(fabs(fn) > tol) {
127 ✗ guess -= fn/dp(guess);
128 ✗ fn = p(guess);
129 }
130 ✗ return guess;
131 ✗ }
132
133 3 Poly integral(Poly const & p) {
134 3 Poly result;
135
136
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 result.reserve(p.size()+1);
137
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 result.push_back(0); // arbitrary const
138
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 3 times.
 15 for(unsigned i = 0; i < p.size(); i++) {
139
1/2
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
 12 result.push_back(p[i]/(i+1));
140 }
141 3 return result;
142
143 ✗ }
144
145 1 Poly derivative(Poly const & p) {
146 1 Poly result;
147
148
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if(p.size() <= 1)
149 ✗ return Poly(0);
150
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 result.reserve(p.size()-1);
151
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 1 times.
 4 for(unsigned i = 1; i < p.size(); i++) {
152
1/2
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 3 result.push_back(i*p[i]);
153 }
154
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 return result;
155 1 }
156
157 3 Poly compose(Poly const & a, Poly const & b) {
158 3 Poly result;
159
160
2/2
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 3 times.
 13 for(unsigned i = a.size(); i > 0; i--) {
161
4/8
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 10 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 10 times.
✗ Branch 15 not taken.
 10 result = Poly(a[i-1]) + result * b;
162 }
163 3 return result;
164
165 ✗ }
166
167 /* This version is backwards - dividing taylor terms
168 Poly divide(Poly const &a, Poly const &b, Poly &r) {
169 Poly c;
170 r = a; // remainder
171
172 const unsigned k = a.size();
173 r.resize(k, 0);
174 c.resize(k, 0);
175
176 for(unsigned i = 0; i < k; i++) {
177 double ci = r[i]/b[0];
178 c[i] += ci;
179 Poly bb = ci*b;
180 std::cout << ci <<"*" << b << ", r= " << r << std::endl;
181 r -= bb.shifted(i);
182 }
183
184 return c;
185 }
186 */
187
188 ✗ Poly divide(Poly const &a, Poly const &b, Poly &r) {
189 ✗ Poly c;
190 ✗ r = a; // remainder
191 ✗ assert(b.size() > 0);
192
193 ✗ const unsigned k = a.degree();
194 ✗ const unsigned l = b.degree();
195 ✗ c.resize(k, 0.);
196
197 ✗ for(unsigned i = k; i >= l; i--) {
198 //assert(i >= 0);
199 ✗ double ci = r.back()/b.back();
200 ✗ c[i-l] += ci;
201 ✗ Poly bb = ci*b;
202 //std::cout << ci <<"*(" << b.shifted(i-l) << ") = "
203 // << bb.shifted(i-l) << " r= " << r << std::endl;
204 ✗ r -= bb.shifted(i-l);
205 ✗ r.pop_back();
206 ✗ }
207 //std::cout << "r= " << r << std::endl;
208 ✗ r.normalize();
209 ✗ c.normalize();
210
211 ✗ return c;
212 ✗ }
213
214 ✗ Poly gcd(Poly const &a, Poly const &b, const double /*tol*/) {
215 ✗ if(a.size() < b.size())
216 ✗ return gcd(b, a);
217 ✗ if(b.size() <= 0)
218 ✗ return a;
219 ✗ if(b.size() == 1)
220 ✗ return a;
221 ✗ Poly r;
222 ✗ divide(a, b, r);
223 ✗ return gcd(b, r);
224 ✗ }
225
226 214560 std::vector<Coord> solve_quadratic(Coord a, Coord b, Coord c)
227 {
228 214560 std::vector<Coord> result;
229
230
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 214559 times.
 214560 if (a == 0) {
231 // linear equation
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (b == 0) return result;
233
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 result.push_back(-c/b);
234 1 return result;
235 }
236
237 214559 Coord delta = b*b - 4*a*c;
238
239
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 214483 times.
 214559 if (delta == 0) {
240 // one root
241
1/2
✓ Branch 2 taken 76 times.
✗ Branch 3 not taken.
 76 result.push_back(-b / (2*a));
242
2/2
✓ Branch 0 taken 160581 times.
✓ Branch 1 taken 53902 times.
 214483 } else if (delta > 0) {
243 // two roots
244 160581 Coord delta_sqrt = sqrt(delta);
245
246 // Use different formulas depending on sign of b to preserve
247 // numerical stability. See e.g.:
248 // http://people.csail.mit.edu/bkph/articles/Quadratics.pdf
249
2/2
✓ Branch 0 taken 77113 times.
✓ Branch 1 taken 83468 times.
 160581 int sign = b >= 0 ? 1 : -1;
250 160581 Coord t = -0.5 * (b + sign * delta_sqrt);
251
1/2
✓ Branch 2 taken 160581 times.
✗ Branch 3 not taken.
 160581 result.push_back(t / a);
252
1/2
✓ Branch 2 taken 160581 times.
✗ Branch 3 not taken.
 160581 result.push_back(c / t);
253 }
254 // no roots otherwise
255
256
1/2
✓ Branch 3 taken 214559 times.
✗ Branch 4 not taken.
 214559 std::sort(result.begin(), result.end());
257 214559 return result;
258 ✗ }
259
260 112032 std::vector<Coord> solve_cubic(Coord a, Coord b, Coord c, Coord d)
261 {
262 // based on:
263 // http://mathworld.wolfram.com/CubicFormula.html
264
265
2/2
✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 102032 times.
 112032 if (a == 0) {
266
1/2
✓ Branch 1 taken 10000 times.
✗ Branch 2 not taken.
 10000 return solve_quadratic(b, c, d);
267 }
268
2/2
✓ Branch 0 taken 10011 times.
✓ Branch 1 taken 92021 times.
 102032 if (d == 0) {
269 // divide by x
270
1/2
✓ Branch 2 taken 10011 times.
✗ Branch 3 not taken.
 10011 std::vector<Coord> result = solve_quadratic(a, b, c);
271
1/2
✓ Branch 2 taken 10011 times.
✗ Branch 3 not taken.
 10011 result.push_back(0);
272
1/2
✓ Branch 3 taken 10011 times.
✗ Branch 4 not taken.
 10011 std::sort(result.begin(), result.end());
273 10011 return result;
274 10011 }
275
276 92021 std::vector<Coord> result;
277
278 // 1. divide everything by a to bring to canonical form
279 92021 b /= a;
280 92021 c /= a;
281 92021 d /= a;
282
283 // 2. eliminate x^2 term: x^3 + 3Qx - 2R = 0
284 92021 Coord Q = (3*c - b*b) / 9;
285 92021 Coord R = (-27 * d + b * (9*c - 2*b*b)) / 54;
286
287 // 3. compute polynomial discriminant
288 92021 Coord D = Q*Q*Q + R*R;
289 92021 Coord term1 = b/3;
290
291
2/2
✓ Branch 0 taken 56787 times.
✓ Branch 1 taken 35234 times.
 92021 if (D > 0) {
292 // only one real root
293 56787 Coord S = cbrt(R + sqrt(D));
294 56787 Coord T = cbrt(R - sqrt(D));
295
1/2
✓ Branch 2 taken 56787 times.
✗ Branch 3 not taken.
 56787 result.push_back(-b/3 + S + T);
296
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 35176 times.
 35234 } else if (D == 0) {
297 // 3 real roots, 2 of which are equal
298 58 Coord rroot = cbrt(R);
299
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 result.reserve(3);
300
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 58 result.push_back(-term1 + 2*rroot);
301
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 58 result.push_back(-term1 - rroot);
302
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 58 result.push_back(-term1 - rroot);
303 } else {
304 // 3 distinct real roots
305
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 35176 times.
 35176 assert(Q < 0);
306 35176 Coord theta = acos(R / sqrt(-Q*Q*Q));
307 35176 Coord rroot = 2 * sqrt(-Q);
308
1/2
✓ Branch 1 taken 35176 times.
✗ Branch 2 not taken.
 35176 result.reserve(3);
309
1/2
✓ Branch 2 taken 35176 times.
✗ Branch 3 not taken.
 35176 result.push_back(-term1 + rroot * cos(theta / 3));
310
1/2
✓ Branch 2 taken 35176 times.
✗ Branch 3 not taken.
 35176 result.push_back(-term1 + rroot * cos((theta + 2*M_PI) / 3));
311
1/2
✓ Branch 2 taken 35176 times.
✗ Branch 3 not taken.
 35176 result.push_back(-term1 + rroot * cos((theta + 4*M_PI) / 3));
312 }
313
314
1/2
✓ Branch 3 taken 92021 times.
✗ Branch 4 not taken.
 92021 std::sort(result.begin(), result.end());
315 92021 return result;
316 92021 }
317
318 81000 std::vector<Coord> solve_quartic(Coord a, Coord b, Coord c, Coord d, Coord e)
319 {
320 // Based on a variation of the Ferrari-Lagrange method, see
321 // "A universal method of solving quartic equations" by S. Shmakov,
322 // International Journal of Pure and Applied Mathematics vol. 71 no. 2.
323
324
2/2
✓ Branch 0 taken 20000 times.
✓ Branch 1 taken 61000 times.
 81000 if (a == 0) {
325
1/2
✓ Branch 1 taken 20000 times.
✗ Branch 2 not taken.
 20000 return solve_cubic(b, c, d, e);
326 }
327
2/2
✓ Branch 0 taken 20000 times.
✓ Branch 1 taken 41000 times.
 61000 if (e == 0) { // Divide by x
328
1/2
✓ Branch 2 taken 20000 times.
✗ Branch 3 not taken.
 20000 auto result = solve_cubic(a, b, c, d);
329
1/2
✓ Branch 2 taken 20000 times.
✗ Branch 3 not taken.
 20000 result.push_back(0);
330
1/2
✓ Branch 3 taken 20000 times.
✗ Branch 4 not taken.
 20000 std::sort(result.begin(), result.end());
331 20000 return result;
332 20000 }
333
334 // Divide out by a so that the leading coefficient is 1.
335 41000 b /= a;
336 41000 c /= a;
337 41000 d /= a;
338 41000 e /= a;
339
340 // Solve the resolvent cubic
341
1/2
✓ Branch 4 taken 41000 times.
✗ Branch 5 not taken.
 41000 auto const resolvent_solutions = solve_cubic(1, -c, b * d - 4 * e, 4 * c * e - sqr(b) * e - sqr(d));
342 // If there are 3 solutions, pick the middle one, else the first one.
343 41000 auto const y = resolvent_solutions[resolvent_solutions.size() == 3];
344
345 // Find the quadratic factors
346
1/2
✓ Branch 2 taken 41000 times.
✗ Branch 3 not taken.
 41000 auto linear_terms = solve_quadratic(1, -b, c - y);
347
1/2
✓ Branch 2 taken 41000 times.
✗ Branch 3 not taken.
 41000 auto constant_terms = solve_quadratic(1, -y, e);
348
5/6
✓ Branch 1 taken 30742 times.
✓ Branch 2 taken 10258 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 30742 times.
✓ Branch 6 taken 10258 times.
✓ Branch 7 taken 30742 times.
 41000 if (linear_terms.size() < 2 || constant_terms.size() < 2) {
349 10258 return {}; // There are no roots
350 }
351
352 {
353 // Reorder constant terms if needed so that they correspond to linear terms
354 30742 auto const current_cross = linear_terms[0] * constant_terms[1] + linear_terms[1] * constant_terms[0];
355 30742 auto const reordered_cross = linear_terms[0] * constant_terms[0] + linear_terms[1] * constant_terms[1];
356
2/2
✓ Branch 2 taken 15466 times.
✓ Branch 3 taken 15276 times.
 30742 if (std::abs(d - reordered_cross) < std::abs(d - current_cross)) {
357 15466 std::swap(constant_terms[0], constant_terms[1]);
358 }
359 }
360
361 30742 std::vector<Coord> result;
362
1/2
✓ Branch 1 taken 30742 times.
✗ Branch 2 not taken.
 30742 result.reserve(4);
363
2/2
✓ Branch 4 taken 61484 times.
✓ Branch 5 taken 30742 times.
 92226 for (size_t i : {0, 1}) {
364
1/2
✓ Branch 4 taken 61484 times.
✗ Branch 5 not taken.
 61484 auto const factor_roots = solve_quadratic(1, linear_terms[i], constant_terms[i]);
365
1/2
✓ Branch 7 taken 61484 times.
✗ Branch 8 not taken.
 61484 result.insert(result.end(), factor_roots.begin(), factor_roots.end());
366 61484 }
367
1/2
✓ Branch 3 taken 30742 times.
✗ Branch 4 not taken.
 30742 std::sort(result.begin(), result.end());
368 30742 return result;
369 41000 }
370
371 } //namespace Geom
372
373 /*
374 Local Variables:
375 mode:c++
376 c-file-style:"stroustrup"
377 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
378 indent-tabs-mode:nil
379 fill-column:99
380 End:
381 */
382 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
383