From p.dalgaard@biostat.ku.dk  Mon Jul 16 17:43:18 2001
From: p.dalgaard@biostat.ku.dk (Peter Dalgaard BSA)
Date: 16 Jul 2001 18:43:18 +0200
Subject: [Rd] polyroot() (PR#751)
In-Reply-To: p.dalgaard@biostat.ku.dk's message of "Mon, 16 Jul 2001 11:35:05 +0200 (MET DST)"
References: <200107160935.LAA10917@pubhealth.ku.dk>
Message-ID: <x2k8189655.fsf@blueberry.kubism.ku.dk>

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


[Rd] polyroot() (PR#751)
Peter Dalgaard BSA