Skip to content
Snippets Groups Projects
Commit 5974fd1c authored by Krueger Jasmin's avatar Krueger Jasmin
Browse files

debugging Hensel validation

parent f8adbe2d
Branches
No related tags found
No related merge requests found
......@@ -721,25 +721,32 @@ function validate_Hensel_a_posteriori(l::Int, f::PolyElem, g::PolyElem, h::PolyE
t = convert_float_to_balls(t)
#compute δ
(q,R) = fastdivrem(s*(igstar*ihstar-f),ihstar, l, validation = true)
δg = t*(igstar*ihstar-f) + q*g
δg = t*(igstar*ihstar-f) + q*igstar
δg = mag_coefficientwise(δg)
#BivPolMod!(δg,l)
δg = trunc(δg, dg+1,l)
#@show δg
δh = mag_coefficientwise(R)
#BivPolMod!(δh,l)
δh = trunc(δh, dh,l)
#@show δh
δgmin = nonzeroMin(δg)
ηg = ϵ*deepcopy(δg)
replace_zeros!(ηg,ϵ*δgmin, degree(g), l)
replace_zeros!(ηg,ϵ*δgmin, dg+1, l)
δhmin = nonzeroMin(δh)
ηh = ϵ*deepcopy(δh)
replace_zeros!(ηh, ϵ*δhmin, degree(h), l)
replace_zeros!(ηh, ϵ*δhmin, dh, l)
δh = δh + ηh
δg = δg + ηg
@show δg
@show δh
@show degree(δg)
@show degree(δh)
invrevh = newtoninversion(ihstar,l)
invrevh = newtoninversion(rev(ihstar,degree(ihstar)+1),l,validation = true)
u = s*igstar + t*ihstar -1
#validation
......@@ -749,12 +756,18 @@ function validate_Hensel_a_posteriori(l::Int, f::PolyElem, g::PolyElem, h::PolyE
(Λg, Λh) = Λ(rg,rh,s,t,u, igstar, ihstar, invrevh)
Λg = trunc(Λg, dg+1,l)
Λh = trunc(Λh, dh, l)
rrg = mag_coefficientwise(Λg*rg + δg)
rrg = trunc(rrg, dg+1, l)
@show degree(rrg)
rrh = mag_coefficientwise(Λh*rh + δh)
rrh = trunc(rrh, dh, l)
@show degree(rrh)
@show (isless_coefficientwise(rrg-rg,ηg))
@show (isless_coefficientwise(rrh-rh,ηh))
if (isless_coefficientwise(rrg-rg,ηg))&&(isless_coefficientwise(rrh-rh,ηh))
gstar = convert_polynomial_to_intervals(gstar,convert_balls_to_Float(rg,rounding = RoundUp))
hstar = convert_polynomial_to_intervals(hstar,convert_balls_to_Float(rh,rounding = RoundUp))
......@@ -787,13 +800,13 @@ function Λ(rg,rh, s,t, ϵ, g, h, invrevh)
reva = rev(a,degree(a)+1)
divabyh = reva*invrevh
Λg = 2*t*irg*irh + ϵ*irg + g*rev(divabyh, degree(divabyh)+1)
Λg = mag_coefficientwise(Λg)
Λg = mag_coefficientwise(-Λg)
#Λh
a2 = 2*s*irh*irg +ϵ*irh
reva2 = rev(a2, degree(a2)+1)
qstar = reva2*invrevh
q = rev(qstar, degree(qstar)+1)
Λh = a2 - h*q
Λh = mag_coefficientwise(Λh)
Λh = mag_coefficientwise(-Λh)
(Λg, Λh)
end
\ No newline at end of file
......@@ -18,7 +18,7 @@ function abs_coefficientwise!(p::PolyElem{T}) where T
end
"""
isless(p::PolyElem{T},q::PolyElem{T}) where T <: Union{FieldElem, AbstractFloat}
isless_coefficientwise(p::PolyElem{S},q::PolyElem{T}) where {S <: Union{FieldElem, AbstractFloat}, T <: Union{FieldElem, AbstractFloat}}
Compare two polynomials and return true if all p_i < q_i ∀ i.
"""
......@@ -69,6 +69,8 @@ function isless_coefficientwise(p::PolyElem{T},q::PolyElem{T}) where T <: arb_po
#print("[",Float64(coeff(coeff(p,i-1),j-1))-Float64(radius(coeff(coeff(p,i-1),j-1))),", ", Float64(coeff(coeff(p,i-1),j-1))+Float64(radius(coeff(coeff(p,i-1),j-1))), "] < [", Float64(coeff(coeff(q,i-1),j-1)) -Float64(radius(coeff(coeff(q,i-1),j-1))),", ", Float64(coeff(coeff(q,i-1),j-1)) +Float64(radius(coeff(coeff(q,i-1),j-1))), "] ?\n")
if ((coeff(coeff(p,i-1),j-1) < coeff(coeff(q,i-1),j-1)) == false)
@show (i-1)
@show (j-1)
@show (coeff(coeff(p,i-1),j-1))
@show coeff(coeff(q,i-1),j-1)
return false
......@@ -430,7 +432,7 @@ function replace_zeros!(p::PolyElem{T}, d, dp = -1, k = -1) where T
else
deg = k
end
@show deg
@show (deg-1)
for i in 0:deg-1
if T <: PolyElem
set_coefficient!(p,i,replace_zeros!(coeff(p,i),d,dp,k))
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment