using Polynomials

coefficients_array = [1, -210.0, 20615.0,-1256850.0,
      53327946.0,-1672280820.0, 40171771630.0, -756111184500.0,
      11310276995381.0, -135585182899530.0,
      1307535010540395.0,     -10142299865511450.0,
      63030812099294896.0,     -311333643161390640.0,
      1206647803780373360.0,     -3599979517947607200.0,
      8037811822645051776.0,      -12870931245150988800.0,
      13803759753640704000.0,      -8752948036761600000.0,
      2432902008176640000.0]

roots_array = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

new_coefficients_array = [1, -210.0 - 2.0^-23, 20615.0,-1256850.0,
      53327946.0,-1672280820.0, 40171771630.0, -756111184500.0,
      11310276995381.0, -135585182899530.0,
      1307535010540395.0,     -10142299865511450.0,
      63030812099294896.0,     -311333643161390640.0,
      1206647803780373360.0,     -3599979517947607200.0,
      8037811822645051776.0,      -12870931245150988800.0,
      13803759753640704000.0,      -8752948036761600000.0,
      2432902008176640000.0]


P = Poly(flipdim(coefficients_array, 1)) #flipdim because Poly takes coefficient from the lowest order, e. g. a_0, a_1, ...
p = poly(roots_array)

P_roots = roots(P)
p_roots = roots(p)

println("Roots for P(x): ", P_roots)
println("Roots for p(x): ", p_roots)

println("Checking |P(z_k)|: ")
for i = 1 : size(P_roots)[1]
      println("P(",i,") = ", P(P_roots[i]))
end