[Rd] polyroot() (PR#751)
Peter Dalgaard BSA
p.dalgaard@biostat.ku.dk
16 Jul 2001 18:43:18 +0200
p.dalgaard@biostat.ku.dk writes:
> 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...
Oops, something went wrong in the copying from the original report.
The test case is
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))))
^
Fell out.
Anyways, the bug is fixed now. Here's the patch if anyone is in a
hurry:
$ cvs diff -u src/appl/cpoly.c
Index: src/appl/cpoly.c
===================================================================
RCS file: /home/rdevel/CVS-ARCHIVE/R/src/appl/cpoly.c,v
retrieving revision 1.23
diff -u -r1.23 cpoly.c
--- src/appl/cpoly.c 26 Apr 2001 15:28:35 -0000 1.23
+++ src/appl/cpoly.c 16 Jul 2001 16:33:59 -0000
@@ -646,7 +646,7 @@
e = hypot(qr[0], qi[0]) * m_re / (a_re + m_re);
for (i=0; i < n; i++)
- e *= (ms + hypot(qr[i], qi[i]));
+ e = e*ms + hypot(qr[i], qi[i]);
return e * (a_re + m_re) - mp * m_re;
}
--
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._