![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
index 52d3ef6a9dcc960526d5fba87a88b4e488f28bc8..95fd0cf3ba8e45c19da4b878537ad7f84ea5576e 100644 (file)
* multi-roots using bernstein method, one approach would be:
sort c
take median and find roots of that
- whenever a segment lies entirely on one side of the median,
+ whenever a segment lies entirely on one side of the median,
find the median of the half and recurse.
in essence we are implementing quicksort on a continuous function
#include <cmath>
#include <map>
-#include "sbasis.h"
-#include "sbasis-to-bezier.h"
-#include "solver.h"
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-to-bezier.h>
+#include <2geom/solver.h>
using namespace std;
namespace Geom{
-Interval bounds_exact(SBasis const &a) {
+/** Find the smallest interval that bounds a
+ \param a sbasis function
+ \returns inteval
+
+*/
+
+#ifdef USE_SBASIS_OF
+OptInterval bounds_exact(SBasisOf<double> const &a) {
+ Interval result = Interval(a.at0(), a.at1());
+ SBasisOf<double> df = derivative(a);
+ vector<double>extrema = roots(df);
+ for (unsigned i=0; i<extrema.size(); i++){
+ result.extendTo(a(extrema[i]));
+ }
+ return result;
+}
+#else
+OptInterval bounds_exact(SBasis const &a) {
Interval result = Interval(a.at0(), a.at1());
SBasis df = derivative(a);
vector<double>extrema = roots(df);
}
return result;
}
+#endif
+
+/** Find a small interval that bounds a
+ \param a sbasis function
+ \returns inteval
+
+*/
+// I have no idea how this works, some clever bounding argument by jfb.
+#ifdef USE_SBASIS_OF
+OptInterval bounds_fast(const SBasisOf<double> &sb, int order) {
+#else
+OptInterval bounds_fast(const SBasis &sb, int order) {
+#endif
+ Interval res(0,0); // an empty sbasis is 0.
-Interval bounds_fast(const SBasis &sb, int order) {
- Interval res;
for(int j = sb.size()-1; j>=order; j--) {
double a=sb[j][0];
double b=sb[j][1];
return res;
}
-Interval bounds_local(const SBasis &sb, const Interval &i, int order) {
- double t0=i.min(), t1=i.max(), lo=0., hi=0.;
+/** Find a small interval that bounds a(t) for t in i to order order
+ \param sb sbasis function
+ \param i domain interval
+ \param order number of terms
+ \return interval
+
+*/
+#ifdef USE_SBASIS_OF
+OptInterval bounds_local(const SBasisOf<double> &sb, const OptInterval &i, int order) {
+#else
+OptInterval bounds_local(const SBasis &sb, const OptInterval &i, int order) {
+#endif
+ double t0=i->min(), t1=i->max(), lo=0., hi=0.;
for(int j = sb.size()-1; j>=order; j--) {
double a=sb[j][0];
double b=sb[j][1];
going from f(a) up to C with slope M takes at least time (C-f(a))/M
From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
Do the same for b: compute some b' such that there are no roots in (b',b].
- -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
+ -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
making things tricky and unpleasant...
*/
return(upper_bound(levels.begin(),levels.end(),x-tol)-levels.begin());
}
+#ifdef USE_SBASIS_OF
static void multi_roots_internal(SBasis const &f,
SBasis const &df,
+#else
+static void multi_roots_internal(SBasis const &f,
+ SBasis const &df,
+#endif
std::vector<double> const &levels,
std::vector<std::vector<double> > &roots,
double htol,
double fa,
double b,
double fb){
-
+
if (f.size()==0){
int idx;
idx=upper_level(levels,0,vtol);
}
return;
}
-////usefull?
+////usefull?
// if (f.size()==1){
// int idxa=upper_level(levels,fa);
// int idxb=upper_level(levels,fb);
}
return;
}
-
+
int idxa=upper_level(levels,fa,vtol);
int idxb=upper_level(levels,fb,vtol);
- Interval bs = bounds_local(df,Interval(a,b));
+ Interval bs = *bounds_local(df,Interval(a,b));
//first times when a level (higher or lower) can be reached from a or b.
double ta_hi,tb_hi,ta_lo,tb_lo;
if (bs.max()>0 && idxb>0)
tb_lo=b+(levels.at(idxb-1)-fb)/bs.max();
}
-
+
double t0,t1;
- t0=std::min(ta_hi,ta_lo);
+ t0=std::min(ta_hi,ta_lo);
t1=std::max(tb_hi,tb_lo);
//hum, rounding errors frighten me! so I add this +tol...
if (t0>t1+htol) return;//no root here.
}
}
+/** Solve f(t)=c for several c at once.
+ \param f sbasis function
+ \param levels vector of 'y' values
+ \param htol, vtol
+ \param a, b left and right bounds
+ \returns a vector of vectors, one for each y giving roots
+
+Effectively computes:
+results = roots(f(y_i)) for all y_i
+
+* algo: -compute f at both ends of the given segment [a,b].
+ -compute bounds m<df(t)<M for df on the segment.
+ let c and C be the levels below and above f(a):
+ going from f(a) down to c with slope m takes at least time (f(a)-c)/m
+ going from f(a) up to C with slope M takes at least time (C-f(a))/M
+ From this we conclude there are no roots before a'=a+min((f(a)-c)/m,(C-f(a))/M).
+ Do the same for b: compute some b' such that there are no roots in (b',b].
+ -if [a',b'] is not empty, repeat the process with [a',(a'+b')/2] and [(a'+b')/2,b'].
+ unfortunately, extra care is needed about rounding errors, and also to avoid the repetition of roots,
+ making things tricky and unpleasant...
+
+TODO: Make sure the code is "rounding-errors proof" and take care about repetition of roots!
+*/
std::vector<std::vector<double> > multi_roots(SBasis const &f,
std::vector<double> const &levels,
double htol,
std::vector<std::vector<double> > roots(levels.size(), std::vector<double>());
SBasis df=derivative(f);
- multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
+ multi_roots_internal(f,df,levels,roots,htol,vtol,a,f(a),b,f(b));
return(roots);
}
//-------------------------------------
-#if 0
-double Laguerre_internal(SBasis const & p,
- double x0,
- double tol,
- bool & quad_root) {
- double a = 2*tol;
- double xk = x0;
- double n = p.size();
- quad_root = false;
- while(a > tol) {
- //std::cout << "xk = " << xk << std::endl;
- Linear b = p.back();
- Linear d(0), f(0);
- double err = fabs(b);
- double abx = fabs(xk);
- for(int j = p.size()-2; j >= 0; j--) {
- f = xk*f + d;
- d = xk*d + b;
- b = xk*b + p[j];
- err = fabs(b) + abx*err;
- }
-
- err *= 1e-7; // magic epsilon for convergence, should be computed from tol
-
- double px = b;
- if(fabs(b) < err)
- return xk;
- //if(std::norm(px) < tol*tol)
- // return xk;
- double G = d / px;
- double H = G*G - f / px;
-
- //std::cout << "G = " << G << "H = " << H;
- double radicand = (n - 1)*(n*H-G*G);
- //assert(radicand.real() > 0);
- if(radicand < 0)
- quad_root = true;
- //std::cout << "radicand = " << radicand << std::endl;
- if(G.real() < 0) // here we try to maximise the denominator avoiding cancellation
- a = - std::sqrt(radicand);
- else
- a = std::sqrt(radicand);
- //std::cout << "a = " << a << std::endl;
- a = n / (a + G);
- //std::cout << "a = " << a << std::endl;
- xk -= a;
- }
- //std::cout << "xk = " << xk << std::endl;
- return xk;
-}
-#endif
void subdiv_sbasis(SBasis const & s,
- std::vector<double> & roots,
+ std::vector<double> & roots,
double left, double right) {
- Interval bs = bounds_fast(s);
- if(bs.min() > 0 || bs.max() < 0)
+ OptInterval bs = bounds_fast(s);
+ if(!bs || bs->min() > 0 || bs->max() < 0)
return; // no roots here
if(s.tailError(1) < 1e-7) {
double t = s[0][0] / (s[0][0] - s[0][1]);
// It is faster to use the bernstein root finder for small degree polynomials (<100?.
+std::vector<double> roots1(SBasis const & s) {
+ std::vector<double> res;
+ double d = s[0][0] - s[0][1];
+ if(d != 0) {
+ double r = s[0][0] / d;
+ if(0 <= r && r <= 1)
+ res.push_back(r);
+ }
+ return res;
+}
+
+/** Find all t s.t s(t) = 0
+ \param a sbasis function
+ \returns vector of zeros (roots)
+
+*/
std::vector<double> roots(SBasis const & s) {
- if(s.size() == 0) return std::vector<double>();
-
- return sbasis_to_bezier(s).roots();
+ switch(s.size()) {
+ case 0:
+ return std::vector<double>();
+ case 1:
+ return roots1(s);
+ default:
+ {
+ Bezier bz;
+ sbasis_to_bezier(bz, s);
+ return bz.roots();
+ }
+ }
}
};
fill-column:99
End:
*/
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :