Merge pull request #69 from ederc/gb-qq

GBs over QQ

Author: ederc

Project.toml

@@ -1,7 +1,7 @@
 name = "AlgebraicSolving"
 uuid = "66b61cbe-0446-4d5d-9090-1ff510639f9d"
 authors = ["ederc <[email protected]>", "Mohab Safey El Din <[email protected]", "Rafael Mohr <[email protected]>"]
-version = "0.5.1"
+version = "0.6.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -24,4 +24,4 @@ Random = "1.6"
 StaticArrays = "1"
 Test = "1.6"
 julia = "1.6"
-msolve_jll = "0.600.600"
+msolve_jll = "0.700.100"

docs/Project.toml

@@ -8,4 +8,3 @@ Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 
 [compat]
 Documenter = "0.26"
-msolve_jll = "0.600.501"

src/algorithms/groebner-bases.jl

@@ -1,4 +1,4 @@
-import msolve_jll: libneogb
+import msolve_jll: libneogb, libmsolve
 
 export groebner_basis, eliminate
 
@@ -11,7 +11,7 @@ the first `eliminate` variables.
 
 At the moment the underlying algorithm is based on variants of Faugère's F4 Algorithm.
 
-**Note**: At the moment only ground fields of characteristic `p`, `p` prime, `p < 2^{31}` are supported.
+**Note**: At the moment only ground fields of characteristic `p`, `p` prime, `p < 2^{31}` and the rationals are supported.
 
 # Arguments
 - `I::Ideal{T} where T <: MPolyRingElem`: input generators.
@@ -21,7 +21,9 @@ At the moment the underlying algorithm is based on variants of Faugère's F4 Alg
 - `nr_thrds::Int=1`: number of threads for parallel linear algebra.
 - `max_nr_pairs::Int=0`: maximal number of pairs per matrix, only bounded by minimal degree if `0`.
 - `la_option::Int=2`: linear algebra option: exact sparse-dense (`1`), exact sparse (`2`, default), probabilistic sparse-dense (`42`), probabilistic sparse(`44`).
-- `complete_reduction::Bool=true`: compute a reduced Gröbner basis for `I`
+- `complete_reduction::Bool=true`: compute a reduced Gröbner basis for `I`.
+- `normalize::Bool=false`: normalize generators of Gröbner basis for `I`, only applicable when working over the rationals.
+- `truncate_lifting::Int=0`: truncates the lifting process to given number of elements, only applicable when working over the rationals.
 - `info_level::Int=0`: info level printout: off (`0`, default), summary (`1`), detailed (`2`).
 
 # Examples
@@ -48,6 +50,8 @@ function eliminate(
         max_nr_pairs::Int=0,
         la_option::Int=2,
         complete_reduction::Bool=true,
+        normalize::Bool=false,
+        truncate_lifting::Int=0,
         info_level::Int=0
         )
     if eliminate <= 0
@@ -57,6 +61,8 @@ function eliminate(
                               max_nr_pairs=max_nr_pairs, la_option=la_option,
                               eliminate=eliminate, intersect=intersect,
                               complete_reduction=complete_reduction,
+                              normalize = normalize,
+                              truncate_lifting = truncate_lifting,
                               info_level=info_level)
     end
 end
@@ -67,7 +73,7 @@ end
 Compute a Groebner basis of the ideal `I` w.r.t. to the degree reverse lexicographical monomial ordering using Faugère's F4 algorithm.
 At the moment the underlying algorithm is based on variants of Faugère's F4 Algorithm.
 
-**Note**: At the moment only ground fields of characteristic `p`, `p` prime, `p < 2^{31}` are supported.
+**Note**: At the moment only ground fields of characteristic `p`, `p` prime, `p < 2^{31}` and the rationals are supported.
 
 # Arguments
 - `I::Ideal{T} where T <: MPolyRingElem`: input generators.
