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

debugging Hensel_validate

parent 0df916ac
No related branches found
No related tags found
No related merge requests found
......@@ -94,6 +94,9 @@ fr = convert_polynomial(fq,[RXY,RX],RDF)
@time (gistar, histar) = henseltruncate(128,fi,gi,hi)
@time (gvstar, hvstar) = henseltruncate(128,fi,gi,hi,validation = true)#@time (grstar, hrstar) = henseltruncate(128,fr,gr,hr)
@time (grstar, hrstar) = henseltruncate(128,fr,gr,hr)
@time (grv,hrv) = hensel_validate(128,fi,gi,hi);
degree(grv)
degree(hrv)
norm1error(grstar,gqstar)
norm1error(hrstar,hqstar)
......
......@@ -142,21 +142,21 @@ function fastdivrem(a::PolyElem, b::PolyElem)
da<db && return (zero(a),a)
m = da-db
invrevb = newtoninversion(rev(b, length(b)), m+1)
@show invrevb
#@show invrevb
reva = deepcopy(rev(a, length(a))) #need to make a copy otherwise the following
#truncation changes the a as well
#reverse seems to work directly on the coefficients
reva = truncate(reva, m+1)
qstar = reva*invrevb
qstar = truncate(qstar, m+1)
@show qstar
#@show qstar
q = rev(qstar, length(qstar))
r = a -b*q
@show a
@show b
@show r
#@show a
#@show b
#@show r
r = truncate(r, db)
@show r
#@show r
q, r, invrevb
end
......@@ -251,7 +251,7 @@ function hensel_step_truncate(deg::Int, f::PolyElem, g::PolyElem, h::PolyElem, s
gstar = g+t*e+q*g
BivPolMod!(gstar, 2deg)
gstar = truncate(gstar, n+1)
@show r
#@show r
hstar = h+r
BivPolMod!(hstar, 2deg)
......@@ -259,7 +259,7 @@ function hensel_step_truncate(deg::Int, f::PolyElem, g::PolyElem, h::PolyElem, s
BivPolMod!(b, 2deg)
sb = s*b
BivPolMod!(sb, 2deg)
@show hstar
#@show hstar
(c,d) = fastdivrem(sb, hstar, 2deg, validation = validation)
sstar = s-d
BivPolMod!(sstar, 2deg)
......@@ -288,7 +288,7 @@ function henseltruncate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem, s::PolyEl
end
BivPolMod!(g,l)
BivPolMod!(h,l)
(g,h)
(g,h,s,t)
end
function henseltruncate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; validation::Bool = false)
......@@ -299,10 +299,15 @@ function henseltruncate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; validatio
end
BivPolMod!(g,l)
BivPolMod!(h,l)
(g,h)
(g,h,s,t)
end
function hensel_validate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; ϵ=1e-10, imax = 20)
dg = degree(g)
dh = degree(h)
df = degree(f)
@show dg
@show dh
#convert input polynomials to real polynomials in Float64
gr = convert_balls_to_Float(g)
hr = convert_balls_to_Float(h)
......@@ -320,36 +325,28 @@ function hensel_validate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; ϵ=1e-10
(q,r) = fastdivrem(s*e,h,l,validation = true)
δ_g = t*e +q*g
BivPolMod!(δ_g,l)
δ_g = truncate(δ_g,dg+1)
abs_coefficientwise!(δ_g)
δ_g = mag_coefficientwise(δ_g)
@show δ_g
δ_h = r
BivPolMod!(δ_h,l)
δ_h = truncate(δ_h,dh+1)
abs_coefficientwise!(δ_h)
δ_h = mag_coefficientwise(δ_h)
@show δ_h
δmin_g = nonzeroMin(δ_g)
δmin_h = nonzeroMin(δ_h)
η_g = parent(f)(0)
η_h = parent(f)(0)
for i in 1:degree(g)+1
for j in 1:l
if coeff(coeff(δ_g,i-1),j-1)== 0
set_coefficient!(η_g[i],j-1,ϵ*δmin_g)
else
set_coefficient!(η_g[i],j-1,ϵ*δ_g)
end
end
end
for i in 1:degree(h)+1
for j in 1:l
if coeff(coeff(δ_h,i-1),j-1)== 0
set_coefficient!(η_h[i],j-1,ϵ*δmin_h)
else
set_coefficient!(η_h[i],j-1,ϵ*δ_h)
end
end
end
η_g = δ_g
η_h = δ_h
replace_zeros!(η_g,δmin_g,dg,l)
replace_zeros!(η_h,δmin_h,dh,l)
η_g = ϵ*η_g
η_h = ϵ*η_h
@show degree(η_g)
@show degree(η_h)
for i in 0:imax
e = g*h -f
......@@ -357,15 +354,22 @@ function hensel_validate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; ϵ=1e-10
(q,r) = fastdivrem(s*e,h,l,validation = true)
gg = g + t*e + q*g
BivPolMod!(gg,l)
gg = truncate(gg,dg+1)
hh = h + r
BivPolMod!(hh,l)
if (radius(hh)-radius(h) < η_h)&&(radius(gg)-radius(g)<η_g) #todo
hh = truncate(hh,dh+1)
#@show isless_coefficientwise(extract_radius(hh)-extract_radius(h),η_h)
#@show extract_radius(hh)-extract_radius(h)
#@show η_h
#@show isless_coefficientwise(extract_radius(gg)-extract_radius(g),η_g)
if (isless_coefficientwise(extract_radius(hh)-extract_radius(h),η_h))&&(isless_coefficientwise(extract_radius(gg)-extract_radius(g),η_g))
return (g,h)
end
g = gg
h = hh
end
print("no convergence")
print("no convergence\n")
(g,h)
end
"""
......@@ -499,8 +503,8 @@ function validate_a_posteriori_inverse(q::PolyElem{T}, b::PolyElem{S}; ϵ = 1e-1
η = mag_coefficientwise(η)
δ = mag_coefficientwise(δ)
@show η
@show δ
#@show η
#@show δ
#truncate in the second variable for bivariate functions
#if T <: PolyElem
......@@ -510,7 +514,7 @@ function validate_a_posteriori_inverse(q::PolyElem{T}, b::PolyElem{S}; ϵ = 1e-1
#replace zero entries in δ by the smallest nonzero entry
δmin = nonzeroMin(δ)
@show δmin
#@show δmin
e = parent(δ)(0)
for i 1:k+1
if T <: PolyElem
......@@ -566,8 +570,8 @@ end
function validate_a_posteriori_inverse_bivariate(q::PolyElem{T}, b::PolyElem{S}, trunc::Int; ϵ = 1e-10, tmax = 20) where {T <: PolyElem, S<: PolyElem}
db = degree(b)
q_rig = convert_float_to_balls(q)
@show q_rig
@show b
#@show q_rig
#@show b
η = 1-mullow(q_rig,b,db+1)
BivPolMod!(η,trunc)
......@@ -577,8 +581,8 @@ function validate_a_posteriori_inverse_bivariate(q::PolyElem{T}, b::PolyElem{S},
η = mag_coefficientwise(η)
δ = mag_coefficientwise(δ)
@show η
@show δ
#@show η
#@show δ
if δ == 0
return (parent(q)(0),0,"exact solution")
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment