1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
using OffsetArrays
using Printf
using Polynomials
BF = BigFloat
setprecision(128)
function emptyArray(type, size)
return OffsetVector{type}(undef, 0:(size - 1))
end
function bairstow(n, p, u, v, iterations)
b = emptyArray(BF, n + 1)
c = emptyArray(BF, n + 1)
b[n] = p[n]
c[n] = 0
c[n - 1] = p[n]
for j in 1:iterations
b[n - 1] = p[n - 1] + u * b[n]
for k in range(n - 2, 0, step=-1)
b[k] = p[k] + u * b[k + 1] + v * b[k + 2]
c[k] = b[k + 1] + u * c[k + 1] + v * c[k + 2]
end
J = c[0] * c[2] - c[1] * c[1]
u += (c[1] * b[1] - c[2] * b[0]) / J
v += (c[1] * b[0] - c[0] * b[1]) / J
end
return (u, v, b)
end
function evalPolynomial(p, x)
res = BF(0)
for i in range(size(p)[1] - 1, 0, step=-1)
res = p[i] + res * x
end
return res
end
function solveQuadratic(p)
if p[1] != 0
delta = p[2]^2 - 4p[1]p[3]
delta += (delta < 0) ? 0im : 0
if p[2] >= 0
b = -p[2] - delta^0.5
return (b / (2p[1]), 2p[3] / b) # <=> P[3]/P[1] = b / (2P[1]) * X <=> X = 2P[3] / b
else
b = -p[2] + delta^0.5
return [2p[3] / b, b / (2p[1])]
end
elseif p[2] != 0
return [-p[3] / p[2]]
else
return []
end
end
function findRoots(p, u, v, iterations)
u, v, b = bairstow(4, p, u, v, iterations)
z1, z2 = solveQuadratic([1, -u, -v])
z3, z4 = solveQuadratic([b[4], b[3], b[2]])
return [z1, z2, z3, z4]
end
w = emptyArray(BF, 5)
w[0] = 1
w[1] = 2
w[2] = 3
w[3] = 4
w[4] = 5
u0 = BF(0.1)
v0 = BF(0.1)
roots_of_w = findRoots(w, u0, v0, 20)
for z in roots_of_w
# println(evalPolynomial(w, z))
res = evalPolynomial(w, z)
@printf("Evaluation at root %.16f + %.16fi: %.16f + %.16fi\n", real(z), imag(z), real(res), imag(res))
end
roots_of_w2 = roots(Polynomial([1, 2, 3, 4, 5]))
function relError(x, y)
return abs((x - y) / x)
end
for (z1, z2) in zip(roots_of_w, roots_of_w2)
error = relError(z1, z2)
@printf("Błąd względny między %.16f + %.16fi oraz %.16f + %.16fi: %.16f\n",
real(z1), imag(z1), real(z2), imag(z2), error)
end
|
