Merge branch 'main' of https://github.com/UCSD-computational-number-theory/DeRham.jl

Author: rybplayer

src/ControlledReduction.jl

@@ -14,6 +14,7 @@ using Memoize
 #include("NemoAdditions.jl")
 
 include("Utils.jl")
+include("GradedExpCache.jl")
 include("FindMonomialBasis.jl")
 include("SlopesPolygon.jl")
 include("PolynomialWithPole.jl")
@@ -31,6 +32,7 @@ include("ExamplePolynomials.jl")
 
 include("ZetaFunction.jl")
 include("PointCounts.jl")
+include("NewtonPolygon.jl")
 
 # TODO: export Zeta Function functions
 

src/DeRham.jl

@@ -0,0 +1,109 @@
+
+abstract type GradedExpCache end
+
+#TODO: convert to SVector, that might help with allocation pressure
+#using StaticArrays
+
+struct PolyExpCache <: GradedExpCache
+    n::Int
+    termorder::Symbol
+    exp_vec_pieces::Vector{Union{Vector{Vector{Int64}},Nothing}}
+    rev_look_pieces::Vector{Union{Dict{Vector{Int64},Int},Nothing}}
+
+    PolyExpCache(n,termorder) = new(n,termorder,[],[])
+end
+
+
+"""
+Gets the exponent vector cache for the degree d,
+returning it as a vector.
+
+This throws an error if the cached degree doesn't exist.
+"""
+function Base.getindex(c::PolyExpCache,d::Int)
+    if c.exp_vec_pieces[d] == nothing
+        errormsg = "This PolyExpCache does not have degree=$d stored. "* 
+                   "This probably means that the PolyExpCache hasn't been "* 
+                   "fully initialized, or there is a bug. If you want to "* 
+                   "initialize this cache to work for degree $d, use the "*
+                   "`get_forward` function."
+        throw(ArgumentError(errormsg))
+    end
+    c.exp_vec_pieces[d]
+end
+
+function Base.getindex(c::PolyExpCache,d::Int,kind::Symbol)
+    if kind == :forward
+        c[d]
+    elseif kind == :reverse
+
+        if c.rev_look_pieces[d] == nothing
+            errormsg = "This PolyExpCache does not have degree=$d stored for "*
+                       "reverse lookup. "* 
+                       "This probably means that the PolyExpCache hasn't been "* 
+                       "fully initialized, or there is a bug. If you want to "* 
+                       "initialize this cache to work for degree $d, use the "*
+                       "`get_forward` function."
+            throw(ArgumentError(errormsg))
+        end
+        c.rev_look_pieces[d]
+    else
+        throw(ArgumentError("Invalid option for indexing PolyExpCache"))
+    end
+end
+
+
+function Base.show(io::IO, c::PolyExpCache)
+    msg = """PolyExpCache for the graded ring k[x_1, ..., x_$(c.n)]
+    Term order: $(c.termorder)
+    Terms cached for degrees: $(cached_degrees(c))
+    Reverse lookups cached for degrees: $(cached_reverse_lookups(c))
+    """
+    print(io, msg)
+end
+
+cached_degrees(c::PolyExpCache) = findall(c.exp_vec_pieces .!= nothing)
+cached_reverse_lookups(c::PolyExpCache) = findall(c.rev_look_pieces .!= nothing)
+
+function get_forward(c::PolyExpCache,d)
+    if c.exp_vec_pieces[d] == nothing
+        generate_degree_forward(c,d)
+    end
+    c.exp_vec_pieces[d]
+end
+
+function generate_degree_forward(c::PolyExpCache,d)
+    l = length(c.exp_vec_pieces) 
+    if l < d
+        append!(c.exp_vec_pieces,fill(nothing,d - l))
+    end
+    c.exp_vec_pieces[d] = gen_exp_vec(c.n,d,c.termorder)
+end
+
+function generate_degree_reverse(c::PolyExpCache,d)
+    evs = get_forward(c,d)
+
+    l = length(c.rev_look_pieces) 
+    if l < d
+        append!(c.rev_look_pieces,fill(nothing,d - l))
+    end
+    c.rev_look_pieces[d] = Dict(evs[i] => i for i in 1:length(evs))
+end
+
+function PolyExpCache(n,termorder,degrees_to_prefill,reverse_to_prefill)
+    c = PolyExpCache(n,termorder)
+    for d in degrees_to_prefill
+        generate_degree_forward(c,d)
+    end
+
+    for d in reverse_to_prefill
+        generate_degree_reverse(c,d)
+    end
+    c
+end
+
+function PolyExpCache(n,termorder,max_degree_prefill)
+    PolyExpCache(n,termorder,1:max_degree_prefill,1:max_degree_prefill)
+end
+
+

