sbasis-roots.cpp revision 8001ba81cb851b38d86650a2fef5817facffb763
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh/** root finding for sbasis functions.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * Copyright 2006 N Hurst
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * Copyright 2007 JF Barraud
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh *
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * It is more efficient to find roots of f(t) = c_0, c_1, ... all at once, rather than iterating.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh *
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * Todo/think about:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * multi-roots using bernstein method, one approach would be:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh sort c
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh take median and find roots of that
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh whenever a segment lies entirely on one side of the median,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh find the median of the half and recurse.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh in essence we are implementing quicksort on a continuous function
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh * the gsl poly roots finder is faster than bernstein too, but we don't use it for 3 reasons:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a) it requires convertion to poly, which is numerically unstable
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh b) it requires gsl (which is currently not a dependency, and would bring in a whole slew of unrelated stuff)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh c) it finds all roots, even complex ones. We don't want to accidently treat a nearly real root as a real root
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshFrom memory gsl poly roots was about 10 times faster than bernstein in the case where all the roots
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshare in [0,1] for polys of order 5. I spent some time working out whether eigenvalue root finding
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshcould be done directly in sbasis space, but the maths was too hard for me. -- njh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshjfbarraud: eigenvalue root finding could be done directly in sbasis space ?
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshnjh: I don't know, I think it should. You would make a matrix whose characteristic polynomial was
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshcorrect, but do it by putting the sbasis terms in the right spots in the matrix. normal eigenvalue
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshroot finding makes a matrix that is a diagonal + a row along the top. This matrix has the property
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshthat its characteristic poly is just the poly whose coefficients are along the top row.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshNow an sbasis function is a linear combination of the poly coeffs. So it seems to me that you
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshshould be able to put the sbasis coeffs directly into a matrix in the right spots so that the
fba63a357654d8b3e84c60007e40aa698cd45d19mikloshcharacteristic poly is the sbasis. You'll still have problems b) and c).
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshWe might be able to lift an eigenvalue solver and include that directly into 2geom. Eigenvalues
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshalso allow you to find intersections of multiple curves but require solving n*m x n*m matrices.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh **/
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
dc4f69a188c203f2fdc65f22d0d57904a8c52dd7miklosh#include <cmath>
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh#include <map>
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh#include <2geom/sbasis.h>
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh#include <2geom/sbasis-to-bezier.h>
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh#include <2geom/solver.h>
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshusing namespace std;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshnamespace Geom{
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322mikloshInterval bounds_exact(SBasis const &a) {
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh Interval result = Interval(a.at0(), a.at1());
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh SBasis df = derivative(a);
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh vector<double>extrema = roots(df);
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh for (unsigned i=0; i<extrema.size(); i++){
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh result.extendTo(a(extrema[i]));
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh }
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh return result;
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh}
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emikloshInterval bounds_fast(const SBasis &sb, int order) {
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh Interval res;
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh for(int j = sb.size()-1; j>=order; j--) {
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double a=sb[j][0];
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double b=sb[j][1];
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double v, t = 0;
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh v = res[0];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (v<0) t = ((b-a)/v+1)*0.5;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (v>=0 || t<0 || t>1) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh res[0] = std::min(a,b);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }else{
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh res[0]=lerp(t, a+v*t, b);
68664e00e2372534b4df2fdc5f54f836bafece18miklosh }
1cda9431ef400135f5e1bd899a94b921bdad0eafmiklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh v = res[1];
68664e00e2372534b4df2fdc5f54f836bafece18miklosh if (v>0) t = ((b-a)/v+1)*0.5;
dc4f69a188c203f2fdc65f22d0d57904a8c52dd7miklosh if (v<=0 || t<0 || t>1) {
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh res[1] = std::max(a,b);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh }else{
a4d12a5147f3d1d6b568a326e39ef5dca384248dmiklosh res[1]=lerp(t, a+v*t, b);
1667116521643e2475184b048e0abb77a2aa9735miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (order>0) res*=pow(.25,order);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return res;
1cda9431ef400135f5e1bd899a94b921bdad0eafmiklosh}
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
68664e00e2372534b4df2fdc5f54f836bafece18mikloshInterval bounds_local(const SBasis &sb, const Interval &i, int order) {
dc4f69a188c203f2fdc65f22d0d57904a8c52dd7miklosh double t0=i.min(), t1=i.max(), lo=0., hi=0.;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh for(int j = sb.size()-1; j>=order; j--) {
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh double a=sb[j][0];
1667116521643e2475184b048e0abb77a2aa9735miklosh double b=sb[j][1];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t = 0;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (lo<0) t = ((b-a)/lo+1)*0.5;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (lo>=0 || t<t0 || t>t1) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh lo = std::min(a*(1-t0)+b*t0+lo*t0*(1-t0),a*(1-t1)+b*t1+lo*t1*(1-t1));
1667116521643e2475184b048e0abb77a2aa9735miklosh }else{
1667116521643e2475184b048e0abb77a2aa9735miklosh lo = lerp(t, a+lo*t, b);
1667116521643e2475184b048e0abb77a2aa9735miklosh }
1667116521643e2475184b048e0abb77a2aa9735miklosh
1667116521643e2475184b048e0abb77a2aa9735miklosh if (hi>0) t = ((b-a)/hi+1)*0.5;
1667116521643e2475184b048e0abb77a2aa9735miklosh if (hi<=0 || t<t0 || t>t1) {
1667116521643e2475184b048e0abb77a2aa9735miklosh hi = std::max(a*(1-t0)+b*t0+hi*t0*(1-t0),a*(1-t1)+b*t1+hi*t1*(1-t1));
1667116521643e2475184b048e0abb77a2aa9735miklosh }else{
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh hi = lerp(t, a+hi*t, b);
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh }
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh }
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh Interval res = Interval(lo,hi);
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh if (order>0) res*=pow(.25,order);
1667116521643e2475184b048e0abb77a2aa9735miklosh return res;
1667116521643e2475184b048e0abb77a2aa9735miklosh}
1667116521643e2475184b048e0abb77a2aa9735miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh//-- multi_roots ------------------------------------
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh// goal: solve f(t)=c for several c at once.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh/* algo: -compute f at both ends of the given segment [a,b].
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh -compute bounds m<df(t)<M for df on the segment.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh let c and C be the levels below and above f(a):
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh going from f(a) down to c with slope m takes at least time (f(a)-c)/m
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh going from f(a) up to C with slope M takes at least time (C-f(a))/M
68664e00e2372534b4df2fdc5f54f836bafece18miklosh From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
68664e00e2372534b4df2fdc5f54f836bafece18miklosh Do the same for b: compute some b' such that there are no roots in (b',b].
68664e00e2372534b4df2fdc5f54f836bafece18miklosh -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
68664e00e2372534b4df2fdc5f54f836bafece18miklosh unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
68664e00e2372534b4df2fdc5f54f836bafece18miklosh making things tricky and unpleasant...
68664e00e2372534b4df2fdc5f54f836bafece18miklosh*/
68664e00e2372534b4df2fdc5f54f836bafece18miklosh//TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
68664e00e2372534b4df2fdc5f54f836bafece18miklosh
68664e00e2372534b4df2fdc5f54f836bafece18miklosh
68664e00e2372534b4df2fdc5f54f836bafece18mikloshstatic int upper_level(vector<double> const &levels,double x,double tol=0.){
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh}
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emikloshstatic void multi_roots_internal(SBasis const &f,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh SBasis const &df,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh std::vector<double> const &levels,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh std::vector<std::vector<double> > &roots,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double htol,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double vtol,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double a,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double fa,
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh double b,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double fb){
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (f.size()==0){
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh int idx;
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh idx=upper_level(levels,0,vtol);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idx<(int)levels.size()&&fabs(levels.at(idx))<=vtol){
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh roots[idx].push_back(a);
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh roots[idx].push_back(b);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return;
1667116521643e2475184b048e0abb77a2aa9735miklosh }
1667116521643e2475184b048e0abb77a2aa9735miklosh////usefull?
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (f.size()==1){
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh// int idxa=upper_level(levels,fa);
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idxb=upper_level(levels,fb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (fa==fb){
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (fa==levels[idxa]){
1667116521643e2475184b048e0abb77a2aa9735miklosh// roots[a]=idxa;
1667116521643e2475184b048e0abb77a2aa9735miklosh// roots[b]=idxa;
1667116521643e2475184b048e0abb77a2aa9735miklosh// }
1667116521643e2475184b048e0abb77a2aa9735miklosh// return;
1667116521643e2475184b048e0abb77a2aa9735miklosh// }
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idx_min=std::min(idxa,idxb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// int idx_max=std::max(idxa,idxb);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if (idx_max==levels.size()) idx_max-=1;
1667116521643e2475184b048e0abb77a2aa9735miklosh// for(int i=idx_min;i<=idx_max; i++){
1667116521643e2475184b048e0abb77a2aa9735miklosh// double t=a+(b-a)*(levels[i]-fa)/(fb-fa);
1667116521643e2475184b048e0abb77a2aa9735miklosh// if(a<t&&t<b) roots[t]=i;
1667116521643e2475184b048e0abb77a2aa9735miklosh// }
1667116521643e2475184b048e0abb77a2aa9735miklosh// return;
1667116521643e2475184b048e0abb77a2aa9735miklosh// }
1667116521643e2475184b048e0abb77a2aa9735miklosh if ((b-a)<htol){
1667116521643e2475184b048e0abb77a2aa9735miklosh //TODO: use different tol for t and f ?
1667116521643e2475184b048e0abb77a2aa9735miklosh //TODO: unsigned idx ? (remove int casts when fixed)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh int idx=std::min(upper_level(levels,fa,vtol),upper_level(levels,fb,vtol));
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh if (idx==(int)levels.size()) idx-=1;
68664e00e2372534b4df2fdc5f54f836bafece18miklosh double c=levels.at(idx);
1cda9431ef400135f5e1bd899a94b921bdad0eafmiklosh if((fa-c)*(fb-c)<=0||fabs(fa-c)<vtol||fabs(fb-c)<vtol){
68664e00e2372534b4df2fdc5f54f836bafece18miklosh roots[idx].push_back((a+b)/2);
68664e00e2372534b4df2fdc5f54f836bafece18miklosh }
68664e00e2372534b4df2fdc5f54f836bafece18miklosh return;
68664e00e2372534b4df2fdc5f54f836bafece18miklosh }
68664e00e2372534b4df2fdc5f54f836bafece18miklosh
68664e00e2372534b4df2fdc5f54f836bafece18miklosh int idxa=upper_level(levels,fa,vtol);
68664e00e2372534b4df2fdc5f54f836bafece18miklosh int idxb=upper_level(levels,fb,vtol);
68664e00e2372534b4df2fdc5f54f836bafece18miklosh
68664e00e2372534b4df2fdc5f54f836bafece18miklosh Interval bs = bounds_local(df,Interval(a,b));
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh //first times when a level (higher or lower) can be reached from a or b.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double ta_hi,tb_hi,ta_lo,tb_lo;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh ta_hi=ta_lo=b+1;//default values => no root there.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh tb_hi=tb_lo=a-1;//default values => no root there.
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh
1667116521643e2475184b048e0abb77a2aa9735miklosh if (idxa<(int)levels.size() && fabs(fa-levels.at(idxa))<vtol){//a can be considered a root.
c53f16f52840e8c0f2be9c1cc3af633c0ba1552emiklosh //ta_hi=ta_lo=a;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh roots[idxa].push_back(a);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh ta_hi=ta_lo=a+htol;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh }else{
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (bs.max()>0 && idxa<(int)levels.size())
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh ta_hi=a+(levels.at(idxa )-fa)/bs.max();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (bs.min()<0 && idxa>0)
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh ta_lo=a+(levels.at(idxa-1)-fa)/bs.min();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idxb<(int)levels.size() && fabs(fb-levels.at(idxb))<vtol){//b can be considered a root.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //tb_hi=tb_lo=b;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh roots[idxb].push_back(b);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh tb_hi=tb_lo=b-htol;
fba63a357654d8b3e84c60007e40aa698cd45d19miklosh }else{
fba63a357654d8b3e84c60007e40aa698cd45d19miklosh if (bs.min()<0 && idxb<(int)levels.size())
fba63a357654d8b3e84c60007e40aa698cd45d19miklosh tb_hi=b+(levels.at(idxb )-fb)/bs.min();
fba63a357654d8b3e84c60007e40aa698cd45d19miklosh if (bs.max()>0 && idxb>0)
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t0,t1;
d27f5758e12c3107ee69e66702043931e0756f6bmiklosh t0=std::min(ta_hi,ta_lo);
d27f5758e12c3107ee69e66702043931e0756f6bmiklosh t1=std::max(tb_hi,tb_lo);
d27f5758e12c3107ee69e66702043931e0756f6bmiklosh //hum, rounding errors frighten me! so I add this +tol...
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (t0>t1+htol) return;//no root here.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (fabs(t1-t0)<htol){
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t0,f(t0),t1,f(t1));
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh }else{
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t,t_left,t_right,ft,ft_left,ft_right;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh t_left =t_right =t =(t0+t1)/2;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh ft_left=ft_right=ft=f(t);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh int idx=upper_level(levels,ft,vtol);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if (idx<(int)levels.size() && fabs(ft-levels.at(idx))<vtol){//t can be considered a root.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh roots[idx].push_back(t);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh //we do not want to count it twice (from the left and from the right)
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh t_left =t-htol/2;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh t_right=t+htol/2;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh ft_left =f(t_left);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh ft_right=f(t_right);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh }
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t0 ,f(t0) ,t_left,ft_left);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,t_right,ft_right,t1 ,f(t1) );
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshstd::vector<std::vector<double> > multi_roots(SBasis const &f,
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh std::vector<double> const &levels,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double htol,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double vtol,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double a,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double b){
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh SBasis df=derivative(f);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return(roots);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh//-------------------------------------
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh#if 0
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshdouble Laguerre_internal(SBasis const & p,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double x0,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double tol,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh bool & quad_root) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double a = 2*tol;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double xk = x0;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double n = p.size();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh quad_root = false;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh while(a > tol) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "xk = " << xk << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh Linear b = p.back();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh Linear d(0), f(0);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double err = fabs(b);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double abx = fabs(xk);
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh for(int j = p.size()-2; j >= 0; j--) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh f = xk*f + d;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh d = xk*d + b;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh b = xk*b + p[j];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh err = fabs(b) + abx*err;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh err *= 1e-7; // magic epsilon for convergence, should be computed from tol
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double px = b;
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh if(fabs(b) < err)
dc4f69a188c203f2fdc65f22d0d57904a8c52dd7miklosh return xk;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //if(std::norm(px) < tol*tol)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh // return xk;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double G = d / px;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double H = G*G - f / px;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "G = " << G << "H = " << H;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double radicand = (n - 1)*(n*H-G*G);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //assert(radicand.real() > 0);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(radicand < 0)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh quad_root = true;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "radicand = " << radicand << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a = - std::sqrt(radicand);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh else
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a = std::sqrt(radicand);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "a = " << a << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh a = n / (a + G);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "a = " << a << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh xk -= a;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh //std::cout << "xk = " << xk << std::endl;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return xk;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh#endif
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshvoid subdiv_sbasis(SBasis const & s,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh std::vector<double> & roots,
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double left, double right) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh Interval bs = bounds_fast(s);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(bs.min() > 0 || bs.max() < 0)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return; // no roots here
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(s.tailError(1) < 1e-7) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double t = s[0][0] / (s[0][0] - s[0][1]);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh roots.push_back(left*(1-t) + t*right);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double middle = (left + right)/2;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh subdiv_sbasis(compose(s, Linear(0, 0.5)), roots, left, middle);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh subdiv_sbasis(compose(s, Linear(0.5, 1.)), roots, middle, right);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh// It is faster to use the bernstein root finder for small degree polynomials (<100?.
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshstd::vector<double> roots1(SBasis const & s) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh std::vector<double> res;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double d = s[0][0] - s[0][1];
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(d != 0) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh double r = s[0][0] / d;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh if(0 <= r and r <= 1)
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh res.push_back(r);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return res;
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9mikloshstd::vector<double> roots(SBasis const & s) {
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh switch(s.size()) {
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh case 0:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return std::vector<double>();
3711b3e25395437ee0a09dbbb2a76d999c4ef322miklosh case 1:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return roots1(s);
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh default:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh return sbasis_to_bezier(s).roots();
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh }
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh}
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh};
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh/*
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh Local Variables:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh mode:c++
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh c-file-style:"stroustrup"
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh indent-tabs-mode:nil
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh fill-column:99
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh End:
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh*/
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
7a7fa095a483e8b652af9f00e5169f62c84f09b9miklosh