GCC Code Coverage Report

Directory: ./
File: src/2geom/sbasis.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 288 323 89.2%
Functions: 24 26 92.3%
Branches: 290 506 57.3%
Line Branch Exec Source
1 /*
2 * sbasis.cpp - S-power basis function class + supporting classes
3 *
4 * Authors:
5 * Nathan Hurst <njh@mail.csse.monash.edu.au>
6 * Michael Sloan <mgsloan@gmail.com>
7 *
8 * Copyright (C) 2006-2007 authors
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it either under the terms of the GNU Lesser General Public
12 * License version 2.1 as published by the Free Software Foundation
13 * (the "LGPL") or, at your option, under the terms of the Mozilla
14 * Public License Version 1.1 (the "MPL"). If you do not alter this
15 * notice, a recipient may use your version of this file under either
16 * the MPL or the LGPL.
17 *
18 * You should have received a copy of the LGPL along with this library
19 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * You should have received a copy of the MPL along with this library
22 * in the file COPYING-MPL-1.1
23 *
24 * The contents of this file are subject to the Mozilla Public License
25 * Version 1.1 (the "License"); you may not use this file except in
26 * compliance with the License. You may obtain a copy of the License at
27 * http://www.mozilla.org/MPL/
28 *
29 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31 * the specific language governing rights and limitations.
32 */
33
34 #include <cmath>
35
36 #include <2geom/sbasis.h>
37 #include <2geom/math-utils.h>
38
39 namespace Geom {
40
41 #ifndef M_PI
42 # define M_PI 3.14159265358979323846
43 #endif
44
45 /** bound the error from term truncation
46 \param tail first term to chop
47 \returns the largest possible error this truncation could give
48 */
49 55246 double SBasis::tailError(unsigned tail) const {
50
1/2
✓ Branch 3 taken 55246 times.
✗ Branch 4 not taken.
 55246 Interval bs = *bounds_fast(*this, tail);
51 110492 return std::max(fabs(bs.min()),fabs(bs.max()));
52 }
53
54 /** test all coefficients are finite
55 */
56 1 bool SBasis::isFinite() const {
57
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
 3 for(unsigned i = 0; i < size(); i++) {
58
1/2
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 if(!(*this)[i].isFinite())
59 ✗ return false;
60 }
61 1 return true;
62 }
63
64 /** Compute the value and the first n derivatives
65 \param t position to evaluate
66 \param n number of derivatives (not counting value)
67 \returns a vector with the value and the n derivative evaluations
68
69 There is an elegant way to compute the value and n derivatives for a polynomial using a variant of horner's rule. Someone will someday work out how for sbasis.
70 */
71 1 std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
72
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 3 std::vector<double> ret(n+1);
73 1 ret[0] = valueAt(t);
74
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 SBasis tmp = *this;
75
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 times.
 6 for(unsigned i = 1; i < n+1; i++) {
76
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 tmp.derive();
77 5 ret[i] = tmp.valueAt(t);
78 }
79 2 return ret;
80 1 }
81
82
83 /** Compute the pointwise sum of a and b (Exact)
84 \param a,b sbasis functions
85 \returns sbasis equal to a+b
86
87 */
88 586958 SBasis operator+(const SBasis& a, const SBasis& b) {
89 586958 const unsigned out_size = std::max(a.size(), b.size());
90 586958 const unsigned min_size = std::min(a.size(), b.size());
91
1/2
✓ Branch 3 taken 586958 times.
✗ Branch 4 not taken.
 586958 SBasis result(out_size, Linear());
92
93
2/2
✓ Branch 0 taken 1035446 times.
✓ Branch 1 taken 586958 times.
 1622404 for(unsigned i = 0; i < min_size; i++) {
94
1/2
✓ Branch 6 taken 1035446 times.
✗ Branch 7 not taken.
 1035446 result[i] = a[i] + b[i];
95 }
96
2/2
✓ Branch 1 taken 41716 times.
✓ Branch 2 taken 586958 times.
 628674 for(unsigned i = min_size; i < a.size(); i++)
97
1/2
✓ Branch 2 taken 41716 times.
✗ Branch 3 not taken.
 41716 result[i] = a[i];
98
2/2
✓ Branch 1 taken 8437 times.
✓ Branch 2 taken 586958 times.
 595395 for(unsigned i = min_size; i < b.size(); i++)
99
1/2
✓ Branch 2 taken 8437 times.
✗ Branch 3 not taken.
 8437 result[i] = b[i];
100
101
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 586958 times.
 586958 assert(result.size() == out_size);
102 586958 return result;
103 ✗ }
104
105 /** Compute the pointwise difference of a and b (Exact)
106 \param a,b sbasis functions
107 \returns sbasis equal to a-b
108
109 */
110 948258 SBasis operator-(const SBasis& a, const SBasis& b) {
111 948258 const unsigned out_size = std::max(a.size(), b.size());
112 948258 const unsigned min_size = std::min(a.size(), b.size());
113
1/2
✓ Branch 3 taken 948258 times.
✗ Branch 4 not taken.
 948258 SBasis result(out_size, Linear());
114
115
2/2
✓ Branch 0 taken 1083384 times.
✓ Branch 1 taken 948258 times.
 2031642 for(unsigned i = 0; i < min_size; i++) {
116
1/2
✓ Branch 6 taken 1083384 times.
✗ Branch 7 not taken.
 1083384 result[i] = a[i] - b[i];
117 }
118
2/2
✓ Branch 1 taken 33017 times.
✓ Branch 2 taken 948258 times.
 981275 for(unsigned i = min_size; i < a.size(); i++)
119
1/2
✓ Branch 2 taken 33017 times.
✗ Branch 3 not taken.
 33017 result[i] = a[i];
120
2/2
✓ Branch 1 taken 498650 times.
✓ Branch 2 taken 948258 times.
 1446908 for(unsigned i = min_size; i < b.size(); i++)
121
1/2
✓ Branch 4 taken 498650 times.
✗ Branch 5 not taken.
 498650 result[i] = -b[i];
122
123
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 948258 times.
 948258 assert(result.size() == out_size);
124 948258 return result;
125 ✗ }
126
127 /** Compute the pointwise sum of a and b and store in a (Exact)
128 \param a,b sbasis functions
129 \returns sbasis equal to a+b
130
131 */
132 19358 SBasis& operator+=(SBasis& a, const SBasis& b) {
133 19358 const unsigned out_size = std::max(a.size(), b.size());
134 19358 const unsigned min_size = std::min(a.size(), b.size());
135 19358 a.resize(out_size);
136
137
2/2
✓ Branch 0 taken 38917 times.
✓ Branch 1 taken 19358 times.
 58275 for(unsigned i = 0; i < min_size; i++)
138
1/2
✓ Branch 3 taken 38917 times.
✗ Branch 4 not taken.
 38917 a[i] += b[i];
139
2/2
✓ Branch 1 taken 19439 times.
✓ Branch 2 taken 19358 times.
 38797 for(unsigned i = min_size; i < b.size(); i++)
140 19439 a[i] = b[i];
141
142
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 19358 times.
 19358 assert(a.size() == out_size);
143 19358 return a;
144 }
145
146 /** Compute the pointwise difference of a and b and store in a (Exact)
147 \param a,b sbasis functions
148 \returns sbasis equal to a-b
149
150 */
151 41446 SBasis& operator-=(SBasis& a, const SBasis& b) {
152 41446 const unsigned out_size = std::max(a.size(), b.size());
153 41446 const unsigned min_size = std::min(a.size(), b.size());
154 41446 a.resize(out_size);
155
156
2/2
✓ Branch 0 taken 147622 times.
✓ Branch 1 taken 41446 times.
 189068 for(unsigned i = 0; i < min_size; i++)
157
1/2
✓ Branch 3 taken 147622 times.
✗ Branch 4 not taken.
 147622 a[i] -= b[i];
158
2/2
✓ Branch 1 taken 90802 times.
✓ Branch 2 taken 41446 times.
 132248 for(unsigned i = min_size; i < b.size(); i++)
159
1/2
✓ Branch 4 taken 90802 times.
✗ Branch 5 not taken.
 90802 a[i] = -b[i];
160
161
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 41446 times.
 41446 assert(a.size() == out_size);
162 41446 return a;
163 }
164
165 /** Compute the pointwise product of a and b (Exact)
166 \param a,b sbasis functions
167 \returns sbasis equal to a*b
168
169 */
170 1196326 SBasis operator*(SBasis const &a, double k) {
171
1/2
✓ Branch 4 taken 1196326 times.
✗ Branch 5 not taken.
 1196326 SBasis c(a.size(), Linear());
172
2/2
✓ Branch 1 taken 1924525 times.
✓ Branch 2 taken 1196326 times.
 3120851 for(unsigned i = 0; i < a.size(); i++)
173
1/2
✓ Branch 4 taken 1924525 times.
✗ Branch 5 not taken.
 1924525 c[i] = a[i] * k;
174 1196326 return c;
175 ✗ }
176
177 /** Compute the pointwise product of a and b and store the value in a (Exact)
178 \param a,b sbasis functions
179 \returns sbasis equal to a*b
180
181 */
182 47409 SBasis& operator*=(SBasis& a, double b) {
183
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 47409 times.
 47409 if (a.isZero()) return a;
184
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 47409 times.
 47409 if (b == 0)
185 ✗ a.clear();
186 else
187
2/2
✓ Branch 7 taken 165023 times.
✓ Branch 8 taken 47409 times.
 212432 for(auto & i : a)
188 165023 i *= b;
189 47409 return a;
190 }
191
192 /** multiply a by x^sh in place (Exact)
193 \param a sbasis function
194 \param sh power
195 \returns a
196
197 */
198 41434 SBasis shift(SBasis const &a, int sh) {
199 41434 size_t n = a.size()+sh;
200
1/2
✓ Branch 3 taken 41434 times.
✗ Branch 4 not taken.
 41434 SBasis c(n, Linear());
201 41434 size_t m = std::max(0, sh);
202
203
2/2
✓ Branch 0 taken 82942 times.
✓ Branch 1 taken 41434 times.
 124376 for(int i = 0; i < sh; i++)
204
1/2
✓ Branch 3 taken 82942 times.
✗ Branch 4 not taken.
 82942 c[i] = Linear(0,0);
205
2/2
✓ Branch 5 taken 165915 times.
✓ Branch 6 taken 41434 times.
 207349 for(size_t i = m, j = std::max(0,-sh); i < n; i++, j++)
206
1/2
✓ Branch 2 taken 165915 times.
✗ Branch 3 not taken.
 165915 c[i] = a[j];
207 41434 return c;
208 ✗ }
209
210 /** multiply a by x^sh (Exact)
211 \param a linear function
212 \param sh power
213 \returns a* x^sh
214
215 */
216 41426 SBasis shift(Linear const &a, int sh) {
217 41426 size_t n = 1+sh;
218
1/2
✓ Branch 3 taken 41426 times.
✗ Branch 4 not taken.
 41426 SBasis c(n, Linear());
219
220
2/2
✓ Branch 0 taken 82926 times.
✓ Branch 1 taken 41426 times.
 124352 for(int i = 0; i < sh; i++)
221
1/2
✓ Branch 3 taken 82926 times.
✗ Branch 4 not taken.
 82926 c[i] = Linear(0,0);
222
1/2
✓ Branch 0 taken 41426 times.
✗ Branch 1 not taken.
 41426 if(sh >= 0)
223
1/2
✓ Branch 1 taken 41426 times.
✗ Branch 2 not taken.
 41426 c[sh] = a;
224 41426 return c;
225 ✗ }
226
227 #if 0
228 SBasis multiply(SBasis const &a, SBasis const &b) {
229 // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)}
230
231 // shift(1, a.Tri*b.Tri)
232 SBasis c(a.size() + b.size(), Linear(0,0));
233 if(a.isZero() || b.isZero())
234 return c;
235 for(unsigned j = 0; j < b.size(); j++) {
236 for(unsigned i = j; i < a.size()+j; i++) {
237 double tri = b[j].tri()*a[i-j].tri();
238 c[i+1/*shift*/] += Linear(-tri);
239 }
240 }
241 for(unsigned j = 0; j < b.size(); j++) {
242 for(unsigned i = j; i < a.size()+j; i++) {
243 for(unsigned dim = 0; dim < 2; dim++)
244 c[i][dim] += b[j][dim]*a[i-j][dim];
245 }
246 }
247 c.normalize();
248 //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
249 return c;
250 }
251 #else
252
253 /** Compute the pointwise product of a and b adding c (Exact)
254 \param a,b,c sbasis functions
255 \returns sbasis equal to a*b+c
256
257 The added term is almost free
258 */
259 992074 SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) {
260
6/6
✓ Branch 1 taken 700605 times.
✓ Branch 2 taken 291469 times.
✓ Branch 4 taken 27538 times.
✓ Branch 5 taken 673067 times.
✓ Branch 6 taken 319007 times.
✓ Branch 7 taken 673067 times.
 992074 if(a.isZero() || b.isZero())
261 319007 return c;
262
1/2
✓ Branch 5 taken 673067 times.
✗ Branch 6 not taken.
 673067 c.resize(a.size() + b.size(), Linear(0,0));
263
2/2
✓ Branch 1 taken 2189064 times.
✓ Branch 2 taken 673067 times.
 2862131 for(unsigned j = 0; j < b.size(); j++) {
264
2/2
✓ Branch 1 taken 4779452 times.
✓ Branch 2 taken 2189064 times.
 6968516 for(unsigned i = j; i < a.size()+j; i++) {
265 4779452 double tri = b[j].tri()*a[i-j].tri();
266
1/2
✓ Branch 3 taken 4779452 times.
✗ Branch 4 not taken.
 4779452 c[i+1/*shift*/] += Linear(-tri);
267 }
268 }
269
2/2
✓ Branch 1 taken 2189064 times.
✓ Branch 2 taken 673067 times.
 2862131 for(unsigned j = 0; j < b.size(); j++) {
270
2/2
✓ Branch 1 taken 4779452 times.
✓ Branch 2 taken 2189064 times.
 6968516 for(unsigned i = j; i < a.size()+j; i++) {
271
2/2
✓ Branch 0 taken 9558904 times.
✓ Branch 1 taken 4779452 times.
 14338356 for(unsigned dim = 0; dim < 2; dim++)
272
1/2
✓ Branch 7 taken 9558904 times.
✗ Branch 8 not taken.
 9558904 c[i][dim] += b[j][dim]*a[i-j][dim];
273 }
274 }
275 673067 c.normalize();
276 //assert(!(0 == c.back()[0] && 0 == c.back()[1]));
277 673067 return c;
278 }
279
280 /** Compute the pointwise product of a and b (Exact)
281 \param a,b sbasis functions
282 \returns sbasis equal to a*b
283
284 */
285 538630 SBasis multiply(SBasis const &a, SBasis const &b) {
286
6/6
✓ Branch 1 taken 532585 times.
✓ Branch 2 taken 6045 times.
✓ Branch 4 taken 43930 times.
✓ Branch 5 taken 488655 times.
✓ Branch 6 taken 49975 times.
✓ Branch 7 taken 488655 times.
 538630 if(a.isZero() || b.isZero()) {
287
1/2
✓ Branch 4 taken 49975 times.
✗ Branch 5 not taken.
 49975 SBasis c(1, Linear(0,0));
288
1/2
✓ Branch 1 taken 49975 times.
✗ Branch 2 not taken.
 49975 return c;
289 49975 }
290
1/2
✓ Branch 6 taken 488655 times.
✗ Branch 7 not taken.
 488655 SBasis c(a.size() + b.size(), Linear(0,0));
291
2/4
✓ Branch 2 taken 488655 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 488655 times.
✗ Branch 6 not taken.
 488655 return multiply_add(a, b, c);
292 488655 }
293 #endif
294 /** Compute the integral of a (Exact)
295 \param a sbasis functions
296 \returns sbasis integral(a)
297
298 */
299 11902 SBasis integral(SBasis const &c) {
300 11902 SBasis a;
301
1/2
✓ Branch 4 taken 11902 times.
✗ Branch 5 not taken.
 11902 a.resize(c.size() + 1, Linear(0,0));
302
1/2
✓ Branch 3 taken 11902 times.
✗ Branch 4 not taken.
 11902 a[0] = Linear(0,0);
303
304
2/2
✓ Branch 1 taken 48006 times.
✓ Branch 2 taken 11902 times.
 59908 for(unsigned k = 1; k < c.size() + 1; k++) {
305 48006 double ahat = -c[k-1].tri()/(2*k);
306
2/4
✓ Branch 1 taken 48006 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 48006 times.
✗ Branch 6 not taken.
 48006 a[k][0] = a[k][1] = ahat;
307 }
308 11902 double aTri = 0;
309
2/2
✓ Branch 1 taken 48006 times.
✓ Branch 2 taken 11902 times.
 59908 for(int k = c.size()-1; k >= 0; k--) {
310 48006 aTri = (c[k].hat() + (k+1)*aTri/2)/(2*k+1);
311
1/2
✓ Branch 1 taken 48006 times.
✗ Branch 2 not taken.
 48006 a[k][0] -= aTri/2;
312
1/2
✓ Branch 1 taken 48006 times.
✗ Branch 2 not taken.
 48006 a[k][1] += aTri/2;
313 }
314 11902 a.normalize();
315 11902 return a;
316 ✗ }
317
318 /** Compute the derivative of a (Exact)
319 \param a sbasis functions
320 \returns sbasis da/dt
321
322 */
323 90736 SBasis derivative(SBasis const &a) {
324 90736 SBasis c;
325
1/2
✓ Branch 4 taken 90736 times.
✗ Branch 5 not taken.
 90736 c.resize(a.size(), Linear(0,0));
326
2/2
✓ Branch 1 taken 13 times.
✓ Branch 2 taken 90723 times.
 90736 if(a.isZero())
327 13 return c;
328
329
2/2
✓ Branch 1 taken 90844 times.
✓ Branch 2 taken 90723 times.
 181567 for(unsigned k = 0; k < a.size()-1; k++) {
330 90844 double d = (2*k+1)*(a[k][1] - a[k][0]);
331
332
1/2
✓ Branch 4 taken 90844 times.
✗ Branch 5 not taken.
 90844 c[k][0] = d + (k+1)*a[k+1][0];
333
1/2
✓ Branch 4 taken 90844 times.
✗ Branch 5 not taken.
 90844 c[k][1] = d - (k+1)*a[k+1][1];
334 }
335 90723 int k = a.size()-1;
336 90723 double d = (2*k+1)*(a[k][1] - a[k][0]);
337
4/4
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 90710 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 7 times.
 90723 if (d == 0 && k > 0) {
338 6 c.pop_back();
339 } else {
340
1/2
✓ Branch 1 taken 90717 times.
✗ Branch 2 not taken.
 90717 c[k][0] = d;
341
1/2
✓ Branch 1 taken 90717 times.
✗ Branch 2 not taken.
 90717 c[k][1] = d;
342 }
343
344 90723 return c;
345 ✗ }
346
347 /** Compute the derivative of this inplace (Exact)
348
349 */
350 5 void SBasis::derive() { // in place version
351
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 4 times.
 5 if(isZero()) return;
352
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
 6 for(unsigned k = 0; k < size()-1; k++) {
353 2 double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
354
355 2 (*this)[k][0] = d + (k+1)*(*this)[k+1][0];
356 2 (*this)[k][1] = d - (k+1)*(*this)[k+1][1];
357 }
358 4 int k = size()-1;
359 4 double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]);
360
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 1 times.
 4 if (d == 0 && k > 0) {
361 1 pop_back();
362 } else {
363 3 (*this)[k][0] = d;
364 3 (*this)[k][1] = d;
365 }
366 }
367
368 /** Compute the sqrt of a
369 \param a sbasis functions
370 \returns sbasis \f[ \sqrt{a} \f]
371
372 It is recommended to use the piecewise version unless you have good reason.
373 TODO: convert int k to unsigned k, and remove cast
374 */
375 4 SBasis sqrt(SBasis const &a, int k) {
376
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 SBasis c;
377
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4 times.
 4 if(a.isZero() || k == 0)
378 ✗ return c;
379
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 c.resize(k, Linear(0,0));
380
1/2
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
 4 c[0] = Linear(std::sqrt(a[0][0]), std::sqrt(a[0][1]));
381
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 4 SBasis r = a - multiply(c, c); // remainder
382
383
3/6
✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 26 times.
✗ Branch 6 not taken.
 26 for(unsigned i = 1; i <= (unsigned)k && i<r.size(); i++) {
384
4/8
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 26 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 26 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 26 times.
✗ Branch 15 not taken.
 26 Linear ci(r[i][0]/(2*c[0][0]), r[i][1]/(2*c[0][1]));
385
1/2
✓ Branch 2 taken 26 times.
✗ Branch 3 not taken.
 26 SBasis cisi = shift(ci, i);
386
6/12
✓ Branch 3 taken 26 times.
✗ Branch 4 not taken.
✓ Branch 9 taken 26 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 26 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 26 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 26 times.
✗ Branch 22 not taken.
 26 r -= multiply(shift((c*2 + cisi), i), SBasis(ci));
387
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 r.truncate(k+1);
388
1/2
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
 26 c += cisi;
389
3/4
✓ Branch 1 taken 26 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 22 times.
 26 if(r.tailError(i) == 0) // if exact
390 4 break;
391
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 4 times.
 26 }
392
393 4 return c;
394 4 }
395
396 /** Compute the recpirocal of a
397 \param a sbasis functions
398 \returns sbasis 1/a
399
400 It is recommended to use the piecewise version unless you have good reason.
401 */
402 ✗ SBasis reciprocal(Linear const &a, int k) {
403 ✗ SBasis c;
404 ✗ assert(!a.isZero());
405 ✗ c.resize(k, Linear(0,0));
406 ✗ double r_s0 = (a.tri()*a.tri())/(-a[0]*a[1]);
407 ✗ double r_s0k = 1;
408 ✗ for(unsigned i = 0; i < (unsigned)k; i++) {
409 ✗ c[i] = Linear(r_s0k/a[0], r_s0k/a[1]);
410 ✗ r_s0k *= r_s0;
411 }
412 ✗ return c;
413 ✗ }
414
415 /** Compute a / b to k terms
416 \param a,b sbasis functions
417 \returns sbasis a/b
418
419 It is recommended to use the piecewise version unless you have good reason.
420 */
421 2 SBasis divide(SBasis const &a, SBasis const &b, int k) {
422
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 SBasis c;
423
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 assert(!a.isZero());
424
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 SBasis r = a; // remainder
425
426 2 k++;
427
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 r.resize(k, Linear(0,0));
428
1/2
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 c.resize(k, Linear(0,0));
429
430
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 times.
 9 for(unsigned i = 0; i < (unsigned)k; i++) {
431
2/4
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
 8 Linear ci(r[i][0]/b[0][0], r[i][1]/b[0][1]); //H0
432
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 c[i] += ci;
433
4/8
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 8 times.
✗ Branch 14 not taken.
 8 r -= shift(multiply(ci,b), i);
434
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 r.truncate(k+1);
435
3/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 7 times.
 8 if(r.tailError(i) == 0) // if exact
436 1 break;
437 }
438
439 4 return c;
440 2 }
441
442 /** Compute a composed with b
443 \param a,b sbasis functions
444 \returns sbasis a(b(t))
445
446 return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
447 */
448 290607 SBasis compose(SBasis const &a, SBasis const &b) {
449
3/6
✓ Branch 6 taken 290607 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 290607 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 290607 times.
✗ Branch 13 not taken.
 290607 SBasis s = multiply((SBasis(Linear(1,1))-b), b);
450
1/2
✓ Branch 1 taken 290607 times.
✗ Branch 2 not taken.
 290607 SBasis r;
451
452
2/2
✓ Branch 1 taken 503419 times.
✓ Branch 2 taken 290607 times.
 794026 for(int i = a.size()-1; i >= 0; i--) {
453
7/14
✓ Branch 7 taken 503419 times.
✗ Branch 8 not taken.
✓ Branch 15 taken 503419 times.
✗ Branch 16 not taken.
✓ Branch 24 taken 503419 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 503419 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 503419 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 503419 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 503419 times.
✗ Branch 37 not taken.
 503419 r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]);
