Validated Hensel

functions.jl

@@ -302,6 +302,72 @@ function henseltruncate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; validatio
     (g,h)
 end
 
+function hensel_validate(l::Int, f::PolyElem, g::PolyElem, h::PolyElem; ϵ=1e-10, imax = 20)
+    #convert input polynomials to real polynomials in Float64
+    gr = convert_balls_to_Float(g)
+    hr = convert_balls_to_Float(h)
+    fr = convert_balls_to_Float(f)
+    #compute Hensel with Float64
+    (g,h,s,t) = henseltruncate(l,fr,gr,hr)
+    #convert resulting polynomials to thin intervals
+    g = convert_float_to_balls(g)
+    h = convert_float_to_balls(h)
+    s = convert_float_to_balls(s)
+    t = convert_float_to_balls(t)
+
+    #compute η
+    e = g*h -f
+    (q,r) = fastdivrem(s*e,h,l,validation = true)
+    δ_g = t*e +q*g
+    BivPolMod!(δ_g,l)
+    abs_coefficientwise!(δ_g)
+    δ_g = mag_coefficientwise(δ_g)
+    @show δ_g
+    δ_h = r
+    BivPolMod!(δ_h,l)
+    abs_coefficientwise!(δ_h)
+    δ_h = mag_coefficientwise(δ_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
+
+    for i in 0:imax
+        e = g*h -f
+        BivPolMod!(e,l)
+        (q,r) = fastdivrem(s*e,h,l,validation = true)
+        gg = g + t*e + q*g
+        BivPolMod!(gg,l)
+        hh = h + r
+        BivPolMod!(hh,l)
+        if (radius(hh)-radius(h) < η_h)&&(radius(gg)-radius(g)<η_g) #todo
+            return (g,h)            
+        end
+        g = gg
+        h = hh        
+    end
+    print("no convergence")
+end
+
 """
     subtract(p::Nemo.AbstractAlgebra.Generic.Poly{QQPolyRingElem}, q::Nemo.AbstractAlgebra.Generic.Poly{Nemo.AbstractAlgebra.Generic.Poly{Float64}})