Merge pull request #96 from rprebet/dimfix

Update I.dim when dimension computed

Author: mohabsafey

.gitignore

@@ -1,3 +1,4 @@
 *~
 docs/build/*
 docs/make-dep.jl
+Manifest.toml

Project.toml

@@ -24,4 +24,4 @@ Random = "1.6"
 StaticArrays = "1"
 Test = "1.6"
 julia = "1.6"
-msolve_jll = "0.700.200"
+msolve_jll = "0.700.501"

src/algorithms/dimension.jl

@@ -18,35 +18,44 @@ julia> dimension(I)
 ```
 """
 function dimension(I::Ideal{T}) where T <: MPolyRingElem
-    
-    gb = get(I.gb, 0, groebner_basis(I, complete_reduction = true))
+
+    !isnothing(I.dim) && return I.dim
+    gb = get!(I.gb, 0) do
+        groebner_basis(I, complete_reduction = true)
+    end
     R = parent(first(gb))
-    res = [trues(ngens(R))]
 
+    res = Set([trues(ngens(R))])
     lead_exps = (_drl_lead_exp).(gb)
-    for lexp in lead_exps 
-        to_del = Int[]
-        new_miss = BitVector[]
-        for (i, mis) in enumerate(res)
-            nz_exps_inds = findall(e -> !iszero(e), lexp)
-            ind_var_inds = findall(mis)
-            if issubset(nz_exps_inds, ind_var_inds)
-                for j in nz_exps_inds
+    for lexp in lead_exps
+        nz_exps = (!iszero).(lexp)
+        nz_exps_ind = findall(nz_exps)
+        next_res = Set{BitVector}()
+        for mis in res
+            if all_lesseq(nz_exps, mis)
+                @inbounds for j in nz_exps_ind
                     new_mis = copy(mis)
                     new_mis[j] = false
-                    push!(new_miss, new_mis)
+                    push!(next_res, new_mis)
                 end
-                push!(to_del, i)
+            else
+                push!(next_res, mis)
             end
         end
-        deleteat!(res, to_del)
-        append!(res, new_miss)
-        unique!(res)
+        res = next_res
     end
 
-    length(res) == 0 && return -1
-    max_length = maximum(mis -> length(findall(mis)), res)
-    return max_length
+    I.dim = isempty(res) ? -1 : maximum(sum, res)
+    return I.dim
+end
+
+function all_lesseq(a::BitVector, b::BitVector)::Bool
+    @inbounds for i in eachindex(a)
+        if a[i] && !b[i]
+            return false
+        end
+    end
+    return true
 end
 
 function _drl_exp_vector(u::Vector{Int})
@@ -55,6 +64,6 @@ end
 
 function _drl_lead_exp(p::MPolyRingElem)
     exps = collect(Nemo.exponent_vectors(p))
-    _, i = findmax((u -> _drl_exp_vector(u)).(exps))
+    _, i = findmax(_drl_exp_vector.(exps))
     return exps[i]
-end
+end

src/algorithms/groebner-bases.jl

@@ -245,6 +245,5 @@ function _core_groebner_basis(
           cglobal(:jl_free), gb_len, gb_exp, gb_cf, jl_ld, field_char)
 
     I.gb[eliminate] = basis
-
-    return basis
+    return I.gb[eliminate]
 end

src/algorithms/solvers.jl

@@ -123,17 +123,17 @@ function _core_msolve(
     jl_dquot      = res_dquot[]
     jl_nb_sols    = nb_sols[]
 
-    # set dimension
-    I.dim = jl_dim
     if jl_dim > 0
         error("Dimension of ideal is greater than zero, no solutions provided.")
     end
+    I.dim = 0
 
     # get rational parametrization
     jl_len  = Base.unsafe_wrap(Array, res_len[], jl_ld)
     nterms  = 0
 
     if jl_dquot == 0
+        I.dim = -1
         C,x = polynomial_ring(QQ,"x")
         I.rat_param = RationalParametrization(Symbol[], ZZRingElem[], C(-1), C(-1), PolyRingElem[])
         I.real_sols = QQFieldElem[]
@@ -392,17 +392,18 @@ QQMPolyRingElem[x1 + 2*x2 + 2*x3 - 1, x1^2 - x1 + 2*x2^2 + 2*x3^2, 2*x1*x2 + 2*x
 
 julia> real_solutions(I)
 4-element Vector{Vector{QQFieldElem}}:
- [5416829397//8589934592, 2708414699//8589934592, -2844258330290649520990905062759917788583//21778071482940061661655974875633165533184]
+ [21667317587//34359738368, 21667317587//68719476736, -2844258330290649520990905062759917788583//21778071482940061661655974875633165533184]
  [1, 0, 0]
- [1945971683//8589934592, 972985841//8589934592, 744426424910260862653434112767010536665//2722258935367507707706996859454145691648]
- [2863311531//8589934592, 0, 3629678580490010276942662479272194255531//10889035741470030830827987437816582766592]
+ [15567773459//68719476736, 15567773459//137438953472, 744426424910260862653434112767010536665//2722258935367507707706996859454145691648]
+ [11453246123//34359738368, 460175067//5846006549323611672814739330865132078623730171904, 3629678580490010276942662479272194255531//10889035741470030830827987437816582766592]
 
 julia> real_solutions(I, interval=true)
 4-element Vector{Vector{Vector{QQFieldElem}}}:
- [[1354207349//2147483648, 2708414699//4294967296], [1354207349//4294967296, 677103675//2147483648], [-355532291286331190123863132844989723573//2722258935367507707706996859454145691648, -1422129165145324760495452531379958894291//10889035741470030830827987437816582766592]]
+ [[10833658793//17179869184, 5416829397//8589934592], [10833658793//34359738368, 5416829397//17179869184], [-355532291286331190123863132844989723573//2722258935367507707706996859454145691648, -1422129165145324760495452531379958894291//10889035741470030830827987437816582766592]]
  [[1, 1], [0, 0], [0, 0]]
- [[972985841//4294967296, 486492921//2147483648], [60811615//536870912, 486492921//4294967296], [93053303113782607831679264095876317083//340282366920938463463374607431768211456, 372213212455130431326717056383505268333//1361129467683753853853498429727072845824]]
- [[1431655765//4294967296, 715827883//2147483648], [-1//4294967296, 1//4294967296], [1814839290245005138471331239636097127765//5444517870735015415413993718908291383296, 907419645122502569235665619818048563883//2722258935367507707706996859454145691648]]
+ [[7783886729//34359738368, 3891943365//17179869184], [7783886729//68719476736, 3891943365//34359738368], [93053303113782607831679264095876317083//340282366920938463463374607431768211456, 372213212455130431326717056383505268333//1361129467683753853853498429727072845824]]
+ [[5726623061//17179869184, 2863311531//8589934592], [-4371663141//1461501637330902918203684832716283019655932542976, 9203501349//2923003274661805836407369665432566039311865085952], [1814839290245005138471331239636097127765//5444517870735015415413993718908291383296, 907419645122502569235665619818048563883//2722258935367507707706996859454145691648]]
+
 ```
 """
 function real_solutions(

src/types.jl

@@ -27,7 +27,7 @@ end
 
 mutable struct Ideal{T <: MPolyRingElem}
     gens::Vector{T}
-    dim::Int
+    dim::Union{Int, Nothing}
     gb::Dict{Int, Vector{T}}
     inter_sols::Vector{Vector{Vector{QQFieldElem}}}
     real_sols::Vector{Vector{QQFieldElem}}
@@ -37,12 +37,11 @@ mutable struct Ideal{T <: MPolyRingElem}
     function Ideal(F::Vector{T}) where {T <: MPolyRingElem}
         I = new{T}()
         I.gens = F
-        I.dim  = -1
         I.gb   = Dict()
-
+        I.dim  = nothing
         return I
     end
-end 
+end
 
 Base.parent(I::Ideal) = Nemo.parent(I.gens[1])
 

test/algorithms/dimension.jl

@@ -2,11 +2,14 @@
     R, (x,y) = polynomial_ring(QQ,["x","y"])
     I =  Ideal([x^2,x*y])
     @test isone(dimension(I))
+    @test isone(I.dim)
 
     R, (x,y,z) = polynomial_ring(GF(101),["x","y","z"])
     I = Ideal([x+2*y+2*z-1, x^2+2*y^2+2*z^2-x, 2*x*y+2*y*z-y])
     @test iszero(dimension(I))
+    @test iszero(I.dim)
 
     I = Ideal([R(1)])
     @test dimension(I) == -1
+    @test I.dim == -1
 end

test/algorithms/groebner-bases.jl

@@ -16,7 +16,7 @@
          z^4 + 38*z^3 + 95*z^2 + 95*z
         ]
     @test G == H
-    
+
     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, eliminate=2, intersect=false)
     H = MPolyRingElem[

test/algorithms/solvers.jl

@@ -5,27 +5,28 @@
          2*x1*x2 + 2*x2*x3 + 2*x3*x4 - x2,
          x2^2 + 2*x1*x3 + 2*x2*x4 - x3])
     sols = Vector{QQFieldElem}[
-        [4862548775//8589934592, 1281562925//8589934592, 2195068207//8589934592, -32714273694608759819673593948262790419337//174224571863520493293247799005065324265472],
-        [3779635503//8589934592, 2638476131//8589934592, 908473689//8589934592, -92633843493479102248442236077889080803361//696898287454081973172991196020261297061888],
+        [9725097551//17179869184, 10252503395//68719476736, 8780272825//34359738368, -32714273694608759819673593948262790419337//174224571863520493293247799005065324265472],
+        [15118542013//34359738368, 10553904527//34359738368, 14535579013//137438953472, -92633843493479102248442236077889080803361//696898287454081973172991196020261297061888],
         [1, 0, 0, 0],
-        [6410479475//8589934592, 2005530653//8589934592, -1585770177//8589934592, 55658687714722201275489356321691134094555601//713623846352979940529142984724747568191373312],
-        [1611414365//8589934592, 673053615//8589934592, 632173751//8589934592, 708759148891639684402860468800417934359477//2787593149816327892691964784081045188247552],
-        [2863311531//8589934592, 0, 0, 14518714321960041107770649917088777022123//43556142965880123323311949751266331066368]
+        [12820958949//17179869184, 16044245225//68719476736, -12686161419//68719476736, 55658687714722201275489356321691134094555601//713623846352979940529142984724747568191373312],
+        [12891314913//68719476736, 10768857835//137438953472, 10114780001//137438953472, 708759148891639684402860468800417934359477//2787593149816327892691964784081045188247552],
+        [11453246123//34359738368, 393751//45671926166590716193865151022383844364247891968, 5717949//182687704666362864775460604089535377456991567872, 14518714321960041107770649917088777022123//43556142965880123323311949751266331066368]
                             ]
     inter_sols = Vector{Vector{QQFieldElem}}[
-        [[2431274387//4294967296, 607818597//1073741824], [320390731//2147483648, 640781463//4294967296], [1097534103//4294967296, 137191763//536870912], [-16357136847304379909836796974131395209669//87112285931760246646623899502532662132736, -4089284211826094977459199243532848802417//21778071482940061661655974875633165533184]],
-        [[1889817751//4294967296, 236227219//536870912], [1319238065//4294967296, 659619033//2147483648], [113559211//1073741824, 454236845//4294967296], [-46316921746739551124221118038944540401681//348449143727040986586495598010130648530944, -2894807609171221945263819877434033775105//21778071482940061661655974875633165533184]],
+        [[4862548775//8589934592, 607818597//1073741824], [5126251697//34359738368, 2563125849//17179869184], [1097534103//4294967296, 4390136413//17179869184], [-16357136847304379909836796974131395209669//87112285931760246646623899502532662132736, -4089284211826094977459199243532848802417//21778071482940061661655974875633165533184]],
+        [[3779635503//8589934592, 7559271007//17179869184], [5276952263//17179869184, 659619033//2147483648], [3633894753//34359738368, 7267789507//68719476736], [-46316921746739551124221118038944540401681//348449143727040986586495598010130648530944, -2894807609171221945263819877434033775105//21778071482940061661655974875633165533184]],
         [[1, 1], [0, 0], [0, 0], [0, 0]],
-        [[3205239737//4294967296, 1602619869//2147483648], [501382663//2147483648, 1002765327//4294967296], [-792885089//4294967296, -24777659//134217728], [3478667982170137579718084770105695880909725//44601490397061246283071436545296723011960832, 27829343857361100637744678160845567047277801//356811923176489970264571492362373784095686656]],
-        [[402853591//2147483648, 805707183//4294967296], [336526807//4294967296, 42065851//536870912], [316086875//4294967296, 79021719//1073741824], [177189787222909921100715117200104483589869//696898287454081973172991196020261297061888, 354379574445819842201430234400208967179739//1393796574908163946345982392040522594123776]],
-        [[1431655765//4294967296, 715827883//2147483648], [-1//4294967296, 1//4294967296], [-1//4294967296, 1//4294967296], [7259357160980020553885324958544388511061//21778071482940061661655974875633165533184, 3629678580490010276942662479272194255531//10889035741470030830827987437816582766592]]
+        [[3205239737//4294967296, 6410479475//8589934592], [2005530653//8589934592, 8022122613//34359738368], [-3171540355//17179869184, -6343080709//34359738368], [3478667982170137579718084770105695880909725//44601490397061246283071436545296723011960832, 27829343857361100637744678160845567047277801//356811923176489970264571492362373784095686656]],
+        [[402853591//2147483648, 6445657457//34359738368], [5384428917//68719476736, 2692214459//34359738368], [316086875//4294967296, 5057390001//68719476736], [177189787222909921100715117200104483589869//696898287454081973172991196020261297061888, 354379574445819842201430234400208967179739//1393796574908163946345982392040522594123776]],
+        [[5726623061//17179869184, 2863311531//8589934592], [-3795019955//365375409332725729550921208179070754913983135744, 3801319971//365375409332725729550921208179070754913983135744], [-2148414955//182687704666362864775460604089535377456991567872, 2159850853//182687704666362864775460604089535377456991567872], [7259357160980020553885324958544388511061//21778071482940061661655974875633165533184, 3629678580490010276942662479272194255531//10889035741470030830827987437816582766592]]
                                             ]
     rat_sols = Vector{QQFieldElem}[[1, 0, 0, 0], [1//3, 0, 0, 1//3]]
 
     @test sols == real_solutions(I)
     @test inter_sols == real_solutions(I, interval=true)
     @test rat_sols == rational_solutions(I)
     @test I.real_sols == real_solutions(I)
+    @test iszero(I.dim)
 
     C, x = polynomial_ring(QQ, "x")
     elim  = 128304*x^8 - 93312*x^7 + 15552*x^6 + 3144*x^5 - 1120*x^4 + 36*x^3 + 15*x^2 - x
@@ -55,11 +56,13 @@
     I = Ideal([x1^2-x2, x1*x3-x4, x2*x4-12, x4^3-x3^2])
     real_solutions(I)
     @test I.rat_param.vars == Symbol[]
+    @test I.dim == -1
 
     I = Ideal([x1^2-x2, x1*x3, x2-12])
 	@test_throws ErrorException real_solutions(I)
     @test_throws ErrorException real_solutions(I, interval=true)
 	@test_throws ErrorException rational_solutions(I)
+    @test isnothing(I.dim)
 
     R, (x, y) = polynomial_ring(QQ,["x","y"])
     # issue 54