# Re: [R] Solving 100th order equation

From: Rolf Turner <r.turner_at_auckland.ac.nz>
Date: Mon, 26 May 2008 11:48:24 +1200

> Rolf Turner <r.turner <at> auckland.ac.nz> 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

I tried doing

> uniroot(p,c(1.9995,2.0005),tol=.Machine\$double.eps)

and got

\$root
[1] 2

\$f.root
[1] 912

\$iter
[1] 5

\$estim.prec
[1] 8.881784e-16

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

######################################################################
```
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html 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 https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.