@@ -77,7 +83,9 @@ At the moment the underlying algorithm is based on variants of Faugère's F4 Alg
 - `la_option::Int=2`: linear algebra option: exact sparse-dense (`1`), exact sparse (`2`, default), probabilistic sparse-dense (`42`), probabilistic sparse(`44`).
 - `eliminate::Int=0`: size of first block of variables to be eliminated.
 - `intersect::Bool=true`: compute the `eliminate`-th elimination ideal.
-- `complete_reduction::Bool=true`: compute a reduced Gröbner basis for `I`
+- `complete_reduction::Bool=true`: compute a reduced Gröbner basis for `I`.
+- `normalize::Bool=false`: normalize generators of Gröbner basis for `I`, only applicable when working over the rationals.
+- `truncate_lifting::Int=0`: truncates the lifting process to given number of elements, only applicable when working over the rationals.
 - `info_level::Int=0`: info level printout: off (`0`, default), summary (`1`), detailed (`2`).
 
 # Examples
@@ -111,6 +119,8 @@ function groebner_basis(
         eliminate::Int=0,
         intersect::Bool=true,
         complete_reduction::Bool=true,
+        normalize::Bool=false,
+        truncate_lifting::Int=0,
         info_level::Int=0
         )
 
@@ -119,6 +129,8 @@ function groebner_basis(
                              max_nr_pairs = max_nr_pairs, la_option = la_option,
                              eliminate = eliminate, intersect = intersect,
                              complete_reduction = complete_reduction,
+                             normalize = normalize,
+                             truncate_lifting = truncate_lifting,
                              info_level = info_level)
     end
 end
@@ -132,6 +144,8 @@ function _core_groebner_basis(
         eliminate::Int=0,
         intersect::Bool=true,
         complete_reduction::Bool=true,
+        normalize::Bool=false,
+        truncate_lifting::Int=0,
         info_level::Int=0
         )
 
@@ -146,11 +160,17 @@ function _core_groebner_basis(
 
     mon_order       = 0
     elim_block_size = eliminate
+    if elim_block_size >= nr_vars
+        error("Number of variables to be eliminated is bigger than number of variables in ring.")
+    end
     reduce_gb       = Int(complete_reduction)
 
     # convert ideal to flattened arrays of ints
-    if !(is_probable_prime(field_char))
-        error("At the moment we only support finite fields.")
+
+    if !(field_char == 0)
+        if !(is_probable_prime(field_char))
+            error("At the moment we only support finite fields or the rationals.")
+        end
     end
 
     # nr_gens might change if F contains zero polynomials
@@ -161,29 +181,52 @@ function _core_groebner_basis(
     gb_exp = Ref(Ptr{Cint}(0))
     gb_cf  = Ref(Ptr{Cvoid}(0))
 
-    nr_terms  = ccall((:export_f4, libneogb), Int,
-        (Ptr{Nothing}, Ptr{Cint}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cvoid}},
-        Ptr{Cint}, Ptr{Cint}, Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cint, Cint,
-        Cint, Cint, Cint, Cint, Cint, Cint, Cint),
-        cglobal(:jl_malloc), gb_ld, gb_len, gb_exp, gb_cf, lens, exps, cfs,
-        field_char, mon_order, elim_block_size, nr_vars, nr_gens, initial_hts,
-        nr_thrds, max_nr_pairs, 0, la_option, reduce_gb, 0, info_level)
+    if field_char == 0
+        nr_terms  = ccall((:export_groebner_qq, libmsolve), Int,
+            (Ptr{Nothing}, Ptr{Cint}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cvoid}},
+            Ptr{Cint}, Ptr{Cint}, Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cint, Cint,
+            Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint),
+            cglobal(:jl_malloc), gb_ld, gb_len, gb_exp, gb_cf, lens, exps, cfs,
+            field_char, mon_order, elim_block_size, nr_vars, nr_gens, initial_hts,
+            nr_thrds, max_nr_pairs, 0, la_option, reduce_gb, 0, truncate_lifting, info_level)
+
+    else
+        nr_terms  = ccall((:export_f4, libneogb), Int,
+            (Ptr{Nothing}, Ptr{Cint}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cvoid}},
+            Ptr{Cint}, Ptr{Cint}, Ptr{Cvoid}, Cint, Cint, Cint, Cint, Cint, Cint,
+            Cint, Cint, Cint, Cint, Cint, Cint, Cint),
+            cglobal(:jl_malloc), gb_ld, gb_len, gb_exp, gb_cf, lens, exps, cfs,
+            field_char, mon_order, elim_block_size, nr_vars, nr_gens, initial_hts,
+            nr_thrds, max_nr_pairs, 0, la_option, reduce_gb, 0, info_level)
+    end
 
     # convert to julia array, also give memory management to julia
     jl_ld   = gb_ld[]
     jl_len  = Base.unsafe_wrap(Array, gb_len[], jl_ld)
     jl_exp  = Base.unsafe_wrap(Array, gb_exp[], nr_terms*nr_vars)
