Re: [R] Solving 100th order equation

From: Rolf Turner <>
Date: Mon, 26 May 2008 11:48:24 +1200

On 26/05/2008, at 10:38 AM, Hans W. Borchers wrote:

> Rolf Turner <r.turner <at>> writes:
>> [...]
>> However uniroot(p,c(1.995,2.005)) gives
>> $root
>> [1] 1.999993
>> $f.root
>> [1] -4.570875e+24

	Whoops!  I read that e+24 as e-24, so scrub all that I said.
	(You'd think, that having seen Bill Venables make a similar
	 error --- which he corrected in the follow-up posting to which
	 I was replying --- I would've been more careful.  Well, you'd
	 think wrong.  Actually it *was* e-24 when I posted; then the
	 Gremlins got in and changed everything. :-) )

>> $iter
>> [1] 4
>> $estim.prec
>> [1] 6.103516e-05
>> What a difference 7.214144e-06 makes! When you're dealing with
>> polynomials of degree 100.
> I don't think R is the right tool to solve this kind of questions
> that belong to
> the realm of Computer Algebra Systems. Yacas is much to weak for
> such high-order
> polybomials, but we can apply more powerful CASystems, for instance
> the free
> Maxima system. Applying
> (%i) nroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, minf, inf)
> immediately shows that there are only *two* real solutions, and then
> (%i) a: realroots(x^100 - 2*x^99 + 10*x^50 + 6*x - 4000, 1e-15);
> (%i) float(a)
> (%o) [x=-1.074126672042147,x=1.999999999999982]
> will provide real solutions with 15 decimals (does not change when
> more decimals
> are used). So the difference that counts is actually much smaller.
> The complex roots -- if needed -- will require special treatment.
> I believe it would be fair to point such questions to Computer
> Algebra mailing
> lists and not try to give the appearance they could be solved
> satisfyingly with
> R. The same as a complicated statistics question in a Matlab or
> Mathematica
> mailing list should be pointed to R-help.
	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

	[1] 2

	[1] 912

	[1] 5

	[1] 8.881784e-16

	Well, 912 is a lot closer to 0 than -4.570875e+24, or even -1.407375e 
	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  
	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''.   
	isn't R.  I don't know if Maxima was using arbitrary precision  
	in the calculations shown.

	Anyhow, I think the real point is that solving 100-th degree  
	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.


			Rolf Turner

Attention:\ This e-mail message is privileged and confid...{{dropped:9}} mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. Received on Sun 25 May 2008 - 23:52:30 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 26 May 2008 - 08:30:41 GMT.

Mailing list information is available at Please read the posting guide before posting to the list.

list of date sections of archive