src/2geom/sbasis-roots.cpp

	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