[R] Solving 100th order equation

[Hans W. Borchers]

Rolf Turner <r.turner <at> auckland.ac.nz> writes:

> 
> 

Dear Rolf Turner,

that is exactly what I say: Use the right tools !

And yes, Maxima is using multi-precision arithmetic. You can validate 
the result applying Mathematica or even the open source 'mpmath' package 
for Python and you will get (both tools agree):

    a = 1.999999999999982236431...

showing how far away the results of a C-like 'float' arithmetic can get 
in your examples below.

By the way, you can compute p(2) by hand because p(x) = x^99(x-2) +-... 
and you get p(2) = 10*2^50 - 3988 .

And yes, sometimes it will be necessary to solve high-order polynomials 
with great precision, probably not in statistics, but even in physics or
astronomy. And as has been mentioned in this thread before, finding roots 
of polynomials corresponds to determining eigenvalues of nxn-matrices, so
for ill-conditioned matrices similar problems will arise.

Very best,  Hans Werner


> 
> 	Be that as it were, it doesn't look like Maxima is doing so wonderfully
> 	either. Consider:
> 
> 	 > library(PolynomF)
> 	 > x <- polynom()
> 	 > p <- x^100-2*x^99+10*x^50+6*x-4000
> 	 > a <- 1.999999999999982
> 	 > p(a)
> 	[1] -1.407375e+14
> 
> 	And notice the + sign in the exponent. 
> 
> 	I tried doing
> 	> uniroot(p,c(1.9995,2.0005),tol=.Machine$double.eps)
> 	and got
> 	$root
> 	[1] 2
> 	$f.root
> 	[1] 912
> 
> Well, 912 is a lot closer to 0 than -4.570875e+24, or even -1.407375e+14
> But it still ain't 0.  Also that ``2'' of course isn't really 2.
> 	
> I set r <- uniroot(p,c(1.9995,2.0005),tol=.Machine$double.eps)$root
> I get p(r) = 912.
> 
> But print(r,digits=22) gives  1.999999999999982.
> And print(a,digits=22) gives the identical result.
> Although p(r) and p(a) are wildly different.
> I guess it's possible that p(a) --- where ``a'' is as it appears,
> written out to 14 decimal places --- really is quite close to zero,
> and that -1.407375e+14 is due to round-off error in the calculation
> process.  But it's also possible that p(a) really is pretty close to
> -1.407375e+14, i.e. Maxima isn't really finding the correct root either.
> 
> 	Or something in between.
> 
> I guess that to get a really ``correct'' answer, and check it properly,
> you need a system that will do ``arbitrary precision arithmetic''.   
> Which isn't R.  I don't know if Maxima was using arbitrary precision  
> arithmetic in the calculations shown.
> 
> Anyhow, I think the real point is that solving 100-th degree polynomials
> is a pretty silly thing to do, whatever package or system or language you use.
> That is, R may be the wrong tool for the job, but it's the wrong job.
> 
> And if you *must* do it, then you should do some careful checking  
> and investigation of the results --- again irrespective of whatever 
> package or system or language you use.
> 
> 		cheers,
> 			Rolf Turner