454 }
455 581214 return r;
456 290607 }
457
458 /** Compute a composed with b to k terms
459 \param a,b sbasis functions
460 \returns sbasis a(b(t))
461
462 return a0 + s(a1 + s(a2 +... where s = (1-u)u; ak =(1 - u)a^0_k + ua^1_k
463 */
464 ✗ SBasis compose(SBasis const &a, SBasis const &b, unsigned k) {
465 ✗ SBasis s = multiply((SBasis(Linear(1,1))-b), b);
466 ✗ SBasis r;
467
468 ✗ for(int i = a.size()-1; i >= 0; i--) {
469 ✗ r = multiply_add(r, s, SBasis(Linear(a[i][0])) - b*a[i][0] + b*a[i][1]);
470 }
471 ✗ r.truncate(k);
472 ✗ return r;
473 ✗ }
474
475 95000 SBasis portion(const SBasis &t, double from, double to) {
476 95000 double fv = t.valueAt(from);
477 95000 double tv = t.valueAt(to);
478
2/4
✓ Branch 4 taken 95000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 95000 times.
✗ Branch 8 not taken.
 95000 SBasis ret = compose(t, Linear(from, to));
479
1/2
✓ Branch 1 taken 95000 times.
✗ Branch 2 not taken.
 95000 ret.at0() = fv;
480
1/2
✓ Branch 1 taken 95000 times.
✗ Branch 2 not taken.
 95000 ret.at1() = tv;
481 95000 return ret;
482 ✗ }
483
484 /*
485 Inversion algorithm. The notation is certainly very misleading. The
486 pseudocode should say:
487
488 c(v) := 0
489 r(u) := r_0(u) := u
490 for i:=0 to k do
491 c_i(v) := H_0(r_i(u)/(t_1)^i; u)
492 c(v) := c(v) + c_i(v)*t^i
493 r(u) := r(u) ? c_i(u)*(t(u))^i
494 endfor
495 */
496
497 //#define DEBUG_INVERSION 1
498
499 /** find the function a^-1 such that a^-1 composed with a to k terms is the identity function
500 \param a sbasis function
501 \returns sbasis a^-1 s.t. a^-1(a(t)) = 1
502
503 The function must have 'unit range'("a00 = 0 and a01 = 1") and be monotonic.
504 */
505 7 SBasis inverse(SBasis a, int k) {
506
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 7 times.
 7 assert(a.size() > 0);
507 7 double a0 = a[0][0];
508
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
 7 if(a0 != 0) {
509 4 a -= a0;
510 }
511 7 double a1 = a[0][1];
512
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 7 times.
 7 assert(a1 != 0); // not invertable.
513
514
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 5 times.
 7 if(a1 != 1) {
515 2 a /= a1;
516 }
517
1/2
✓ Branch 3 taken 7 times.
✗ Branch 4 not taken.
 7 SBasis c(k, Linear()); // c(v) := 0
518
6/6
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 6 times.
 7 if(a.size() >= 2 && k == 2) {
519
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 1 c[0] = Linear(0,1);
520
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
 1 Linear t1(1+a[1][0], 1-a[1][1]); // t_1
521
3/6
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
 1 c[1] = Linear(-a[1][0]/t1[0], -a[1][1]/t1[1]);
522
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 4 times.
 6 } else if(a.size() >= 2) { // non linear
523
1/2
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 SBasis r = Linear(0,1); // r(u) := r_0(u) := u
524
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 Linear t1(1./(1+a[1][0]), 1./(1-a[1][1])); // 1./t_1
525 2 Linear one(1,1);
526 2 Linear t1i = one; // t_1^0
527
2/4
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 SBasis one_minus_a = SBasis(one) - a;
528
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 SBasis t = multiply(one_minus_a, a); // t(u)
529
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 SBasis ti(one); // t(u)^0
530 #ifdef DEBUG_INVERSION
531 std::cout << "a=" << a << std::endl;
532 std::cout << "1-a=" << one_minus_a << std::endl;
533 std::cout << "t1=" << t1 << std::endl;
534 //assert(t1 == t[1]);
535 #endif
536
537 //c.resize(k+1, Linear(0,0));
538
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 times.
 13 for(unsigned i = 0; i < (unsigned)k; i++) { // for i:=0 to k do
539 #ifdef DEBUG_INVERSION
540 std::cout << "-------" << i << ": ---------" <<std::endl;
541 std::cout << "r=" << r << std::endl
542 << "c=" << c << std::endl
543 << "ti=" << ti << std::endl
544 << std::endl;
545 #endif
546
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
 12 if(r.size() <= i) // ensure enough space in the remainder, probably not needed
547 ✗ r.resize(i+1, Linear(0,0));
548
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
 12 Linear ci(r[i][0]*t1i[0], r[i][1]*t1i[1]); // c_i(v) := H_0(r_i(u)/(t_1)^i; u)
549 #ifdef DEBUG_INVERSION
550 std::cout << "t1i=" << t1i << std::endl;
551 std::cout << "ci=" << ci << std::endl;
552 #endif
553
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 12 times.
 36 for(int dim = 0; dim < 2; dim++) // t1^-i *= 1./t1
554 24 t1i[dim] *= t1[dim];
555
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 c[i] = ci; // c(v) := c(v) + c_i(v)*t^i
556 // change from v to u parameterisation
557
3/6
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 12 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 12 times.
✗ Branch 13 not taken.
 12 SBasis civ = one_minus_a*ci[0] + a*ci[1];
558 // r(u) := r(u) - c_i(u)*(t(u))^i
559 // We can truncate this to the number of final terms, as no following terms can
560 // contribute to the result.
561
2/4
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
 12 r -= multiply(civ,ti);
562
1/2
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
 12 r.truncate(k);
563
3/4
✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 11 times.
 12 if(r.tailError(i) == 0)
564 1 break; // yay!
565
2/4
✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 11 times.
✗ Branch 6 not taken.
 11 ti = multiply(ti,t);
566
2/2
✓ Branch 1 taken 11 times.
✓ Branch 2 taken 1 times.
 12 }
567 #ifdef DEBUG_INVERSION
568 std::cout << "##########################" << std::endl;
569 #endif
570 2 } else
571
2/4
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 c = Linear(0,1); // linear
572
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 c -= a0; // invert the offset
573
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 7 c /= a1; // invert the slope
574 7 return c;
575 ✗ }
576
577 /** Compute the sine of a to k terms
578 \param b linear function
579 \returns sbasis sin(a)
580
581 It is recommended to use the piecewise version unless you have good reason.
582 */
583 46 SBasis sin(Linear b, int k) {
584
1/2
✓ Branch 3 taken 46 times.
✗ Branch 4 not taken.
 46 SBasis s(k+2, Linear());
585
1/2
✓ Branch 5 taken 46 times.
✗ Branch 6 not taken.
 46 s[0] = Linear(std::sin(b[0]), std::sin(b[1]));
586
1/2
✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
 46 double tr = s[0].tri();
587 46 double t2 = b.tri();
588
1/2
✓ Branch 5 taken 46 times.
✗ Branch 6 not taken.
 46 s[1] = Linear(std::cos(b[0])*t2 - tr, -std::cos(b[1])*t2 + tr);
589
590 46 t2 *= t2;
591
2/2
✓ Branch 0 taken 184 times.
✓ Branch 1 taken 46 times.
 230 for(int i = 0; i < k; i++) {
592
1/2
✓ Branch 3 taken 184 times.
✗ Branch 4 not taken.
 368 Linear bo(4*(i+1)*s[i+1][0] - 2*s[i+1][1],
593
3/6
✓ Branch 1 taken 184 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 184 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 184 times.
✗ Branch 10 not taken.
 368 -2*s[i+1][0] + 4*(i+1)*s[i+1][1]);
594
1/2
✓ Branch 3 taken 184 times.
✗ Branch 4 not taken.
 184 bo -= s[i]*(t2/(i+1));
595
596
597
1/2
✓ Branch 3 taken 184 times.
✗ Branch 4 not taken.
 184 s[i+2] = bo/double(i+2);
598 }
599
600 46 return s;
601 ✗ }
602
603 /** Compute the cosine of a
604 \param b linear function
605 \returns sbasis cos(a)
606
607 It is recommended to use the piecewise version unless you have good reason.
608 */
609 23 SBasis cos(Linear bo, int k) {
610 23 return sin(Linear(bo[0] + M_PI/2,
611 23 bo[1] + M_PI/2),
612
1/2
✓ Branch 4 taken 23 times.
✗ Branch 5 not taken.
 46 k);
613 }
614
615 /** compute fog^-1.
616 \param f,g sbasis functions
617 \returns sbasis f(g^-1(t)).
618
619 ("zero" = double comparison threshold. *!*we might divide by "zero"*!*)
620 TODO: compute order according to tol?
621 TODO: requires g(0)=0 & g(1)=1 atm... adaptation to other cases should be obvious!
622 */
623 23400 SBasis compose_inverse(SBasis const &f, SBasis const &g, unsigned order, double zero){
624
1/2
✓ Branch 3 taken 23400 times.
✗ Branch 4 not taken.
 23400 SBasis result(order, Linear(0.)); //result
625
1/2
✓ Branch 2 taken 23400 times.
✗ Branch 3 not taken.
 23400 SBasis r=f; //remainder
626
4/8
✓ Branch 5 taken 23400 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 23400 times.
✗ Branch 9 not taken.
✓ Branch 15 taken 23400 times.
✗ Branch 16 not taken.
✓ Branch 19 taken 23400 times.
✗ Branch 20 not taken.
 23400 SBasis Pk=Linear(1)-g,Qk=g,sg=Pk*Qk;
627
1/2
✓ Branch 1 taken 23400 times.
✗ Branch 2 not taken.
 23400 Pk.truncate(order);
628
1/2
✓ Branch 1 taken 23400 times.
✗ Branch 2 not taken.
 23400 Qk.truncate(order);
629
1/2
✓ Branch 3 taken 23400 times.
✗ Branch 4 not taken.
 23400 Pk.resize(order,Linear(0.));
630
1/2
✓ Branch 3 taken 23400 times.
✗ Branch 4 not taken.
 23400 Qk.resize(order,Linear(0.));
631
1/2
✓ Branch 3 taken 23400 times.
✗ Branch 4 not taken.
 23400 r.resize(order,Linear(0.));
632
633 23400 int vs = valuation(sg,zero);
634
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23400 times.
 23400 if (vs == 0) { // to prevent infinite loop
635 ✗ return result;
636 }
637
638
2/2
✓ Branch 0 taken 46800 times.
✓ Branch 1 taken 23400 times.
 70200 for (unsigned k=0; k<order; k+=vs){
639
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double p10 = Pk.at(k)[0];// we have to solve the linear system:
640
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double p01 = Pk.at(k)[1];//
641
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double q10 = Qk.at(k)[0];// p10*a + q10*b = r10
642
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double q01 = Qk.at(k)[1];// &
643
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double r10 = r.at(k)[0];// p01*a + q01*b = r01
644
1/2
✓ Branch 1 taken 46800 times.
✗ Branch 2 not taken.
 46800 double r01 = r.at(k)[1];//
645 double a,b;
646 46800 double det = p10*q01-p01*q10;
647
648 //TODO: handle det~0!!
649
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 46800 times.
 46800 if (fabs(det)<zero){
650 ✗ a=b=0;
651 }else{
652 46800 a=( q01*r10-q10*r01)/det;
653 46800 b=(-p01*r10+p10*r01)/det;
654 }
655
1/2
✓ Branch 3 taken 46800 times.
✗ Branch 4 not taken.
 46800 result[k] = Linear(a,b);
656
5/10
✓ Branch 3 taken 46800 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 46800 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 46800 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 46800 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 46800 times.
✗ Branch 18 not taken.
 46800 r=r-Pk*a-Qk*b;
657
658
2/4
✓ Branch 2 taken 46800 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 46800 times.
✗ Branch 6 not taken.
 46800 Pk=Pk*sg;
659
2/4
✓ Branch 2 taken 46800 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 46800 times.
✗ Branch 6 not taken.
 46800 Qk=Qk*sg;
660
661
1/2
✓ Branch 3 taken 46800 times.
✗ Branch 4 not taken.
 46800 Pk.resize(order,Linear(0.)); // truncates if too high order, expands with zeros if too low
662
1/2
✓ Branch 3 taken 46800 times.
✗ Branch 4 not taken.
 46800 Qk.resize(order,Linear(0.));
663
1/2
✓ Branch 3 taken 46800 times.
✗ Branch 4 not taken.
 46800 r.resize(order,Linear(0.));
664
665 }
666 23400 result.normalize();
667 23400 return result;
668 23400 }
669
670 }
671
672 /*
673 Local Variables:
674 mode:c++
675 c-file-style:"stroustrup"
676 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
677 indent-tabs-mode:nil
678 fill-column:99
679 End:
680 */
681 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
682