src/GradedExpCache.jl

@@ -0,0 +1,64 @@
+"""
+    vp(a, p)
+    Returns the p-adic valuation of a 
+    
+    INPUTS: 
+    * "a" -- integer
+    * "p" -- integer, a prime number 
+"""
+function vp(a, p)
+    @assert is_prime(p)
+    if a == 0 
+        return Inf 
+    end 
+
+    e = 0 
+    while mod(a, p) == 0
+        a = div(a, p)
+        e = e + 1
+    end
+
+    return e
+end 
+
+"""
+    newton_polygon(f, p)
+    Returns the Newton polygon of f 
+
+    INPUTS: 
+    * "f" -- list, coefficients of a polynomial, constant term first
+    * "p" -- integer, a prime number
+"""
+function newton_polygon(f, p)
+    @assert is_prime(p)
+
+    pts = []
+    for (i, a) in enumerate(f)
+        if a != 0
+            push!(pts, (i-1, vp(a, p)))
+        end 
+    end 
+
+    # compute the lower convex hull of points
+    hull = Tuple{Int64, Int64}[]
+    for pt in pts
+        while length(hull) >= 2
+            (x1, y1) = hull[end-1]
+            (x2, y2) = hull[end]
+            (x3, y3) = pt
+            # Compute cross product to decide whether pt is "below" the last segment.
+            cp = (x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)
+            if cp <= 0
+                pop!(hull)
+            else
+                break
+            end
+        end
+        push!(hull, pt)
+    end
+    return SlopesPolygon(hull)
+   #return hull
+end 
+
+
+

src/NewtonPolygon.jl

@@ -103,18 +103,22 @@ end
 """
 Given an array of vertices
 """
-function SlopesPolygon(vertices::Array{Tuple{Int,Int}})
+#function SlopesPolygon(vertices::Array{Tuple{Int,Int}})
+function SlopesPolygon(vertices::Vector{Tuple{Int64, Int64}})
     #TODO implement
 
     n = length(vertices)-1
     slopelengths = zeros(Int,n)
     slopes = zeros(Rational{Int},n)
-    for i = 2:n
+    for i = 2:n+1
         slopevec = vertices[i] .- vertices[i-1]
