[Rd] BUG: polyroot() (PR#751)
Peter Dalgaard BSA
p.dalgaard@biostat.ku.dk
07 Dec 2000 00:58:27 +0100
mavip1@inet.polyu.edu.hk writes:
> I have found that the polyroot()
> function in R-1.1.1(both solaris
> and Win32 version) gives totally
> incorrect result. Here is the offending
> code:
>
> # Polyroot bug report:
> # from R-1.1.1
> > sort(abs(polyroot(c(1, -2,1,0,0,0,0,0,0,0,0,0,-2,5,-2,0,0,0,0,0,0,0,0,0,1,-2,1))))
> [1] 0.8758259 0.9486499 0.9731015 1.5419189 1.7466214 1.7535362 1.7589484
> [8] 2.0216317 2.4421509 2.5098488 2.6615572 2.7717724 2.8444048 2.9147685
> [15] 2.9302278 2.9757159 3.2037151 3.2052962 3.2498510 3.2815881 3.3250881
> [22] 4.6593700 4.7048556 5.2738258 5.7369920 5.9017094
> # from Splus 2000
> > sort(abs(polyroot(c(1, -2,1,0,0,0,0,0,0,0,0,0,-2,5,-2,0,0,0,0,0,0,0,0,0,1,-2,1))))
> [1] 0.8038812 0.8038812 0.8758259 0.8758259 0.9235722
> [6] 0.9235722 0.9439788 0.9439788 0.9536692 0.9536692
> [11] 0.9582403 0.9582403 0.9596028 1.0420978 1.0435796
> [16] 1.0435796 1.0485816 1.0485816 1.0593458 1.0593458
> [21] 1.0827524 1.0827524 1.1417794 1.1417794 1.2439649
> [26] 1.2439649
>
> >From mathematical theory I know that the Splus result
> is correct while the R result is wrong. I would have
> checked the source if I had the time.
There seems to be a rather bad numerical instability with the R code.
If we define
polyeval<-function(coef,x){v<-0;for (i in rev(coef)) v<-v*x+i; v}
(Horner's scheme), we get
> polyeval(1:5,polyroot(1:5))
[1] -2.220446e-16-2.664833e-16i -2.220446e-16-1.554475e-16i
[3] -2.220446e-16+1.436839e-16i -2.220446e-16-5.076235e-16i
which is kind of OK except that it looks odd that the real part is the
same for all 4 roots. Increasing the degree gives
> polyeval(10:1,polyroot(10:1))
[1] 9.230968e-08+5.257143e-08i 9.230977e-08+5.257158e-08i
[3] 9.230961e-08+5.257162e-08i 9.230956e-08+5.257151e-08i
[5] 9.230974e-08+5.257145e-08i 9.230975e-08+5.257161e-08i
[7] 9.230949e-08+5.257178e-08i 9.230974e-08+5.257161e-08i
[9] 9.230969e-08+5.257184e-08i
which is fairly OK, but it does look odd that we can't get closer to
zero than that since there are no multiple roots, and again it is
peculiar that the values are so similar. The next items in the
sequence are -ahem- interesting:
> polyeval(11:1,polyroot(11:1))
[1] 2.305607e-06+1.791577e-05i 2.309550e-06+1.791714e-05i
[3] 2.307622e-06+1.792101e-05i 2.304150e-06+1.791869e-05i
[5] 2.307087e-06+1.791537e-05i 2.309716e-06+1.791873e-05i
[7] 2.306030e-06+1.792098e-05i 2.304727e-06+1.792011e-05i
[9] 2.308569e-06+1.791588e-05i 2.308994e-06+1.792018e-05i
> polyeval(12:1,polyroot(12:1))
[1] 70.94521-90.6597i 70.94521-90.6597i 70.94521-90.6597i 70.94521-90.6597i
[5] 70.94521-90.6597i 70.94521-90.6597i 70.94521-90.6597i 70.94521-90.6597i
[9] 70.94521-90.6597i 70.94521-90.6597i 70.94521-90.6597i
> polyeval(13:1,polyroot(13:1))
[1] -9.076647e-05-0.0003440688i -1.253803e-04-0.0003691299i
[3] -8.894859e-05-0.0003941022i -7.310528e-05-0.0003638818i
[5] -1.144849e-04-0.0003471283i -1.221901e-04-0.0003814892i
[7] -7.857251e-05-0.0003869614i -7.284137e-05-0.0003760067i
[9] -1.029218e-04-0.0003424790i -1.136443e-04-0.0003910680i
[11] -1.016283e-04-0.0003956270i -1.225930e-04-0.0003567609i
> polyeval(14:1,polyroot(14:1))
[1] 11.28833+ 28.72917i 11.28833+ 28.72917i
[3] 5809.74030- 427.90119i -19407.81458+ 3462.91985i
[5] -2823.47547- 494.63498i -40909.65435-10789.57874i
[7] -8336.95783+ 5241.26935i 344.71455- 1520.94876i
[9] 5357.44530-24665.94094i 57313.79892-43530.37668i
[11] 11464.67537+29484.08287i 30010.79343+24624.25551i
[13] 25799.32156+53755.92101i
> polyeval(15:1,polyroot(15:1))
[1] 4.110828e+01+5.028631e+01i 4.123746e+01+5.033469e+01i
[3] 4.113423e+01+5.043357e+01i 4.106175e+01+5.031924e+01i
[5] 2.130747e+06+7.226414e+05i -1.088699e+08+1.096129e+09i
[7] 1.785940e+09+1.407187e+09i 9.750674e+10-5.059480e+10i
[9] -3.102572e+09-9.212065e+10i 7.213065e+09-1.023845e+10i
[11] 1.682049e+10+2.202050e+10i 1.889123e+10-1.046838e+11i
[13] -6.220741e+10-4.563469e+10i -1.477005e+10-3.393510e+10i
Root finding in high-degree polynomials is known to be nasty business,
but this looks more than ordinarily bad.
--
O__ ---- Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._