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