-    ptr     = reinterpret(Ptr{Int32}, gb_cf[])
-    jl_cf   = Base.unsafe_wrap(Array, ptr, nr_terms)
 
-    if intersect == true
-        basis = _convert_finite_field_array_to_abstract_algebra(
-            jl_ld, jl_len, jl_cf, jl_exp, R, eliminate)
+    #  coefficient handling depending on field characteristic
+    if field_char == 0
+        ptr     = reinterpret(Ptr{BigInt}, gb_cf[])
+        jl_cf   = [QQFieldElem(unsafe_load(ptr, i)) for i in 1:nr_terms]
     else
-        basis = _convert_finite_field_array_to_abstract_algebra(
-            jl_ld, jl_len, jl_cf, jl_exp, R, 0)
+        ptr     = reinterpret(Ptr{Int32}, gb_cf[])
+        jl_cf   = Base.unsafe_wrap(Array, ptr, nr_terms)
+    end
+
+    #  shall eliminated variables be removed?
+    if !intersect
+        eliminate = 0
     end
 
+    #  convert to basis
+    if field_char == 0
+        basis = _convert_rational_array_to_abstract_algebra(
+            jl_ld, jl_len, jl_cf, jl_exp, R, normalize, eliminate)
+    else
+        basis = _convert_finite_field_array_to_abstract_algebra(
+            jl_ld, jl_len, jl_cf, jl_exp, R, eliminate)
+    end
     ccall((:free_f4_julia_result_data, libneogb), Nothing ,
           (Ptr{Nothing}, Ptr{Ptr{Cint}}, Ptr{Ptr{Cint}},
            Ptr{Ptr{Cvoid}}, Int, Int),

src/interfaces/nemo.jl

@@ -124,3 +124,54 @@ function _convert_finite_field_array_to_abstract_algebra(
 
     return basis
 end
+
+function _convert_rational_array_to_abstract_algebra(
+        bld::Int32,
+        blen::Vector{Int32},
+        bcf::Vector{QQFieldElem},
+        bexp::Vector{Int32},
+        R::MPolyRing,
+        normalize::Bool=false,
+        eliminate::Int=0
+        )
+
+    #  Note: Over QQ msolve already returns the eliminated basis
+    #  only in order to save memory due to possible huge
+    #  coefficient sizes.
+    if characteristic(R) != 0
+        error("We assume QQ as base field here.")
+    end
+
+    nr_gens = bld
+    nr_vars = nvars(R)
+    CR      = coefficient_ring(R)
+
+    basis = (typeof(R(0)))[]
+
+    len   = 0
+
+    for i in 1:nr_gens
+        if bcf[len+1] == 0
+            push!(basis, R(0))
+        else
+            g  = MPolyBuildCtx(R)
+            lc = bcf[len+1]
+
+            if normalize && lc != 1
+                for j in 1:blen[i]
+                    push_term!(g, bcf[len+j]/lc,
+                               convert(Vector{Int}, bexp[(len+j-1)*nr_vars+1:(len+j)*nr_vars]))
+                end
+            else
+                for j in 1:blen[i]
+                    push_term!(g, bcf[len+j],
+                               convert(Vector{Int}, bexp[(len+j-1)*nr_vars+1:(len+j)*nr_vars]))
+                end
+            end
+            push!(basis, finish(g))
+        end
+        len +=  blen[i]
+    end
+
+    return basis
+end

test/algorithms/groebner-bases.jl

@@ -1,8 +1,4 @@
 @testset "Algorithms -> Gröbner bases" begin
-    R, (x,y,z) = polynomial_ring(QQ,["x","y","z"], internal_ordering=:degrevlex)
-    F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
-    #= not a finite field =#
-    @test_throws ErrorException groebner_basis(Ideal(F))
     R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"], internal_ordering=:degrevlex)
     I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     G = groebner_basis(I)
@@ -38,6 +34,29 @@
     I = Ideal([R(0)])
     G = groebner_basis(I)
     @test G == [R(0)]
+
+    R, (x1, x2) = polynomial_ring(QQ, ["x1", "x2"])
+    I = Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
+    G = groebner_basis(I)
+    H = MPolyRingElem[
+        3*x1 - 2*x2 - 5
+        4*x2^2 + 20*x2 + 3802951800684688204490109616153
+        ]
+    @test G == H
+    J = Ideal([3*x1^2 + ZZRingElem(2)^100, 2*x1*x2 + 5*x1 + ZZRingElem(2)^100])
+    G = groebner_basis(J, normalize=true)
+    H = MPolyRingElem[
+        x1 - 2//3*x2 - 5//3
+        x2^2 + 5*x2 + 3802951800684688204490109616153//4
+        ]
+    @test G == H
+    R, (x,y,z) = polynomial_ring(QQ,["x","y","z"], internal_ordering=:degrevlex)
+    I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
+    G = eliminate(I, 2)
+    H = MPolyRingElem[
+        84*z^4 - 40*z^3 + z^2 + z
+    ]
+    @test G == H
 end
 
 @testset "Algorithms -> Sig Gröbner bases" begin

test/interfaces/nemo.jl

@@ -9,13 +9,19 @@
     cmp = (Int32[1], BigInt[1], Int32[1, 0], 1) 
     @test AlgebraicSolving._convert_to_msolve(I.gens) == cmp
     for _GF in [GF, AlgebraicSolving.Nemo.Native.GF]
-      R, (x,y,z) = polynomial_ring(_GF(2147483659),["x","y","z"], internal_ordering=:degrevlex)
-      F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
-      # prime is bigger than 2^31, should throw an error
-      @test_throws ErrorException AlgebraicSolving._convert_to_msolve(F)
-      R, (x,y,z) = polynomial_ring(_GF(101),["x","y","z"], internal_ordering=:degrevlex)
-      F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
-      res = AlgebraicSolving._convert_to_msolve(F)
-      @test AlgebraicSolving._convert_finite_field_array_to_abstract_algebra(Int32(3), res[1], res[2], res[3], R) == F
+        R, (x,y,z) = polynomial_ring(_GF(2147483659),["x","y","z"], internal_ordering=:degrevlex)
+        F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
+        # prime is bigger than 2^31, should throw an error
+        @test_throws ErrorException AlgebraicSolving._convert_to_msolve(F)
+        R, (x,y,z) = polynomial_ring(_GF(101),["x","y","z"], internal_ordering=:degrevlex)
+        F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
+        @show F
+        res = AlgebraicSolving._convert_to_msolve(F)
+        @test AlgebraicSolving._convert_finite_field_array_to_abstract_algebra(Int32(3), res[1], res[2], res[3], R) == F
     end
+    R, (x,y,z) = polynomial_ring(QQ,["x","y","z"], internal_ordering=:degrevlex)
+    F = [x^2+1-3, x*y-z, x*z^2-3*y^2]
+    res = AlgebraicSolving._convert_to_msolve(F)
+    res_qq =AlgebraicSolving.QQFieldElem[res[2][i] for i in 1:2:length(res[2])]
+    @test AlgebraicSolving._convert_rational_array_to_abstract_algebra(Int32(3), res[1], res_qq, res[3], R) == F
 end