@@ -15,7 +15,7 @@ to identify roots of polynomials with suspected multiplicities over
15
15
16
16
Example:
17
17
18
- ```jldoctest
18
+ ```
19
19
julia> using Polynomials
20
20
21
21
julia> p = fromroots([sqrt(2), sqrt(2), sqrt(2), 1, 1])
@@ -30,7 +30,7 @@ julia> roots(p)
30
30
1.4142350577588885 + 3.72273772278647e-5im
31
31
32
32
julia> Polynomials.Multroot.multroot(p)
33
- (values = [0.9999999999999993, 1.4142135623730958], multiplicities = [2, 3], κ = 5.218455674370637 , ϵ = 2.8736226244218195e-16)
33
+ (values = [0.9999999999999993, 1.4142135623730958], multiplicities = [2, 3], κ = 5.218455674370636 , ϵ = 2.8736226244218195e-16)
34
34
```
35
35
36
36
The algorithm has two stages. First it uses `pejorative_manifold` to
@@ -60,7 +60,7 @@ multiplicity structure:
60
60
recursively, this allows the residual tolerance to scale up to match
61
61
increasing numeric errors.
62
62
63
- Returns a named tuple with
63
+ Returns a named tuple with
64
64
65
65
* `values`: the identified roots
66
66
* `multiplicities`: the corresponding multiplicities
@@ -78,15 +78,15 @@ is misidentified.
78
78
"""
79
79
function multroot (p:: Polynomials.StandardBasisPolynomial{T} ; verbose= false ,
80
80
kwargs... ) where {T}
81
-
81
+
82
82
z, l = pejorative_manifold (p; kwargs... )
83
83
z̃ = pejorative_root (p, z, l)
84
84
κ, ϵ = stats (p, z̃, l)
85
-
85
+
86
86
verbose && show_stats (κ, ϵ)
87
87
88
88
(values = z̃, multiplicities = l, κ = κ, ϵ = ϵ)
89
-
89
+
90
90
end
91
91
92
92
# The multiplicity structure, `l`, gives rise to a pejorative manifold `Πₗ = {Gₗ(z) | z∈ Cᵐ}`.
@@ -99,11 +99,11 @@ function pejorative_manifold(p::Polynomials.StandardBasisPolynomial{T};
99
99
u, v, w, θ′, κ = Polynomials. ngcd (p, derivative (p),
100
100
atol= ρ* norm (p), satol = θ* norm (p),
101
101
rtol = zT, srtol = zT)
102
-
102
+
103
103
zs = roots (v)
104
104
nrts = length (zs)
105
105
ls = ones (Int, nrts)
106
-
106
+
107
107
while ! Polynomials. isconstant (u)
108
108
109
109
normp = 1 + norm (u, 2 )
@@ -127,7 +127,7 @@ function pejorative_manifold(p::Polynomials.StandardBasisPolynomial{T};
127
127
zs, ls # estimate for roots, multiplicities
128
128
129
129
end
130
-
130
+
131
131
"""
132
132
pejorative_root(p, zs, ls; kwargs...)
133
133
@@ -158,7 +158,7 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int};
158
158
159
159
# storage
160
160
a = p[2 : end ]. / p[1 ] # a ~ (p[n-1], p[n-2], ..., p[0])/p[n]
161
- W = Diagonal ([min (1 , 1 / abs (aᵢ)) for aᵢ in a])
161
+ W = Diagonal ([min (1 , 1 / abs (aᵢ)) for aᵢ in a])
162
162
J = zeros (S, m, n)
163
163
G = zeros (S, 1 + m)
164
164
Δₖ = zeros (S, n)
@@ -167,37 +167,37 @@ function pejorative_root(p, zs::Vector{S}, ls::Vector{Int};
167
167
cvg = false
168
168
δₖ₀ = - Inf
169
169
for ctr in 1 : maxsteps
170
-
170
+
171
171
evalJ! (J, zₖs, ls)
172
172
evalG! (G, zₖs, ls)
173
-
173
+
174
174
Δₖ .= (W* J) \ (W* (view (G, 2 : 1 + m) .- a)) # weighted least squares
175
175
176
176
δₖ₁ = norm (Δₖ, 2 )
177
177
Δ = δₖ₀ - δₖ₁
178
178
179
- if ctr > 10 && δₖ₁ >= δₖ₀
179
+ if ctr > 10 && δₖ₁ >= δₖ₀
180
180
δₖ₀ < δₖ₁ && @warn " Increasing Δ, terminating search"
181
181
cvg = true
182
182
break
183
183
end
184
184
185
185
zₖs .- = Δₖ
186
186
187
- if δₖ₁^ 2 <= (δₖ₀ - δₖ₁) * τ
187
+ if δₖ₁^ 2 <= (δₖ₀ - δₖ₁) * τ
188
188
cvg = true
189
189
break
190
190
end
191
191
192
192
δₖ₀ = δₖ₁
193
193
end
194
194
verbose && show_stats (stats (p, zₖs, ls)... )
195
-
195
+
196
196
if cvg
197
197
return zₖs
198
198
else
199
199
@info ("""
200
- The multiplicity count may be in error: the initial guess for the roots failed
200
+ The multiplicity count may be in error: the initial guess for the roots failed
201
201
to converge to a pejorative root.
202
202
""" )
203
203
return (zₘs)
@@ -242,7 +242,7 @@ function evalJ!(J, zs::Vector{T}, ls::Vector) where {T}
242
242
for (k, zₖ) in enumerate (zs)
243
243
k == j && continue
244
244
for i in n- 1 : - 1 : 1
245
- J[1 + i,j] -= zₖ * J[i,j]
245
+ J[1 + i,j] -= zₖ * J[i,j]
246
246
end
247
247
end
248
248
end
0 commit comments