From: Rolf Turner <r.turner_at_auckland.ac.nz>

Date: Mon, 26 May 2008 11:48:24 +1200

*>> 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.
*

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

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> 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

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 $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.
*