#!/usr/bin/env julia

# Jacek Poziemski 272389

import Pkg
Pkg.add("Polynomials")
using Polynomials

coefficients::Vector{Float64} = [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]

polynomial::Polynomial = Polynomial(reverse(coefficients))

polynomialroots::Vector{Float64} = roots(polynomial)

wilkinson::Polynomial = fromroots(polynomialroots)

println("P(n):")

for k in 1:20
    z_k::Float64 = polynomialroots[k]
    Pz_k::Float64 = abs(polynomial(z_k))
    pz_k::Float64 = abs(wilkinson(z_k))
    error::Float64 = abs(z_k - k)
    println("k: $k; z_k: $z_k; P(z_k): $Pz_k; p(z_k): $pz_k; error: $error")
end 

println()

wilkinsoncoefficients::Vector{Float64} = [1, -210.0 - 2.0^(-23.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]

polynomialb::Polynomial = Polynomial(reverse(wilkinsoncoefficients))

polynomialbroots::Vector{ComplexF64} = roots(polynomialb)

wilkinsonb::Polynomial = fromroots(polynomialbroots)

println("p(n):")

for k in 1:20
    z_k::ComplexF64 = polynomialbroots[k]
    Pz_k::Float64 = abs(polynomialb(z_k))
    pz_k::Float64 = abs(wilkinsonb(z_k))
    error::Float64 = abs(z_k - k)
    println("k: $k; z_k: $z_k; P(z_k): $Pz_k; p(z_k): $pz_k; error: $error")
end