Wilkinson Polynomial Condition Number

> library(polynom)
> p <- polynomial(1)
> for (i in 1:20) {
+     p <- p * polynomial(c(-i, 1))
+ }
> wilk <- as.function(p)
> print(p)
2432902008176640000 - 8752948036761600000*x + 13803759753640699904*x^2 -
1.2870931245151e+19*x^3 + 8037811822645052416*x^4 - 3599979517947607040*x^5 +
1206647803780372992*x^6 - 311333643161390720*x^7 + 63030812099294896*x^8 -
10142299865511450*x^9 + 1307535010540395*x^10 - 135585182899530*x^11 +
11310276995381*x^12 - 756111184500*x^13 + 40171771630*x^14 - 1672280820*x^15 +
53327946*x^16 - 1256850*x^17 + 20615*x^18 - 210*x^19 + x^20

> p <- polynomial(1)
> delta <- 10^(-10)
> for (i in 1:20) {
+     p <- p * polynomial(c(-i, 1))
+ }
> coefs <- coef(p)
> coefs[20] <- coefs[20] + delta
> q <- polynomial(coefs)
> roots <- solve(q)
> num <- (as.real(roots) - 1:20)/(1:20)
> denom <- delta/210
> print(abs(round(num/denom)))
 [1]           0           2          64         688        6023       21113
 [7]      740827    13409515   112938531   604315082  2382813387  6615177797
[13] 17187871827 26228890219 46778830844 45372718091 25027432627 14144780190
[19]  3418301551   473250162
> print(which.max(abs(round(num/denom))))
[1] 15

> p <- polynomial(1)
> delta <- 10^(-10)
> for (i in 1:20) {
+     p <- p * polynomial(c(-i, 1))
+ }
> coefs <- coef(p)
> coefs[16] <- coefs[16] + delta
> q <- polynomial(coefs)
> roots <- solve(q)
> num <- (as.real(roots) - 1:20)/(1:20)
> denom <- delta/210
> print(abs(round(num/denom)))
 [1]          0          6        275       4502      37013     188876
 [7]     798259    4087109   23229309  108196501  371447471  948653199
[13] 1770324181 2524013797 2641062547 2022974705 1124210309  411491811
[19]   92387368    9377744
> print(which.max(abs(round(num/denom))))
[1] 15

> p <- polynomial(1)
> for (i in 1:4) {
+     p <- p * polynomial(c(-i, 1))
+ }
> p
24 - 50*x + 35*x^2 - 10*x^3 + x^4
> wilk <- as.function(p)
> print(wilk(1:4))
[1] 0 0 0 0

> p <- polynomial(1)
> for (i in 1:16) {
+     p <- p * polynomial(c(-i, 1))
+ }
> p
20922789888000 - 70734282393600*x + 102992244837120*x^2 - 87077748875904*x^3 +
48366009233424*x^4 - 18861567058880*x^5 + 5374523477960*x^6 - 1146901283528*x^7
+ 185953177553*x^8 - 23057159840*x^9 + 2185031420*x^10 - 156952432*x^11 +
8394022*x^12 - 323680*x^13 + 8500*x^14 - 136*x^15 + x^16
> wilk <- as.function(p)
> print(wilk(1:16))
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

> p <- polynomial(1)
> for (i in 1:19) {
+     p <- p * polynomial(c(-i, 1))
+ }
> p
-121645100408832000 + 431565146817638400*x - 668609730341153280*x^2 +
610116075740491776*x^3 - 371384787345228032*x^4 + 161429736530119008*x^5 -
52260903362512720*x^6 + 12953636989943900*x^7 - 2503858755467550*x^8 +
381922055502195*x^9 - 46280647751910*x^10 + 4465226757381*x^11 -
342252511900*x^12 + 20692933630*x^13 - 973941900*x^14 + 34916946*x^15 -
920550*x^16 + 16815*x^17 - 190*x^18 + x^19
> wilk <- as.function(p)
> print(wilk(1:19))
 [1]        0        0        0     8192    38400    82944   150528   393216
 [9]   839808  1280000  1874048  3317760  5494528  7375872  9734400 14680064
[17] 21381376 26873856 33362176