Thanks Peter. That seems to work for me ( I tried bunch of 4th order polynomials, and I got the correct solutions)

On Thu, Sep 2, 2010 at 11:11 AM, Peter Vanroose <peter_vanroose@yahoo.co.uk> wrote:
> Is there any way around it ? I could scale down points by 10,
> and again scale back the results. But I want to get a generic
> solution as I use this program in a library and I would have
> no way of telling automatically when to scale and when not to scale.

Yes, of course, I see the problem.

A "generic" way would be to compare the first and last nonzero coefficient, and scale by a factor which makes them approximately equal (i.e., scale by the n-the root of that factor), then scale back the result.

I was thinking of adding a method "solve_robust()" to the vnl_rnpoly_solve class which does something like this, but I'm not feeling comfortable enough to do this right away for the general m-dimensional case. But for your particular (generic) use, it should not be too difficult: something like:

double factor = coeffs[4]/coeffs[0];
// take 4th root:
factor = sqrt(factor); factor = sqrt(factor);
double mfactor = 1.0;
for (int i=0; i<=4; ++i) coeffs[i]/=mfactor, mfactor*=factor;
...
vnl_rnpoly_solve solver(l);
vcl_vector<vnl_vector<double>*> re = solver.real();
vcl_vector<vnl_vector<double>*> im = solver.imag();
...
(*(*rp)) *= factor; // globally scale the vnl_vector back

I've put a working implementation of exactly this in the last SVN version of
core/vnl/algo/tests/test_rnpoly_roots.cxx


--      Peter.