worked on a posteriori validation for Hensel

functions.jl

@@ -293,6 +293,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 gstar
     #@show r
     hstar = h+r
     BivPolMod!(hstar, 2deg)
@@ -697,7 +698,12 @@ function validate_a_posteriori_inverse_bivariate(q::PolyElem{T}, b::PolyElem{S},
     (r, tmax, string("no convergence in ", tmax, " steps"))
 end
 
-function validate_Hensel(l::Int, f::PolyElem, g::PolyElem, h::PolyElem)
+"""
+    validate_Hensel(l::Int, f::PolyElem, g::PolyElem, h::PolyElem)
+
+Compute a numeric approximation of Hensel and then validate the results, returning interval valued polynomial results for g and h
+"""
+function validate_Hensel_a_posteriori(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; ϵ = 1e-10, tmax = 20)
     #convert f,g,h to polynomials in Float64
     fr = convert_balls_to_Float(f)
     gr = convert_balls_to_Float(g)
@@ -706,34 +712,88 @@ function validate_Hensel(l::Int, f::PolyElem, g::PolyElem, h::PolyElem)
     dh = degree(h)
     #compute numerical solution
     (gstar, hstar, s, t )= henseltruncate(l,fr,gr,hr,validation = false)
-    gstar = convert_float_to_balls(gstar)
-    hstar = convert_float_to_balls(hstar)
+    @show degree(gstar)
+    @show degree(hstar)
+    @show gstar
+    igstar = convert_float_to_balls(gstar)
+    ihstar = convert_float_to_balls(hstar)
+    s = convert_float_to_balls(s)
+    t = convert_float_to_balls(t)
     #compute δ
-    (q,R) = fastdivrem(s*(gstar*hstar-f),hstar)
-    δg = t*e + q*g
-    BivPolMod!(δg,l)
-    δh = R
+    (q,R) = fastdivrem(s*(igstar*ihstar-f),ihstar, l, validation = true)
+    δg = t*(igstar*ihstar-f) + q*g
+    δg = mag_coefficientwise(δg)
+    #BivPolMod!(δg,l)
+    δg = trunc(δg, dg+1,l)
+    δh = mag_coefficientwise(R)
+    #BivPolMod!(δh,l)
+    δh = trunc(δh, dh,l)
     δgmin = nonzeroMin(δg)
-    ηg = deepcopy(δg)
-    replace_zeros!(ηg,δgmin, degree(g), l)
+    ηg = ϵ*deepcopy(δg)
+    replace_zeros!(ηg,ϵ*δgmin, degree(g), l)
     δhmin = nonzeroMin(δh)
-    ηh = deepcopy(δh)
-    replace_zeros!(ηh, δhmin, degree(h), l)
+    ηh = ϵ*deepcopy(δh)
+    replace_zeros!(ηh, ϵ*δhmin, degree(h), l)
     δh = δh + ηh
     δg = δg + ηg
+    @show δg
+    @show δh
+
+    invrevh = newtoninversion(ihstar,l)
+    u = s*igstar + t*ihstar -1
 
     #validation
     rg = parent(g)(0)
     rh = parent(h)(0)
     for t in 1:tmax
-        (Λg, Λh) = Λ(rg,rh)
-        rrg = Λg*rg + δg
+        (Λ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)
-        rrh = Λh*rh + δh
-        rrh = trunc(rrh, dh+1, l)
+        @show degree(rrg)
+        rrh = mag_coefficientwise(Λh*rh + δh)
+        rrh = trunc(rrh, dh, l)
+        @show degree(rrh)
         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)) 
+            return(gstar,hstar)       
         end  
+        rg = rrg
+        rh = rrh
     end
+    print("no convergence after ", tmax, " steps\n")
+    (rg,rh)
+end
+
+"""
+    Lambda(rg,rh, s,t, ϵ, g, h, invrevh)
+
+Helper function for the validation of Hensel.
+# Arguments
+- `rg,rh` current iteratives of the radii of g and handle
+- `s,t,g,h` our initial numerical approximations
+- `invrevh` rev(h)^-1 so it can be computed only once and reused in following iterations, validated 
+- `ϵ` s*g+t*h = 1 + ϵ
+"""
+function Λ(rg,rh, s,t, ϵ, g, h, invrevh)
+    rrg = convert_balls_to_Float(rg, rounding = RoundUp)
+    rrh = convert_balls_to_Float(rh, rounding = RoundUp)
+    irg = convert_polynomial_to_intervals(parent(rrg)(0),rrg)
+    irh = convert_polynomial_to_intervals(parent(rrh)(0),rrh)
+    #Λg
+    a = (2*s*irg + 1 + ϵ)*irh
+    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)
+    #Λ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)
+    (Λg, Λh)
 end
\ No newline at end of file

practicalfunctions.jl

@@ -431,7 +431,7 @@ function replace_zeros!(p::PolyElem{T}, d, dp = -1, k = -1) where T
         deg = k        
     end
     @show deg
-    for i in 0:deg
+    for i in 0:deg-1
         if T <: PolyElem
             set_coefficient!(p,i,replace_zeros!(coeff(p,i),d,dp,k))
         elseif coeff(p,i) == 0