-        push!(slopes, slopevec[2] // slopevec[1])
-        push!(slopelengths,slopevec[2])
+        #push!(slopes, slopevec[2] // slopevec[1])
+        #push!(slopelengths,slopevec[2])
+        slopes[i-1] = slopevec[2] // slopevec[1]
+        slopelengths[i-1] = slopevec[2]
     end
-
+    println(slopelengths)
+    println(slopes)
     SlopesPolygon(slopes,slopelengths,values(slopes,slopelengths)...)
 end
 

src/SlopesPolygon.jl

@@ -28,7 +28,7 @@ partials = [ derivative(polynomial, i) for i in 1:(n+1) ]
 # TODO: Ensure partials are not all zero.
 
 """
-t   gen_exp_vec(n, d, order)
+   gen_exp_vec(n, d, order)
 
 Returns all nonnegative integer lists of length n who entires sum to d
 

src/Utils.jl

@@ -33,6 +33,48 @@ end
 
 default_params() = ZetaFunctionParams(false,false,:costachunks,:invlex,true,false,false,false)
 
+"""
+Give the minimal PolyExpCache needed for controlled reduction
+as in Prop 1.15 of Costa's thesis.
+
+n in this function is the number of variables minus one,
+i.e. the weight of the affine Monsky-Washnitzer homology module.
+
+the d we need are
+
+FORWARD:
+h*d - n - 1 for h in 1:n (if such values are positive)
+d*n - n 
+d*n - n - 1
+d*n - n - length(S),
+d*n - n - length(S) + 1,
+d*n - n + d - length(S) + 1,
+d
+
+REVERSE:
+d
+n*d - n
+
+"""
+function controlled_reduction_cache(n,d,S,termorder)
+    degsforward = [d*n - n,
+                   d*n - n - 1,
+                   d*n - n - length(S),
+                   d*n - n - length(S) + 1,
+                   d*n - n + d - length(S),
+                   d]
+    for h in 1:n
+        if 0 < h*d - n - 1
+            append!(degsforward,h*d - n - 1)
+        end
+    end
+    #append!(degsforward,[h*d - n - 1 for h in 1:n])
+
+    degsreverse = [d,n*d - n]
+
+    PolyExpCache(n+1,termorder,degsforward,degsreverse)
+end
+
 """
     compute_frobenius_matrix(n,d,Reductions,T)
 
@@ -44,15 +86,18 @@ INPUTS:
 * "N" -- integer, series precision
 * "Reductions" -- output of computeReductionOfTransformLA
 * "T" -- output of computeT
+* "Basis" -- ?????
+* "params" -- the ControlledReductionParamaters
+* "cache" -- the GradedExpCache used for this controlled reduction
 """
-function compute_frobenius_matrix(n, p, d, N_m, Reductions, T, Basis, params)
+function compute_frobenius_matrix(n, p, d, N_m, Reductions, T, Basis, params, cache)
     verbose = params.verbose
     termorder = params.termorder
     (9 < verbose) && println("Terms after controlled reduction: $Reductions")
     R = parent(T[1,1])
     frob_mat_temp = []
     denomArray = []
-    ev = gen_exp_vec(n+1,d*n-n-1,termorder)
+    ev = cache[d*n-n-1]#gen_exp_vec(n+1,d*n-n-1,termorder)
     VS = matrix_space(R,length(ev),1)
     for i in 1:length(Reductions)
         e = Basis[i][2] # pole order of basis element 
@@ -281,9 +326,14 @@ function zeta_function(f; S=[-1], verbose=false, givefrobmat=false, algorithm=:c
     n = nvars(PR) - 1
     d = total_degree(f)
     R = coefficient_ring(PR)
+    if S == [-1]
+        S = collect(0:n)
+    end
 
     params = ZetaFunctionParams(verbose,givefrobmat,algorithm,termorder,vars_reversed,fastevaluation,always_use_bigints,use_gpu)
 
+    cache = controlled_reduction_cache(n,d,S,termorder)
+
     (9 < verbose) && println("Working with a degree $d hypersurface in P^$n")
 
     basis = compute_monomial_bases(f, R, PR, params.termorder) # basis of cohomology 
@@ -317,9 +367,6 @@ function zeta_function(f; S=[-1], verbose=false, givefrobmat=false, algorithm=:c
     precisionringpoly, pvars = polynomial_ring(precisionring, ["x$i" for i in 0:n])
 
     #S = SmallestSubsetSmooth.smallest_subset_s_smooth(fLift,n)
-    if S == [-1]
-        S = collect(0:n)
-    end
 
     if (0 < verbose)
         println("Starting linear algebra problem")
@@ -409,7 +456,7 @@ function zeta_function(f; S=[-1], verbose=false, givefrobmat=false, algorithm=:c
     #end
     #TODO: check which algorithm we're using
     (2 < verbose) && println("Pseudo inverse matrix:\n$pseudo_inverse_mat")
-    Reductions = reducetransform(FBasis, N_m, S, fLift, pseudo_inverse_mat, p,  params) 
+    Reductions = reducetransform(FBasis, N_m, S, fLift, pseudo_inverse_mat, p,  params, cache) 
     (1 < verbose) && println(Reductions)
     #return Reductions
     #if (1 < verbose)
@@ -420,9 +467,10 @@ function zeta_function(f; S=[-1], verbose=false, givefrobmat=false, algorithm=:c
     #        println()
     #    end 
     #end 
-    ev = gen_exp_vec(n+1,n*d-n-1,termorder)
+    ev = cache[n*d - n - 1] 
+    #ev = gen_exp_vec(n+1,n*d-n-1,termorder)
     (9 < verbose) && println(convert_p_to_m([Reductions[1][1][1],Reductions[2][1][1]],ev))
-    FM = compute_frobenius_matrix(n, p, d, N_m, Reductions, T, Basis, params)
+    FM = compute_frobenius_matrix(n, p, d, N_m, Reductions, T, Basis, params, cache)
     # display(FM)
     (9 < verbose) && println("The Frobenius matrix is $FM")
 

src/ZetaFunction.jl

@@ -10,3 +10,8 @@ n;p;f;zeta_function
 3;19;x1^4*x2 + x1*x3^4 + x2^4*x3;ZZRingElem[]
 3;11;5*x1^6*x2 + 3*x1^6*x3 + 7*x1^5*x2^2 + 9*x1^5*x2*x3 + x1^5*x3^2 + 8*x1^4*x2^3 + x1^4*x2^2*x3 + 5*x1^4*x2*x3^2 + 6*x1^3*x2^4 + 6*x1^3*x2^3*x3 + 9*x1^3*x2^2*x3^2 + 8*x1^3*x2*x3^3 + 5*x1^3*x3^4 + 5*x1^2*x2^5 + 2*x1^2*x2^4*x3 + 7*x1^2*x2^3*x3^2 + 4*x1^2*x2*x3^4 + 10*x1^2*x3^5 + 7*x1*x2^6 + 3*x1*x2^5*x3 + 5*x1*x2^4*x3^2 + 8*x1*x2^3*x3^3 + 2*x1*x2^2*x3^4 + 3*x1*x2*x3^5 + 5*x1*x3^6 + 2*x2^7 + 4*x2^6*x3 + 2*x2^5*x3^2 + 2*x2^4*x3^3 + 9*x2^3*x3^4 + 7*x2^2*x3^5 + 4*x2*x3^6 + 2*x3^7;ZZRingElem[4177248169415651,-1518999334332964,655931530734689,-200859416110144,34237400473320,-8455600419926,858292959524,-267305524607,68887149485,-21921295814,5806049601,2472191414,-1995260839,889046048,-336993701,104438702,-30635791,7347488,-1499069,168854,36051,-12374,3535,-1247,364,-326,120,-64,19,-4,1]
 3;17;x1^3 + x2^3 + x3^3;ZZRingElem[]
+4;11;x1^3 + x1^2*x2 + x1^2*x3 + x1^2*x4 + x1*x2^2 + x1*x2*x3 + x1*x2*x4 + x1*x3^2 + x1*x3*x4 + x1*x4^2 + x2^3 + x2^2*x3 + x2^2*x4 + x2*x3^2 + x2*x3*x4 + x2*x4^2 + x3^3 + x3^2*x4 + x3*x4^2 + x4^3;ZZRingElem[]
+3;11;x1^3*x2 + 6*x1^3*x3 + 3*x1^2*x2^2 + 5*x1^2*x3^2 + 5*x1*x2^3 + 10*x1*x2^2*x3 + 9*x1*x3^3 + 8*x2^4 + 7*x2^3*x3 + 7*x2^2*x3^2 + x2*x3^3 + 3*x3^4;ZZRingElem[]
+4;13;9*x1^3 + 9*x1^2*x2 + 3*x1^2*x3 + 2*x1^2*x4 + 7*x1*x2^2 + x1*x2*x3 + 12*x1*x2*x4 + 3*x1*x3^2 + 6*x1*x3*x4 + 2*x1*x4^2 + 9*x2^3 + 12*x2^2*x3 + 3*x2^2*x4 + 3*x2*x3^2 + 3*x2*x3*x4 + 3*x2*x4^2 + 12*x3^3 + 3*x3^2*x4 + 9*x3*x4^2;ZZRingElem[]
+4;13;x1^3 + x1^2*x2 + x1^2*x3 + x1^2*x4 + x1*x2^2 + x1*x2*x3 + x1*x2*x4 + x1*x3^2 + x1*x3*x4 + x1*x4^2 + x2^3 + x2^2*x3 + x2^2*x4 + x2*x3^2 + x2*x3*x4 + x2*x4^2 + x3^3 + x3^2*x4 + x3*x4^2 + x4^3;ZZRingElem[]
+3;13;x1^3 + x1*x3^2 + x2^3;ZZRingElem[13,-2,1]