[Rd] polyroot() (PR#751)

[p.dalgaard@biostat.ku.dk]

In a bug report from Nov.28 2000, Li Dongfeng 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] 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] 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.
----- (end Dongfeng)

This was analyzed somewhat, and it was found that polyroot() does
work, sort of, for lower degree polynomials, but becomes numerically
unstable when the degree is large. The exact result is
compiler-dependent, but the problem is still with us:

> 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] 0.8025898 0.8828790 0.9270844 0.9412183 0.9413746 0.9498326 0.9506305
 [8] 0.9999350 1.0687862 1.0759356 1.2420437 1.2883753 1.3198206 1.4519195
[15] 1.4768297 1.6330082 1.6838551 1.6970280 1.7250806 1.7471799 1.8787228
[22] 1.8989308 1.9040403 1.9357768 1.9907273

It has been nagging us from the bug repository ever since.

Last night I got the bright idea of feeding that polynomial to the
TOMS-419 routine off of Netlib, and got results that agree with the
ones stated for S-PLUS. Now this *is* the algorithm that the one in R
is supposed to be a C transcription of, so there would seem to be a
transcription error somewhere...

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