1
/** root finding for sbasis functions.
2
* Copyright 2006 N Hurst
3
* Copyright 2007 JF Barraud
5
* It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
8
* multi-roots using bernstein method, one approach would be:
10
take median and find roots of that
11
whenever a segment lies entirely on one side of the median,
12
find the median of the half and recurse.
14
in essence we are implementing quicksort on a continuous function
16
* the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
18
a) it requires convertion to poly, which is numerically unstable
20
b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
22
c) it finds all roots, even complex ones. We don't want to accidently treat a nearly real root as a real root
24
From memory gsl poly roots was about 10 times faster than bernstein in the case where all the roots
25
are in [0,1] for polys of order 5. I spent some time working out whether eigenvalue root finding
26
could be done directly in sbasis space, but the maths was too hard for me. -- njh
28
jfbarraud: eigenvalue root finding could be done directly in sbasis space ?
30
njh: I don't know, I think it should. You would make a matrix whose characteristic polynomial was
31
correct, but do it by putting the sbasis terms in the right spots in the matrix. normal eigenvalue
32
root finding makes a matrix that is a diagonal + a row along the top. This matrix has the property
33
that its characteristic poly is just the poly whose coefficients are along the top row.
35
Now an sbasis function is a linear combination of the poly coeffs. So it seems to me that you
36
should be able to put the sbasis coeffs directly into a matrix in the right spots so that the
37
characteristic poly is the sbasis. You'll still have problems b) and c).
39
We might be able to lift an eigenvalue solver and include that directly into 2geom. Eigenvalues
40
also allow you to find intersections of multiple curves but require solving n*m x n*m matrices.
48
#include "sbasis-to-bezier.h"
55
Interval bounds_exact(SBasis const &a) {
56
Interval result = Interval(a.at0(), a.at1());
57
SBasis df = derivative(a);
58
vector<double>extrema = roots(df);
59
for (unsigned i=0; i<extrema.size(); i++){
60
result.extendTo(a(extrema[i]));
65
Interval bounds_fast(const SBasis &sb, int order) {
67
for(int j = sb.size()-1; j>=order; j--) {
73
if (v<0) t = ((b-a)/v+1)*0.5;
74
if (v>=0 || t<0 || t>1) {
75
res[0] = std::min(a,b);
77
res[0]=lerp(t, a+v*t, b);
81
if (v>0) t = ((b-a)/v+1)*0.5;
82
if (v<=0 || t<0 || t>1) {
83
res[1] = std::max(a,b);
85
res[1]=lerp(t, a+v*t, b);
88
if (order>0) res*=pow(.25,order);
92
Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
93
double t0=i.min(), t1=i.max(), lo=0., hi=0.;
94
for(int j = sb.size()-1; j>=order; j--) {
99
if (lo<0) t = ((b-a)/lo+1)*0.5;
100
if (lo>=0 || t<t0 || t>t1) {
101
lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
103
lo = lerp(t, a+lo*t, b);
106
if (hi>0) t = ((b-a)/hi+1)*0.5;
107
if (hi<=0 || t<t0 || t>t1) {
108
hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
110
hi = lerp(t, a+hi*t, b);
113
Interval res = Interval(lo,hi);
114
if (order>0) res*=pow(.25,order);
118
//-- multi_roots ------------------------------------
119
// goal: solve f(t)=c for several c at once.
120
/* algo: -compute f at both ends of the given segment [a,b].
121
-compute bounds m<df(t)<M for df on the segment.
122
let c and C be the levels below and above f(a):
123
going from f(a) down to c with slope m takes at least time (f(a)-c)/m
124
going from f(a) up to C with slope M takes at least time (C-f(a))/M
125
From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
126
Do the same for b: compute some b' such that there are no roots in (b',b].
127
-if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
128
unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
129
making things tricky and unpleasant...
131
//TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
134
static int upper_level(vector<double> const &levels,double x,double tol=0.){
135
return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
138
static void multi_roots_internal(SBasis const &f,
140
std::vector<double> const &levels,
141
std::vector<std::vector<double> > &roots,
151
idx=upper_level(levels,0,vtol);
152
if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
153
roots[idx].push_back(a);
154
roots[idx].push_back(b);
160
// int idxa=upper_level(levels,fa);
161
// int idxb=upper_level(levels,fb);
163
// if (fa==levels[idxa]){
169
// int idx_min=std::min(idxa,idxb);
170
// int idx_max=std::max(idxa,idxb);
171
// if (idx_max==levels.size()) idx_max-=1;
172
// for(int i=idx_min;i<=idx_max; i++){
173
// double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
174
// if(a<t&&t<b) roots[t]=i;
179
//TODO: use different tol for t and f ?
180
//TODO: unsigned idx ? (remove int casts when fixed)
181
int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
182
if (idx==(int)levels.size()) idx-=1;
183
double c=levels.at(idx);
184
if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
185
roots[idx].push_back((a+b)/2);
190
int idxa=upper_level(levels,fa,vtol);
191
int idxb=upper_level(levels,fb,vtol);
193
Interval bs = bounds_local(df,Interval(a,b));
195
//first times when a level (higher or lower) can be reached from a or b.
196
double ta_hi,tb_hi,ta_lo,tb_lo;
197
ta_hi=ta_lo=b+1;//default values => no root there.
198
tb_hi=tb_lo=a-1;//default values => no root there.
200
if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
202
roots[idxa].push_back(a);
205
if (bs.max()>0 && idxa<(int)levels.size())
206
ta_hi=a+(levels.at(idxa )-fa)/bs.max();
207
if (bs.min()<0 && idxa>0)
208
ta_lo=a+(levels.at(idxa-1)-fa)/bs.min();
210
if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
212
roots[idxb].push_back(b);
215
if (bs.min()<0 && idxb<(int)levels.size())
216
tb_hi=b+(levels.at(idxb )-fb)/bs.min();
217
if (bs.max()>0 && idxb>0)
218
tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
222
t0=std::min(ta_hi,ta_lo);
223
t1=std::max(tb_hi,tb_lo);
224
//hum, rounding errors frighten me! so I add this +tol...
225
if (t0>t1+htol) return;//no root here.
227
if (fabs(t1-t0)<htol){
228
multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
230
double t,t_left,t_right,ft,ft_left,ft_right;
231
t_left =t_right =t =(t0+t1)/2;
232
ft_left=ft_right=ft=f(t);
233
int idx=upper_level(levels,ft,vtol);
234
if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
235
roots[idx].push_back(t);
236
//we do not want to count it twice (from the left and from the right)
242
multi_roots_internal(f,df,levels,roots,htol,vtol,t0 ,f(t0) ,t_left,ft_left);
243
multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1 ,f(t1) );
247
std::vector<std::vector<double> > multi_roots(SBasis const &f,
248
std::vector<double> const &levels,
254
std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
256
SBasis df=derivative(f);
257
multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
261
//-------------------------------------
264
double Laguerre_internal(SBasis const & p,
273
//std::cout << "xk = " << xk << std::endl;
276
double err = fabs(b);
277
double abx = fabs(xk);
278
for(int j = p.size()-2; j >= 0; j--) {
282
err = fabs(b) + abx*err;
285
err *= 1e-7; // magic epsilon for convergence, should be computed from tol
290
//if(std::norm(px) < tol*tol)
293
double H = G*G - f / px;
295
//std::cout << "G = " << G << "H = " << H;
296
double radicand = (n - 1)*(n*H-G*G);
297
//assert(radicand.real() > 0);
300
//std::cout << "radicand = " << radicand << std::endl;
301
if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
302
a = - std::sqrt(radicand);
304
a = std::sqrt(radicand);
305
//std::cout << "a = " << a << std::endl;
307
//std::cout << "a = " << a << std::endl;
310
//std::cout << "xk = " << xk << std::endl;
315
void subdiv_sbasis(SBasis const & s,
316
std::vector<double> & roots,
317
double left, double right) {
318
Interval bs = bounds_fast(s);
319
if(bs.min() > 0 || bs.max() < 0)
320
return; // no roots here
321
if(s.tailError(1) < 1e-7) {
322
double t = s[0][0] / (s[0][0] - s[0][1]);
323
roots.push_back(left*(1-t) + t*right);
326
double middle = (left + right)/2;
327
subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
328
subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
331
// It is faster to use the bernstein root finder for small degree polynomials (<100?.
333
std::vector<double> roots(SBasis const & s) {
334
if(s.size() == 0) return std::vector<double>();
336
return sbasis_to_bezier(s).roots();
344
c-file-style:"stroustrup"
345
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
350
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
added
removed
src/2geom/sbasis-roots.cpp
