From f2d5b93c65bc682114b4f33e31d38db9c25a69a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Fri, 18 Apr 2025 15:33:34 +0200 Subject: [PATCH 01/33] add current files --- src/connectivity/Cannytools.jl | 309 +++++++++++++++++++++++++++++++++ src/connectivity/computeRM.jl | 90 ++++++++++ 2 files changed, 399 insertions(+) create mode 100644 src/connectivity/Cannytools.jl create mode 100644 src/connectivity/computeRM.jl diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl new file mode 100644 index 0000000..27b970b --- /dev/null +++ b/src/connectivity/Cannytools.jl @@ -0,0 +1,309 @@ +# Return the polynomials in F, but injected in the polynomial ring with newvarias_S as new variables +function change_ringvar(F::Vector{P}, newvarias_S::Vector{Symbol}) where {P <: MPolyRingElem} + R = parent(first(F)) + # Locate variables of R in newvarias + to_varias = Vector{Int}(undef,0) + for v in newvarias_S + ind = findfirst(x->x==v, R.S) + push!(to_varias, typeof(ind)==Nothing ? length(R.S)+1 : ind) + end + + ind_novarias = setdiff(eachindex(R.S), to_varias) + newR, newvarias = polynomial_ring(base_ring(R), newvarias_S) + + res = typeof(first(F))[] + ctx = MPolyBuildCtx(newR) + + for f in F + for (e, c) in zip(exponent_vectors(f), coefficients(f)) + @assert(all([ e[i]==0 for i in ind_novarias ]), "Occurence of old variable.s found!") + push!(e, 0) + push_term!(ctx, c, [e[i] for i in to_varias ]) + end + push!(res, finish(ctx)) + end + + return res +end + +function change_ringvar(F::Vector{P}, newvarias_S::Vector{Symbol}) where {P <: PolyRingElem} + R = parent(first(F)) + to_varias = [ v==R.S ? 1 : 2 for v in newvarias_S] + newR, = polynomial_ring(base_ring(R), newvarias_S) + + res = typeof(zero(newR))[] + ctx = MPolyBuildCtx(newR) + + LcF = [ filter!(t->t[2]!=0, collect(enumerate(coefficients(f)))) for f in F ] + for f in LcF + for (ef, c) in f + e = [ef-1, 0] + push_term!(ctx, c, [e[i] for i in to_varias ]) + end + push!(res, finish(ctx)) + end + + return res +end + +function change_ringvar(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}) + return first(change_ringvar([f], newvarias_S)) +end + +function change_ringvar(F::Vector{P}, ind_newvarias::Union{Vector{I}, UnitRange{I}}) where {P <: MPolyRingElem, I <: Int64} + R = parent(first(F)) + return change_ringvar(F, [R.S[i] for i in ind_newvarias]) +end + +function change_ringvar(f::MPolyRingElem, ind_newvarias::Union{Vector{I}, UnitRange{I}}) where {I <: Int64} + R = parent(f) + return first(change_ringvar([f], [R.S[i] for i in ind_newvarias])) +end + +# Return the polynomials in F, but injected in the polynomial ring with the variables occuring in F +function change_ringvar(F::Vector{P}) where {P <: MPolyRingElem} + union_varias = Set{Symbol}() + for f in F + union!(union_varias, map(Symbol, vars(f)) ) + end + return change_ringvar(F, collect(union_varias)) +end + +function change_ringvar(f::MPolyRingElem) + return first(change_ringvar([f])) +end + +# Return the polynomials in F, but injected in the polynomial ring with newvarias_S as new variables +function change_ringvar_mod(F::Vector{P}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) where {P <: MPolyRingElem} + R = parent(first(F)) + # Locate variables of R in newvarias + to_varias = Vector{Int}(undef,0) + for v in newvarias_S + ind = findfirst(x->x==v, oldvarias_S) + push!(to_varias, typeof(ind)==Nothing ? length(oldvarias_S)+1 : ind) + end + + ind_novarias = setdiff(eachindex(oldvarias_S), to_varias) + newR, newvarias = polynomial_ring(base_ring(R), newvarias_S) + + res = typeof(first(F))[] + ctx = MPolyBuildCtx(newR) + + for f in F + for (e, c) in zip(exponent_vectors(f), coefficients(f)) + @assert(all([ e[i]==0 for i in ind_novarias ]), "Occurence of old variable.s found!") + push!(e, 0) + push_term!(ctx, c, [e[i] for i in to_varias ]) + end + push!(res, finish(ctx)) + end + + return res +end + +function change_ringvar_mod(F::Vector{P}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) where {P <: PolyRingElem} + R = parent(first(F)) + A, (x,) = polynomial_ring(base_ring(R), oldvarias_S) + return change_ringvar_mod([ evaluate(f, x) for f in F ], newvarias_S, oldvarias_S) +end + +function change_ringvar_mod(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) + return first(change_ringvar_mod([f], newvarias_S, oldvarias_S)) +end + +#function change_ringvar_mod(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}) +# return first(change_ringvar_mod([f], newvarias_S, oldvarias_S)) +#end + +function MPolyBuild(F::Vector{Vector{P}}, newvarias_S::Vector{Symbol}, idx::Int) where {P <: RingElem} + A = parent(first(first(F))) + R, = polynomial_ring(A, newvarias_S) + to_varias = [ i==idx ? 1 : 2 for i in eachindex(newvarias_S) ] + + res = typeof(zero(R))[] + ctx = MPolyBuildCtx(R) + + LcF = [ filter!(t->t[2]!=0, collect(enumerate(f))) for f in F ] + for f in LcF + for (ef, c) in f + e = [ef-1, 0] + push_term!(ctx, c, [e[i] for i in to_varias ]) + end + push!(res, finish(ctx)) + end + + return res +end + +function MPolyBuild(f::Vector{P}, newvarias_S::Vector{Symbol}, idx::Int) where {P <: RingElem} + return first(MPolyBuild([f], newvarias_S, idx)) +end + +function computepolarproj(j::Int, V::AlgebraicSolving.Ideal, dimV::Int, varbs; dimproj=j-1, characteristic=0, output="minors", verb=0, nr_thrds=1) + #println("$j $V $dimV") + # Compute the set of points x where pi_j(T_x(V)) has dimension < dimproj + @assert output in ["minors", "groebner", "real", "parametric", "interval"] "Wrong output parameter" + R = parent(V) + n, hV = R.nvars, V.gens + c = n - dimV + + JW = transpose([ derivative(f, k) for k=j+1:n, f in hV ]) + sizeminors = c + dimproj - (j-1) + hW = vcat(hV, compute_minors(sizeminors, JW, R)) + filter!(!iszero, hW) + output_functions = Dict( + "minors" => x -> x, + "groebner" => x -> groebner_basis(x, info_level=verb,nr_thrds=nr_thrds), + "real" => x -> real_solutions(x, info_level=verb,nr_thrds=nr_thrds), + "interval" => x -> real_solutions(x, interval=true, info_level=verb,nr_thrds=nr_thrds), + "parametric" => x -> rational_parametrization(x, info_level=verb,nr_thrds=nr_thrds) + ) + + return output_functions[output](AlgebraicSolving.Ideal(hW)) +end + +function computepolarphi(j::Int, V::AlgebraicSolving.Ideal, phi::Vector{T} where T<:MPolyRingElem, dimV::Int, varbs; dimproj=j-1, characteristic=0, output="minors", verb=0, nr_thrds=1) + # Compute the set of points x where pi_j(T_x(V)) has dimension < dimproj + @assert output in ["minors", "groebner", "real", "parametric", "interval"] "Wrong output parameter" + R = parent(V) + n, hV = nvars(R), V.gens + c = n - dimV + + if j==0 + JW = transpose([ derivative(f, k) for k=1:n, f in hV ]) + sizeminors = c + else + JW = transpose([ derivative(f, k) for k=1:n, f in vcat(hV,phi) ]) + sizeminors = c + j + dimproj - (j-1) + end + hW = vcat(hV, compute_minors(sizeminors, JW, R)) + output_functions = Dict( + "minors" => x -> x, + "groebner" => x -> groebner_basis(x, info_level=verb,nr_thrds=nr_thrds), + "real" => x -> real_solutions(x, info_level=verb,nr_thrds=nr_thrds), + "interval" => x -> real_solutions(x, interval=true, info_level=verb,nr_thrds=nr_thrds), + "parametric" => x -> rational_parametrization(x, info_level=verb,nr_thrds=nr_thrds) + ) + + return output_functions[output](AlgebraicSolving.Ideal(hW)) +end + +function compute_minors(p, A, R) + #Computes the p-minors of a matrix A + n, m = size(A) + rowsmins = collect(combinations(1:n, p)) + colsmins = collect(combinations(1:m, p)) + mins = Vector{eltype(A)}(undef, length(rowsmins) * length(colsmins)) + k = 1 + for rowsmin in rowsmins + for colsmin in colsmins + mins[k] = detmpoly(A[rowsmin, colsmin], R) + k += 1 + end + end + + return mins +end + +function combinations(a, n) + # Helper function to recursively generate combinations + function _combinations(a, n, start, chosen) + if length(chosen) == n + return [chosen] + elseif start > length(a) + return Vector{Int}([]) + else + # Include the current element and recurse + include_current = _combinations(a, n, start + 1, [chosen; a[start]]) + # Exclude the current element and recurse + exclude_current = _combinations(a, n, start + 1, chosen) + return vcat(include_current, exclude_current) + end + end + return _combinations(a, n, 1, Vector{Int}([])) +end + +function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) + # Get the size of the matrix + n = size(A, 1) + if n != size(A, 2) + throw(ArgumentError("Matrix must be square")) + end + + if n == 1 + return A[1, 1] + end + + # Initialize the determinant polynomial + detA = zero(R) + + # Compute the determinant polynomial + for j = 1:n + submatrix = A[2:end, [i for i = 1:n if i != j]] + detA += (-1)^(1+j)*A[1, j] * detmpoly(submatrix, R) + end + + return detA +end + +function contfrac_convergents(x::Rational{Int}) + q = Rational{Int}[] + fpart, ipart = modf(x) + push!(q, ipart) + na, da, nb, db = Int(1), Int(0), Int(ipart), Int(1) + while true + fpart == 0 && return q + x = inv(fpart) + fpart, ipart = modf(x) + na, da, nb, db = nb, db, na+ipart*nb, da+ipart*db + push!(q, nb//db) + end +end + +function modfQQ(x::QQFieldElem) + ipart, fpart = [ f(numerator(x), denominator(x)) for f in [div, rem]] + return QQFieldElem(fpart, denominator(x)), ipart +end + +function small_mid_point(a::QQFieldElem,b::QQFieldElem) + a == b && return a + if b < a + a, b = b, a + end + x = (a+b)/ QQFieldElem(2) + fpart, ipart = modfQQ(x) + println("\n",fpart," , ", ipart) + q1, q2 = QQFieldElem(1//0), QQFieldElem(ipart,1) + while true + #println("(",q1," ; ",q2,")") + if numerator(fpart) == 0 || (q2 > a && q2 < b) + return q2 + end + x = inv(fpart) + fpart, ipart = modfQQ(x) + println(fpart," , ", ipart) + q1, q2 = q2, QQFieldElem(numerator(q1) + ipart*numerator(q2), denominator(q1) + ipart*denominator(q2)) + end +end + +function MidRationalPoints(S::Vector{Vector{QQFieldElem}}) + # S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] + # such that the [l_i, r_i] are rational and disjoint open intervals. + # + # It orders the [l_i,r_i], removes repetitions and computes the simplest + # rational points between these intervals. + if isempty(S) + return QQFieldElem[] + end + S1 = sort(S, lt=(x, y) -> x[2] <= y[1]) + n = unique!(S1) |> length + ratioP = Vector{QQFieldElem}(undef, n-1) + for i in 1:(n- 1) + ri, li1 = S1[i][2], S1[i+1][1] + @assert ri < li1 "The intervals are not disjoint." + #ratioP[i] = small_mid_point(ri, li1) + ratioP[i] = simplest_between(ri, li1) + end + return ratioP +end + diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl new file mode 100644 index 0000000..504fc95 --- /dev/null +++ b/src/connectivity/computeRM.jl @@ -0,0 +1,90 @@ +#using Pkg +#Pkg.activate("AlgebraicSolving.jl") +#using Revise +#using AlgebraicSolving +#using Nemo + +export computeRM, computepolarproj, computepolarphi +# DEBUG +export change_ringvar, compute_minors, detmpoly, change_ringvar_mod, MPolyBuild + +include("Cannytools.jl") + +function computeRM(V::Ideal{T} where T <: QQMPolyRingElem, dimV::Int; Q=Vector{Vector{QQFieldElem}}([[]]), C=Vector{Vector{QQFieldElem}}([]), v=0) + # V is a variety of dimension dimV + # Q are base points with rational coefficients + # C are query points with rational coefficients + println(Q) + A = parent(V) + varias, nvarias = gens(A), nvars(A) + hV = V.gens + #println(Q) + e = length(first(Q)) + fixvarias, newvarias = varias[1:e], varias[e+1:nvarias] + + R = Vector{AlgebraicSolving.Ideal{QQMPolyRingElem}}(undef,0) + + for q in Q + ## Fq ## + if e > 0 + hFq = change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], e+1:nvarias) + Fq = AlgebraicSolving.Ideal(hFq) + else + Fq = V + end + + if dimV - e <= 1 + curve = change_ringvar(Fq.gens, A.S) + push!(R, AlgebraicSolving.Ideal(vcat(curve, [fixvarias[j] - q[j] for j in 1:e]))) + else + ## sing(Fq) ## + v>0 && println("Compute first the singular points") + singFq = computepolarproj(0, Fq, dimV-e, newvarias, output="real", verb=max(v-1,0), nr_thrds=Threads.nthreads()) + @assert(length(singFq)==0, "Non-emtpy real sing locus!") + + ## K(pi_1,Fq) ## + v>0 && println("First critical points") + K1Fq = computepolarproj(1, Fq, dimV-e, newvarias, output="interval", verb=max(v-1,0), nr_thrds=Threads.nthreads()) + + ## K(pi_2, Fq) ## + v>0 && println("Second critical points") + K2Fq = computepolarproj(2, Fq, dimV-e, newvarias, output="minors", verb=max(v-1,0), nr_thrds=Threads.nthreads()) + polar = change_ringvar(K2Fq.gens, A.S) + push!(R, AlgebraicSolving.Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) + + ## Points with vertical tg in K(pi_2, Fq) ## + v>0 && println("Vertical tg points") + K1WmFq = computepolarproj(2, K2Fq, 1, newvarias, output="interval", dimproj=0, verb=max(v-1,0), nr_thrds=Threads.nthreads()) + + ## New base points ## + K1W = vcat(K1Fq, K1WmFq) + K1WRat = MidRationalPoints(getindex(K1W, 1)) + # Heuristic to be proven + #K1WRat = K1WRat[2:end-1] + ########## + append!(K1WRat, [ c[e+1] for c in C if c[1:e]==q && !(c[e+1] in K1WRat) ]) |> sort! + newQ = [ vcat(q, [kv]) for kv in K1WRat ] + + if !isempty(newQ) + RFq = computeRM(V, dimV, Q=newQ, C=C) + R = vcat(R, RFq) + end + end + end + + return R +end + +#= +## Test ## +R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) +#h = x1^2+x2^2+x3^2+x4^2-1 +h = (x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) +h = evaluate(h,[x1,x2,x3],[x1+rand(-10:10)*x2+rand(-10:10)*x3+rand(-10:10)*x4,x2+rand(-10:10)*x3+rand(-10:10)*x4,x3+rand(-10:10)*x4]) +V = AlgebraicSolving.Ideal([h]) + +carte = (computeRM(V, 3, [Vector{QQFieldElem}(undef,0)])) +#println(carte) +=# + + From 85d734c485cde27f0f1a481d3d895bcf5b992eef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Fri, 18 Apr 2025 17:18:53 +0200 Subject: [PATCH 02/33] simplify computepolar + remove useless functions + remove useless checks --- src/AlgebraicSolving.jl | 2 + src/connectivity/Cannytools.jl | 215 +++------------------------------ src/connectivity/computeRM.jl | 142 ++++++++++++---------- 3 files changed, 96 insertions(+), 263 deletions(-) diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index 5cab42f..9c66519 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -18,6 +18,8 @@ include("algorithms/paramcurve.jl") include("algorithms/hilbert.jl") #= siggb =# include("siggb/siggb.jl") +#= connectivity =# +include("connectivity/computeRM.jl") #= examples =# include("examples/katsura.jl") include("examples/cyclic.jl") diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 27b970b..16c274a 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -1,191 +1,22 @@ -# Return the polynomials in F, but injected in the polynomial ring with newvarias_S as new variables -function change_ringvar(F::Vector{P}, newvarias_S::Vector{Symbol}) where {P <: MPolyRingElem} - R = parent(first(F)) - # Locate variables of R in newvarias - to_varias = Vector{Int}(undef,0) - for v in newvarias_S - ind = findfirst(x->x==v, R.S) - push!(to_varias, typeof(ind)==Nothing ? length(R.S)+1 : ind) - end - - ind_novarias = setdiff(eachindex(R.S), to_varias) - newR, newvarias = polynomial_ring(base_ring(R), newvarias_S) - - res = typeof(first(F))[] - ctx = MPolyBuildCtx(newR) - - for f in F - for (e, c) in zip(exponent_vectors(f), coefficients(f)) - @assert(all([ e[i]==0 for i in ind_novarias ]), "Occurence of old variable.s found!") - push!(e, 0) - push_term!(ctx, c, [e[i] for i in to_varias ]) - end - push!(res, finish(ctx)) - end - - return res -end - -function change_ringvar(F::Vector{P}, newvarias_S::Vector{Symbol}) where {P <: PolyRingElem} - R = parent(first(F)) - to_varias = [ v==R.S ? 1 : 2 for v in newvarias_S] - newR, = polynomial_ring(base_ring(R), newvarias_S) - - res = typeof(zero(newR))[] - ctx = MPolyBuildCtx(newR) - - LcF = [ filter!(t->t[2]!=0, collect(enumerate(coefficients(f)))) for f in F ] - for f in LcF - for (ef, c) in f - e = [ef-1, 0] - push_term!(ctx, c, [e[i] for i in to_varias ]) - end - push!(res, finish(ctx)) - end - - return res -end - -function change_ringvar(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}) - return first(change_ringvar([f], newvarias_S)) -end - -function change_ringvar(F::Vector{P}, ind_newvarias::Union{Vector{I}, UnitRange{I}}) where {P <: MPolyRingElem, I <: Int64} - R = parent(first(F)) - return change_ringvar(F, [R.S[i] for i in ind_newvarias]) -end - -function change_ringvar(f::MPolyRingElem, ind_newvarias::Union{Vector{I}, UnitRange{I}}) where {I <: Int64} - R = parent(f) - return first(change_ringvar([f], [R.S[i] for i in ind_newvarias])) -end - -# Return the polynomials in F, but injected in the polynomial ring with the variables occuring in F -function change_ringvar(F::Vector{P}) where {P <: MPolyRingElem} - union_varias = Set{Symbol}() - for f in F - union!(union_varias, map(Symbol, vars(f)) ) - end - return change_ringvar(F, collect(union_varias)) -end - -function change_ringvar(f::MPolyRingElem) - return first(change_ringvar([f])) -end - -# Return the polynomials in F, but injected in the polynomial ring with newvarias_S as new variables -function change_ringvar_mod(F::Vector{P}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) where {P <: MPolyRingElem} - R = parent(first(F)) - # Locate variables of R in newvarias - to_varias = Vector{Int}(undef,0) - for v in newvarias_S - ind = findfirst(x->x==v, oldvarias_S) - push!(to_varias, typeof(ind)==Nothing ? length(oldvarias_S)+1 : ind) - end - - ind_novarias = setdiff(eachindex(oldvarias_S), to_varias) - newR, newvarias = polynomial_ring(base_ring(R), newvarias_S) - - res = typeof(first(F))[] - ctx = MPolyBuildCtx(newR) - - for f in F - for (e, c) in zip(exponent_vectors(f), coefficients(f)) - @assert(all([ e[i]==0 for i in ind_novarias ]), "Occurence of old variable.s found!") - push!(e, 0) - push_term!(ctx, c, [e[i] for i in to_varias ]) - end - push!(res, finish(ctx)) - end - - return res -end - -function change_ringvar_mod(F::Vector{P}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) where {P <: PolyRingElem} - R = parent(first(F)) - A, (x,) = polynomial_ring(base_ring(R), oldvarias_S) - return change_ringvar_mod([ evaluate(f, x) for f in F ], newvarias_S, oldvarias_S) -end - -function change_ringvar_mod(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}, oldvarias_S::Vector{Symbol}) - return first(change_ringvar_mod([f], newvarias_S, oldvarias_S)) -end - -#function change_ringvar_mod(f::Union{MPolyRingElem, PolyRingElem}, newvarias_S::Vector{Symbol}) -# return first(change_ringvar_mod([f], newvarias_S, oldvarias_S)) -#end - -function MPolyBuild(F::Vector{Vector{P}}, newvarias_S::Vector{Symbol}, idx::Int) where {P <: RingElem} - A = parent(first(first(F))) - R, = polynomial_ring(A, newvarias_S) - to_varias = [ i==idx ? 1 : 2 for i in eachindex(newvarias_S) ] - - res = typeof(zero(R))[] - ctx = MPolyBuildCtx(R) - - LcF = [ filter!(t->t[2]!=0, collect(enumerate(f))) for f in F ] - for f in LcF - for (ef, c) in f - e = [ef-1, 0] - push_term!(ctx, c, [e[i] for i in to_varias ]) - end - push!(res, finish(ctx)) - end - - return res -end - -function MPolyBuild(f::Vector{P}, newvarias_S::Vector{Symbol}, idx::Int) where {P <: RingElem} - return first(MPolyBuild([f], newvarias_S, idx)) -end - -function computepolarproj(j::Int, V::AlgebraicSolving.Ideal, dimV::Int, varbs; dimproj=j-1, characteristic=0, output="minors", verb=0, nr_thrds=1) - #println("$j $V $dimV") - # Compute the set of points x where pi_j(T_x(V)) has dimension < dimproj - @assert output in ["minors", "groebner", "real", "parametric", "interval"] "Wrong output parameter" +# Compute the set of points x where phi_j(T_x(V)) has dimension < dimproj +function computepolar( + j::Int, # j-th first coordinate images of phi + V::Ideal{P}; # input ideal + phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) + dimproj = j-1, # maximum dimension of tangent space of phi + v=0, # verbosity level + ) where (P <: MPolyRingElem) + V.dim == -1 && dimension(V) R = parent(V) - n, hV = R.nvars, V.gens - c = n - dimV + n = nvars(R) + c = n - V.dim + nphi = length(phi) - JW = transpose([ derivative(f, k) for k=j+1:n, f in hV ]) - sizeminors = c + dimproj - (j-1) - hW = vcat(hV, compute_minors(sizeminors, JW, R)) - filter!(!iszero, hW) - output_functions = Dict( - "minors" => x -> x, - "groebner" => x -> groebner_basis(x, info_level=verb,nr_thrds=nr_thrds), - "real" => x -> real_solutions(x, info_level=verb,nr_thrds=nr_thrds), - "interval" => x -> real_solutions(x, interval=true, info_level=verb,nr_thrds=nr_thrds), - "parametric" => x -> rational_parametrization(x, info_level=verb,nr_thrds=nr_thrds) - ) + JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in vcat(V.gens, phi[1:min(j,nphi)])]) + sizeminors = c + min(nphi,j) + min(dimproj,j-1) - (j-1) + minors = compute_minors(sizeminors, JW, R) - return output_functions[output](AlgebraicSolving.Ideal(hW)) -end - -function computepolarphi(j::Int, V::AlgebraicSolving.Ideal, phi::Vector{T} where T<:MPolyRingElem, dimV::Int, varbs; dimproj=j-1, characteristic=0, output="minors", verb=0, nr_thrds=1) - # Compute the set of points x where pi_j(T_x(V)) has dimension < dimproj - @assert output in ["minors", "groebner", "real", "parametric", "interval"] "Wrong output parameter" - R = parent(V) - n, hV = nvars(R), V.gens - c = n - dimV - - if j==0 - JW = transpose([ derivative(f, k) for k=1:n, f in hV ]) - sizeminors = c - else - JW = transpose([ derivative(f, k) for k=1:n, f in vcat(hV,phi) ]) - sizeminors = c + j + dimproj - (j-1) - end - hW = vcat(hV, compute_minors(sizeminors, JW, R)) - output_functions = Dict( - "minors" => x -> x, - "groebner" => x -> groebner_basis(x, info_level=verb,nr_thrds=nr_thrds), - "real" => x -> real_solutions(x, info_level=verb,nr_thrds=nr_thrds), - "interval" => x -> real_solutions(x, interval=true, info_level=verb,nr_thrds=nr_thrds), - "parametric" => x -> rational_parametrization(x, info_level=verb,nr_thrds=nr_thrds) - ) - - return output_functions[output](AlgebraicSolving.Ideal(hW)) + return Ideal(vcat(V.gens, minors)) end function compute_minors(p, A, R) @@ -246,20 +77,6 @@ function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) return detA end -function contfrac_convergents(x::Rational{Int}) - q = Rational{Int}[] - fpart, ipart = modf(x) - push!(q, ipart) - na, da, nb, db = Int(1), Int(0), Int(ipart), Int(1) - while true - fpart == 0 && return q - x = inv(fpart) - fpart, ipart = modf(x) - na, da, nb, db = nb, db, na+ipart*nb, da+ipart*db - push!(q, nb//db) - end -end - function modfQQ(x::QQFieldElem) ipart, fpart = [ f(numerator(x), denominator(x)) for f in [div, rem]] return QQFieldElem(fpart, denominator(x)), ipart diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 504fc95..418a4fb 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,75 +4,89 @@ #using AlgebraicSolving #using Nemo -export computeRM, computepolarproj, computepolarphi +export computeRM, computepolarproj, computepolarphi, computepolar # DEBUG -export change_ringvar, compute_minors, detmpoly, change_ringvar_mod, MPolyBuild +export compute_minors, detmpoly, change_ringvar_mod, MPolyBuild include("Cannytools.jl") -function computeRM(V::Ideal{T} where T <: QQMPolyRingElem, dimV::Int; Q=Vector{Vector{QQFieldElem}}([[]]), C=Vector{Vector{QQFieldElem}}([]), v=0) - # V is a variety of dimension dimV - # Q are base points with rational coefficients - # C are query points with rational coefficients - println(Q) - A = parent(V) - varias, nvarias = gens(A), nvars(A) - hV = V.gens - #println(Q) - e = length(first(Q)) - fixvarias, newvarias = varias[1:e], varias[e+1:nvarias] - - R = Vector{AlgebraicSolving.Ideal{QQMPolyRingElem}}(undef,0) - - for q in Q - ## Fq ## - if e > 0 - hFq = change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], e+1:nvarias) - Fq = AlgebraicSolving.Ideal(hFq) - else - Fq = V - end - - if dimV - e <= 1 - curve = change_ringvar(Fq.gens, A.S) - push!(R, AlgebraicSolving.Ideal(vcat(curve, [fixvarias[j] - q[j] for j in 1:e]))) - else - ## sing(Fq) ## - v>0 && println("Compute first the singular points") - singFq = computepolarproj(0, Fq, dimV-e, newvarias, output="real", verb=max(v-1,0), nr_thrds=Threads.nthreads()) - @assert(length(singFq)==0, "Non-emtpy real sing locus!") - - ## K(pi_1,Fq) ## - v>0 && println("First critical points") - K1Fq = computepolarproj(1, Fq, dimV-e, newvarias, output="interval", verb=max(v-1,0), nr_thrds=Threads.nthreads()) - - ## K(pi_2, Fq) ## - v>0 && println("Second critical points") - K2Fq = computepolarproj(2, Fq, dimV-e, newvarias, output="minors", verb=max(v-1,0), nr_thrds=Threads.nthreads()) - polar = change_ringvar(K2Fq.gens, A.S) - push!(R, AlgebraicSolving.Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) - - ## Points with vertical tg in K(pi_2, Fq) ## - v>0 && println("Vertical tg points") - K1WmFq = computepolarproj(2, K2Fq, 1, newvarias, output="interval", dimproj=0, verb=max(v-1,0), nr_thrds=Threads.nthreads()) - - ## New base points ## - K1W = vcat(K1Fq, K1WmFq) - K1WRat = MidRationalPoints(getindex(K1W, 1)) - # Heuristic to be proven - #K1WRat = K1WRat[2:end-1] - ########## - append!(K1WRat, [ c[e+1] for c in C if c[1:e]==q && !(c[e+1] in K1WRat) ]) |> sort! - newQ = [ vcat(q, [kv]) for kv in K1WRat ] - - if !isempty(newQ) - RFq = computeRM(V, dimV, Q=newQ, C=C) - R = vcat(R, RFq) +function computeRM( + V::Ideal{T} where T <: QQMPolyRingElem; # input ideal + Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients + C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients + v::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) + ) + # Some preprocessing + V.dim == -1 && dimension(V) + println(Q) + isempty(Q) && push!(Q,[]) + # Get data + A = parent(V) + varias, nvarias = gens(A), nvars(A) + e = length(first(Q)) + fixvarias = varias[1:e] + + RM = Vector{Ideal{QQMPolyRingElem}}(undef,0) + for q in Q + ## Fq ## + if e > 0 + hFq = change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], A.S[e+1:end]) + Fq = Ideal(hFq) + else + Fq = V + end + # Genericity assumption (can be checked) + Fq.dim = V.dim - e + + if V.dim - e <= 1 + curve = change_ringvar(Fq.gens, A.S) + push!(RM, Ideal(vcat(curve, [fixvarias[j] - q[j] for j in 1:e]))) + else + ## sing(Fq) ## + if checks + v>0 && println("Compute first the singular points") + singFq = computepolar(0, Fq, v=max(v-1,0)) + @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), "Non-emtpy real sing locus!") + end + + ## K(pi_1,Fq) ## + v>0 && println("First critical points") + K1Fq = computepolar(1, Fq, v=max(v-1,0)) + K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## K(pi_2, Fq) ## + v>0 && println("Second critical points") + K2Fq = computepolar(2, Fq, v=max(v-1,0)) + if checks + @assert(isone(dimension(K2Fq)), "Non-generic polar variety") + end + K2Fq.dim = 1 + polar = change_ringvar(K2Fq.gens, A.S) + push!(RM, Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) + + ## Points with vertical tg in K(pi_2, Fq) ## + v>0 && println("Vertical tg points") + K1WmFq = computepolar(2, K2Fq, dimproj=0, v=max(v-1,0)) + K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## New base points ## + K1W = vcat(K1Fq, K1WmFq) + K1WRat = MidRationalPoints(getindex(K1W, 1)) + # Heuristic to be proven + #K1WRat = K1WRat[2:end-1] + ########## + append!(K1WRat, [ c[e+1] for c in C if c[1:e]==q && !(c[e+1] in K1WRat) ]) |> sort! + newQ = [ vcat(q, [kv]) for kv in K1WRat ] + + if !isempty(newQ) + RMFq = computeRM(V, Q=newQ, C=C) + append!(RM, RMFq) + end end - end - end + end - return R + return RM end #= @@ -81,7 +95,7 @@ R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) #h = x1^2+x2^2+x3^2+x4^2-1 h = (x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) h = evaluate(h,[x1,x2,x3],[x1+rand(-10:10)*x2+rand(-10:10)*x3+rand(-10:10)*x4,x2+rand(-10:10)*x3+rand(-10:10)*x4,x3+rand(-10:10)*x4]) -V = AlgebraicSolving.Ideal([h]) +V = Ideal([h]) carte = (computeRM(V, 3, [Vector{QQFieldElem}(undef,0)])) #println(carte) From 7f0c5614e40be5760791d989a8fbcb2ff9832425 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Fri, 18 Apr 2025 18:25:17 +0200 Subject: [PATCH 03/33] remove query points --- src/connectivity/computeRM.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 418a4fb..1d8d9f8 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -17,13 +17,19 @@ function computeRM( v::Int=0, # verbosity level checks::Bool=false # perform checks (dimension, regularity, etc.) ) + # Some base cases + A = parent(V) + varias = gens(A) + if length(varias)<=2 + return [V] + end # Some preprocessing V.dim == -1 && dimension(V) println(Q) + println(C) + println() isempty(Q) && push!(Q,[]) - # Get data - A = parent(V) - varias, nvarias = gens(A), nvars(A) + # Base points e = length(first(Q)) fixvarias = varias[1:e] @@ -76,11 +82,12 @@ function computeRM( # Heuristic to be proven #K1WRat = K1WRat[2:end-1] ########## - append!(K1WRat, [ c[e+1] for c in C if c[1:e]==q && !(c[e+1] in K1WRat) ]) |> sort! + Cq = [c for c in C if c[1:e]==q] + append!(K1WRat, getindex.(Cq, e+1)) |> sort! newQ = [ vcat(q, [kv]) for kv in K1WRat ] if !isempty(newQ) - RMFq = computeRM(V, Q=newQ, C=C) + RMFq = computeRM(V, Q=newQ, C=Cq) append!(RM, RMFq) end end From 9477d47954c0b1a1cf5f219325cd5c3c3c6ebb82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Tue, 22 Apr 2025 20:18:42 +0200 Subject: [PATCH 04/33] remove custom small_mid_point + debut canny param --- src/connectivity/Cannytools.jl | 30 +-------- src/connectivity/computeRM.jl | 113 ++++++++++++++++++--------------- 2 files changed, 64 insertions(+), 79 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 16c274a..49dd3c6 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -77,32 +77,6 @@ function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) return detA end -function modfQQ(x::QQFieldElem) - ipart, fpart = [ f(numerator(x), denominator(x)) for f in [div, rem]] - return QQFieldElem(fpart, denominator(x)), ipart -end - -function small_mid_point(a::QQFieldElem,b::QQFieldElem) - a == b && return a - if b < a - a, b = b, a - end - x = (a+b)/ QQFieldElem(2) - fpart, ipart = modfQQ(x) - println("\n",fpart," , ", ipart) - q1, q2 = QQFieldElem(1//0), QQFieldElem(ipart,1) - while true - #println("(",q1," ; ",q2,")") - if numerator(fpart) == 0 || (q2 > a && q2 < b) - return q2 - end - x = inv(fpart) - fpart, ipart = modfQQ(x) - println(fpart," , ", ipart) - q1, q2 = q2, QQFieldElem(numerator(q1) + ipart*numerator(q2), denominator(q1) + ipart*denominator(q2)) - end -end - function MidRationalPoints(S::Vector{Vector{QQFieldElem}}) # S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] # such that the [l_i, r_i] are rational and disjoint open intervals. @@ -118,8 +92,8 @@ function MidRationalPoints(S::Vector{Vector{QQFieldElem}}) for i in 1:(n- 1) ri, li1 = S1[i][2], S1[i+1][1] @assert ri < li1 "The intervals are not disjoint." - #ratioP[i] = small_mid_point(ri, li1) - ratioP[i] = simplest_between(ri, li1) + eps = (li1-ri)//1000 + ratioP[i] = simplest_between(ri + eps, li1 - eps) end return ratioP end diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 1d8d9f8..9257763 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,13 +4,22 @@ #using AlgebraicSolving #using Nemo -export computeRM, computepolarproj, computepolarphi, computepolar +export roadmap, computepolarproj, computepolarphi, computepolar # DEBUG -export compute_minors, detmpoly, change_ringvar_mod, MPolyBuild - +export compute_minors, detmpoly, change_ringvar_mod, MPolyBuild, small_mid_point include("Cannytools.jl") -function computeRM( +function roadmap( + V::Ideal{P}; # input ideal + Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients + C::Vector{Vector{P}}=Vector{Vector{P}}(), # query points with rational coefficients + v::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + CQ = real_solutions(C, info_level=max(v-1,0), nr_thrds=Threads.nthreads()) +end + +function roadmap( V::Ideal{T} where T <: QQMPolyRingElem; # input ideal Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients @@ -36,60 +45,62 @@ function computeRM( RM = Vector{Ideal{QQMPolyRingElem}}(undef,0) for q in Q ## Fq ## - if e > 0 - hFq = change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], A.S[e+1:end]) - Fq = Ideal(hFq) + Fq = e > 0 ? Ideal(change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], A.S[e+1:end])) : V + # Genericity assumption (can be checked) + if checks + @assert(dimension(Fq) == V.dim - e, "Non-generic polar variety") else - Fq = V + Fq.dim = V.dim - e end - # Genericity assumption (can be checked) - Fq.dim = V.dim - e - if V.dim - e <= 1 + # Terminal case (dim <=1) + if Fq.dim <= 1 curve = change_ringvar(Fq.gens, A.S) push!(RM, Ideal(vcat(curve, [fixvarias[j] - q[j] for j in 1:e]))) + continue + end + + ## sing(Fq) ## + if checks + v>0 && println("Compute first the singular points") + singFq = computepolar(0, Fq, v=max(v-1,0)) + @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), "Non-emtpy real sing locus!") + end + + ## K(pi_1,Fq) ## + v>0 && println("First critical points") + K1Fq = computepolar(1, Fq, v=max(v-1,0)) + K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## K(pi_2, Fq) ## + v>0 && println("Second critical points") + K2Fq = computepolar(2, Fq, v=max(v-1,0)) + if checks + @assert(isone(dimension(K2Fq)), "Non-generic polar variety") else - ## sing(Fq) ## - if checks - v>0 && println("Compute first the singular points") - singFq = computepolar(0, Fq, v=max(v-1,0)) - @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), "Non-emtpy real sing locus!") - end - - ## K(pi_1,Fq) ## - v>0 && println("First critical points") - K1Fq = computepolar(1, Fq, v=max(v-1,0)) - K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) - - ## K(pi_2, Fq) ## - v>0 && println("Second critical points") - K2Fq = computepolar(2, Fq, v=max(v-1,0)) - if checks - @assert(isone(dimension(K2Fq)), "Non-generic polar variety") - end K2Fq.dim = 1 - polar = change_ringvar(K2Fq.gens, A.S) - push!(RM, Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) - - ## Points with vertical tg in K(pi_2, Fq) ## - v>0 && println("Vertical tg points") - K1WmFq = computepolar(2, K2Fq, dimproj=0, v=max(v-1,0)) - K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) - - ## New base points ## - K1W = vcat(K1Fq, K1WmFq) - K1WRat = MidRationalPoints(getindex(K1W, 1)) - # Heuristic to be proven - #K1WRat = K1WRat[2:end-1] - ########## - Cq = [c for c in C if c[1:e]==q] - append!(K1WRat, getindex.(Cq, e+1)) |> sort! - newQ = [ vcat(q, [kv]) for kv in K1WRat ] - - if !isempty(newQ) - RMFq = computeRM(V, Q=newQ, C=Cq) - append!(RM, RMFq) - end + end + polar = change_ringvar(K2Fq.gens, A.S) + push!(RM, Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) + + ## Points with vertical tg in K(pi_2, Fq) ## + v>0 && println("Vertical tg points") + K1WmFq = computepolar(2, K2Fq, dimproj=0, v=max(v-1,0)) + K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## New base points ## + K1W = vcat(K1Fq, K1WmFq) + K1WRat = MidRationalPoints(getindex(K1W, 1)) + # Heuristic to be proven + #K1WRat = K1WRat[2:end-1] + ########## + Cq = [c for c in C if c[1:e]==q] + append!(K1WRat, getindex.(Cq, e+1)) |> sort! + newQ = [ vcat(q, [kv]) for kv in K1WRat ] + + if !isempty(newQ) + RMFq = roadmap(V, Q=newQ, C=Cq) + append!(RM, RMFq) end end From 0d811458d37af7b337415f7c8b6217f37d404f6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 23 Apr 2025 18:33:23 +0200 Subject: [PATCH 05/33] roadmap query ideal --- src/algorithms/paramcurve.jl | 2 +- src/connectivity/computeRM.jl | 24 +++++++++++++----------- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/src/algorithms/paramcurve.jl b/src/algorithms/paramcurve.jl index 7868ff5..f74a061 100644 --- a/src/algorithms/paramcurve.jl +++ b/src/algorithms/paramcurve.jl @@ -1,4 +1,4 @@ -export rational_curve_parametrization +export rational_curve_parametrization, add_genvars @doc Markdown.doc""" rational_curve_parametrization(I::Ideal{T} where T <: MPolyRingElem, ) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 9257763..1979775 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -9,16 +9,6 @@ export roadmap, computepolarproj, computepolarphi, computepolar export compute_minors, detmpoly, change_ringvar_mod, MPolyBuild, small_mid_point include("Cannytools.jl") -function roadmap( - V::Ideal{P}; # input ideal - Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients - C::Vector{Vector{P}}=Vector{Vector{P}}(), # query points with rational coefficients - v::Int=0, # verbosity level - checks::Bool=false # perform checks (dimension, regularity, etc.) -) where (P <: QQMPolyRingElem) - CQ = real_solutions(C, info_level=max(v-1,0), nr_thrds=Threads.nthreads()) -end - function roadmap( V::Ideal{T} where T <: QQMPolyRingElem; # input ideal Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients @@ -95,7 +85,7 @@ function roadmap( #K1WRat = K1WRat[2:end-1] ########## Cq = [c for c in C if c[1:e]==q] - append!(K1WRat, getindex.(Cq, e+1)) |> sort! + append!(K1WRat, unique(getindex.(Cq, e+1))) |> sort! newQ = [ vcat(q, [kv]) for kv in K1WRat ] if !isempty(newQ) @@ -107,6 +97,18 @@ function roadmap( return RM end +function roadmap( + V::Ideal{P}, # input ideal + C::Ideal{P}; # query points as an ideal + v::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + println("plop") + @assert(parent(V)==parent(C), "Equations for variety and query points must be in the same ring") + CQ = real_solutions(C, info_level=max(v-1,0), nr_thrds=Threads.nthreads()) + return roadmap(V, C=CQ, v=v, checks=checks) +end + #= ## Test ## R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) From c65a16932851ba2deea860872f78fbb56efb3a5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 23 Apr 2025 19:51:46 +0200 Subject: [PATCH 06/33] MidRationalPoints: improve + add query points --- src/connectivity/Cannytools.jl | 57 ++++++++++++++++++++++++++-------- src/connectivity/computeRM.jl | 4 +-- 2 files changed, 45 insertions(+), 16 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 49dd3c6..6a7f90f 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -77,24 +77,55 @@ function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) return detA end -function MidRationalPoints(S::Vector{Vector{QQFieldElem}}) - # S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] +function MidRationalPoints(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T <: QQFieldElem} + # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] # such that the [l_i, r_i] are rational and disjoint open intervals. + # * Q is a list of rationals that intersects no [l_i,r_i] # - # It orders the [l_i,r_i], removes repetitions and computes the simplest - # rational points between these intervals. - if isempty(S) - return QQFieldElem[] + # It orders the [l_i,r_i], and compute a list ratioP such that + # strictly between each of these intervals there is: + # - either at least one element of Q + # - or the simplest rational number + isempty(S) && return Q + + S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q) + ratioP = T[] + qidx = 1 + qlen = length(Q1) + + # Handle left gap before first interval + while qidx <= qlen && Q1[qidx] < S1[1][1] + push!(ratioP, Q1[qidx]) + qidx += 1 end - S1 = sort(S, lt=(x, y) -> x[2] <= y[1]) - n = unique!(S1) |> length - ratioP = Vector{QQFieldElem}(undef, n-1) - for i in 1:(n- 1) + + # Loop through gaps between sorted disjoint intervals + for i in 1:(length(S1) - 1) ri, li1 = S1[i][2], S1[i+1][1] - @assert ri < li1 "The intervals are not disjoint." - eps = (li1-ri)//1000 - ratioP[i] = simplest_between(ri + eps, li1 - eps) + @assert ri < li1 "Intervals are not disjoint." + inserted = false + while qidx <= qlen && Q1[qidx] < li1 + push!(ratioP, Q1[qidx]) + inserted = true + qidx += 1 + end + if !inserted + eps = (li1 - ri)//1000 # for open interval + # We choose the simplest in absolute value + if -ri > li1 # this means ri is negative and the largest in absolute value + push!(ratioP, -simplest_between(-ri - eps, -li1 + eps)) + else + push!(ratioP, simplest_between(ri + eps, li1 - eps)) + end + end + end + + # Append remaining right-side Q points + while qidx <= qlen + push!(ratioP, Q1[qidx]) + qidx += 1 end + return ratioP end diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 1979775..de7a387 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,9 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolarproj, computepolarphi, computepolar -# DEBUG -export compute_minors, detmpoly, change_ringvar_mod, MPolyBuild, small_mid_point +export roadmap, computepolar include("Cannytools.jl") function roadmap( From 9b2691aaa21d63ff6dca9dd866da3ae91ae7352d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 23 Apr 2025 21:18:22 +0200 Subject: [PATCH 07/33] avoid intermediate fibers when already intermediate query points + big bug --- src/algorithms/paramcurve.jl | 2 +- src/algorithms/solvers.jl | 2 +- src/connectivity/Cannytools.jl | 1 + src/connectivity/computeRM.jl | 41 ++++++++++++++++------------------ 4 files changed, 22 insertions(+), 24 deletions(-) diff --git a/src/algorithms/paramcurve.jl b/src/algorithms/paramcurve.jl index f74a061..a4e7640 100644 --- a/src/algorithms/paramcurve.jl +++ b/src/algorithms/paramcurve.jl @@ -1,4 +1,4 @@ -export rational_curve_parametrization, add_genvars +export rational_curve_parametrization, add_genvars, change_ringvar @doc Markdown.doc""" rational_curve_parametrization(I::Ideal{T} where T <: MPolyRingElem, ) diff --git a/src/algorithms/solvers.jl b/src/algorithms/solvers.jl index ebf4539..af5a162 100644 --- a/src/algorithms/solvers.jl +++ b/src/algorithms/solvers.jl @@ -272,7 +272,7 @@ end rational_solutions(I::Ideal{T} where T <: MPolyRingElem, ) Given an ideal `I` with a finite solution set over the complex numbers, return -the rational roots of the ideal. +the rational roots of the ideal. # Arguments - `I::Ideal{T} where T <: MPolyRingElem`: input generators. diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 6a7f90f..2cf35df 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -105,6 +105,7 @@ function MidRationalPoints(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T <: @assert ri < li1 "Intervals are not disjoint." inserted = false while qidx <= qlen && Q1[qidx] < li1 + @assert(Q1[qidx] > ri, "A query point is singular") push!(ratioP, Q1[qidx]) inserted = true qidx += 1 diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index de7a387..c508184 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar +export roadmap, computepolar, MidRationalPoints include("Cannytools.jl") function roadmap( @@ -22,9 +22,6 @@ function roadmap( end # Some preprocessing V.dim == -1 && dimension(V) - println(Q) - println(C) - println() isempty(Q) && push!(Q,[]) # Base points e = length(first(Q)) @@ -50,18 +47,19 @@ function roadmap( ## sing(Fq) ## if checks - v>0 && println("Compute first the singular points") + v>0 && println("Check real quasi-smoothness") singFq = computepolar(0, Fq, v=max(v-1,0)) - @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), "Non-emtpy real sing locus!") + @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), + "Non-empty real sing locus!") end ## K(pi_1,Fq) ## - v>0 && println("First critical points") + v>0 && println("V-critical points") K1Fq = computepolar(1, Fq, v=max(v-1,0)) K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## - v>0 && println("Second critical points") + v>0 && println("Polar variety") K2Fq = computepolar(2, Fq, v=max(v-1,0)) if checks @assert(isone(dimension(K2Fq)), "Non-generic polar variety") @@ -72,20 +70,20 @@ function roadmap( push!(RM, Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) ## Points with vertical tg in K(pi_2, Fq) ## - v>0 && println("Vertical tg points") + v>0 && println("W-critical points with vertical tangent") K1WmFq = computepolar(2, K2Fq, dimproj=0, v=max(v-1,0)) K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) - ## New base points ## + ## New base and query points ## + Cq = isempty(q) ? C : [ c[2:end] for c in C if c[1] == q[e]] K1W = vcat(K1Fq, K1WmFq) - K1WRat = MidRationalPoints(getindex(K1W, 1)) - # Heuristic to be proven - #K1WRat = K1WRat[2:end-1] + # Heuristic to be proven (Reeb's th) + #K1W = K1W[2:end-1] ########## - Cq = [c for c in C if c[1:e]==q] - append!(K1WRat, unique(getindex.(Cq, e+1))) |> sort! - newQ = [ vcat(q, [kv]) for kv in K1WRat ] + K1WRat = MidRationalPoints(first.(K1W), unique(first.(Cq))) + newQ = vcat.(Ref(q), K1WRat) + # Recursively compute roadmap of possible fibers if !isempty(newQ) RMFq = roadmap(V, Q=newQ, C=Cq) append!(RM, RMFq) @@ -96,13 +94,12 @@ function roadmap( end function roadmap( - V::Ideal{P}, # input ideal - C::Ideal{P}; # query points as an ideal - v::Int=0, # verbosity level - checks::Bool=false # perform checks (dimension, regularity, etc.) + V::Ideal{P}, # input ideal + C::Ideal{P}; # ideal defining query points + v::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) ) where (P <: QQMPolyRingElem) - println("plop") - @assert(parent(V)==parent(C), "Equations for variety and query points must be in the same ring") + @assert(parent(V)==parent(C), "Equations for variety and query points must live the same ring") CQ = real_solutions(C, info_level=max(v-1,0), nr_thrds=Threads.nthreads()) return roadmap(V, C=CQ, v=v, checks=checks) end From 8c96896a4aadd205c8d59e3989d83d7c1a214c46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 23 Apr 2025 22:23:37 +0200 Subject: [PATCH 08/33] no evaluation during recursive steps of roadmap --- src/connectivity/Cannytools.jl | 5 +++-- src/connectivity/computeRM.jl | 27 +++++++++++++++------------ 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 2cf35df..57f1af9 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -3,13 +3,14 @@ function computepolar( j::Int, # j-th first coordinate images of phi V::Ideal{P}; # input ideal phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) + dim = -1, # specify dimension of V dimproj = j-1, # maximum dimension of tangent space of phi v=0, # verbosity level ) where (P <: MPolyRingElem) - V.dim == -1 && dimension(V) + dim == -1 && V.dim == -1 && dimension(V) R = parent(V) n = nvars(R) - c = n - V.dim + c = n - max(dim, V.dim) nphi = length(phi) JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in vcat(V.gens, phi[1:min(j,nphi)])]) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index c508184..e15943b 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -25,13 +25,12 @@ function roadmap( isempty(Q) && push!(Q,[]) # Base points e = length(first(Q)) - fixvarias = varias[1:e] RM = Vector{Ideal{QQMPolyRingElem}}(undef,0) for q in Q ## Fq ## - Fq = e > 0 ? Ideal(change_ringvar([evaluate(h, fixvarias, q) for h in V.gens], A.S[e+1:end])) : V # Genericity assumption (can be checked) + Fq = fbr(V, q) if checks @assert(dimension(Fq) == V.dim - e, "Non-generic polar variety") else @@ -40,47 +39,46 @@ function roadmap( # Terminal case (dim <=1) if Fq.dim <= 1 - curve = change_ringvar(Fq.gens, A.S) - push!(RM, Ideal(vcat(curve, [fixvarias[j] - q[j] for j in 1:e]))) + push!(RM, Fq) continue end ## sing(Fq) ## if checks v>0 && println("Check real quasi-smoothness") - singFq = computepolar(0, Fq, v=max(v-1,0)) + singFq = fbr(computepolar(e, V, v=max(v-1,0)), q) @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), "Non-empty real sing locus!") end ## K(pi_1,Fq) ## v>0 && println("V-critical points") - K1Fq = computepolar(1, Fq, v=max(v-1,0)) + K1Fq = fbr(computepolar(e+1, V, v=max(v-1,0)), q) K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## v>0 && println("Polar variety") - K2Fq = computepolar(2, Fq, v=max(v-1,0)) + K2Fqmins = computepolar(e+2, V, v=max(v-1,0)) + K2Fq = fbr(K2Fqmins, q) if checks @assert(isone(dimension(K2Fq)), "Non-generic polar variety") else K2Fq.dim = 1 end - polar = change_ringvar(K2Fq.gens, A.S) - push!(RM, Ideal(vcat(polar, [fixvarias[j] - q[j] for j in 1:e]))) + push!(RM, K2Fq) ## Points with vertical tg in K(pi_2, Fq) ## v>0 && println("W-critical points with vertical tangent") - K1WmFq = computepolar(2, K2Fq, dimproj=0, v=max(v-1,0)) + K1WmFq = fbr(computepolar(e+2, K2Fqmins, dim=e+1, dimproj=e+0, v=max(v-1,0)), q) K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) ## New base and query points ## - Cq = isempty(q) ? C : [ c[2:end] for c in C if c[1] == q[e]] + Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] K1W = vcat(K1Fq, K1WmFq) # Heuristic to be proven (Reeb's th) #K1W = K1W[2:end-1] ########## - K1WRat = MidRationalPoints(first.(K1W), unique(first.(Cq))) + K1WRat = MidRationalPoints(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) newQ = vcat.(Ref(q), K1WRat) # Recursively compute roadmap of possible fibers @@ -104,6 +102,11 @@ function roadmap( return roadmap(V, C=CQ, v=v, checks=checks) end +function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) + vars = gens(parent(I)) + return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) +end + #= ## Test ## R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) From 5de7bd0758b4bba33e027cd4d90e74c83183b344 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Thu, 24 Apr 2025 19:58:33 +0200 Subject: [PATCH 09/33] correct bug, improve output --- src/connectivity/Cannytools.jl | 12 ++- src/connectivity/computeRM.jl | 140 +++++++++++++++++---------------- 2 files changed, 80 insertions(+), 72 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 57f1af9..456a324 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -3,21 +3,25 @@ function computepolar( j::Int, # j-th first coordinate images of phi V::Ideal{P}; # input ideal phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) - dim = -1, # specify dimension of V dimproj = j-1, # maximum dimension of tangent space of phi v=0, # verbosity level + only_mins = false # return only minors without eqns of V ) where (P <: MPolyRingElem) - dim == -1 && V.dim == -1 && dimension(V) + V.dim == -1 && dimension(V) R = parent(V) n = nvars(R) - c = n - max(dim, V.dim) + c = n - V.dim nphi = length(phi) JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in vcat(V.gens, phi[1:min(j,nphi)])]) sizeminors = c + min(nphi,j) + min(dimproj,j-1) - (j-1) minors = compute_minors(sizeminors, JW, R) - return Ideal(vcat(V.gens, minors)) + if only_mins + return minors + else + return vcat(V.gens, minors) + end end function compute_minors(p, A, R) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index e15943b..f26f65d 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,86 +4,82 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints +export roadmap, computepolar, MidRationalPoints, fbr include("Cannytools.jl") function roadmap( - V::Ideal{T} where T <: QQMPolyRingElem; # input ideal - Q::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # base points with rational coefficients - C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients - v::Int=0, # verbosity level - checks::Bool=false # perform checks (dimension, regularity, etc.) - ) + V::Ideal{T} where T <: QQMPolyRingElem; # input ideal + q::Vector{QQFieldElem}=QQFieldElem[], # single base point with rational coefficients + C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients + v::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) # Some base cases A = parent(V) varias = gens(A) if length(varias)<=2 - return [V] + return [V.gens] end # Some preprocessing V.dim == -1 && dimension(V) - isempty(Q) && push!(Q,[]) # Base points - e = length(first(Q)) - - RM = Vector{Ideal{QQMPolyRingElem}}(undef,0) - for q in Q - ## Fq ## - # Genericity assumption (can be checked) - Fq = fbr(V, q) - if checks - @assert(dimension(Fq) == V.dim - e, "Non-generic polar variety") - else - Fq.dim = V.dim - e - end + e = length(q) - # Terminal case (dim <=1) - if Fq.dim <= 1 - push!(RM, Fq) - continue - end + RM = Vector{Tuple{Vector{QQMPolyRingElem}, Vector{QQFieldElem}}}(undef, 0) + ## Fq ## + # Genericity assumption (can be checked) + if checks + @assert(dimension(fbr(V,q)) == V.dim - e, "Non-generic coordinates") + end - ## sing(Fq) ## - if checks - v>0 && println("Check real quasi-smoothness") - singFq = fbr(computepolar(e, V, v=max(v-1,0)), q) - @assert(isempty(real_solutions(singFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads())), - "Non-empty real sing locus!") - end + # Terminal case (dim <=1) + if V.dim - e <= 1 + push!(RM, ([], q)) + return RM + end - ## K(pi_1,Fq) ## - v>0 && println("V-critical points") - K1Fq = fbr(computepolar(e+1, V, v=max(v-1,0)), q) - K1Fq = real_solutions(K1Fq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) - - ## K(pi_2, Fq) ## - v>0 && println("Polar variety") - K2Fqmins = computepolar(e+2, V, v=max(v-1,0)) - K2Fq = fbr(K2Fqmins, q) - if checks - @assert(isone(dimension(K2Fq)), "Non-generic polar variety") - else - K2Fq.dim = 1 - end - push!(RM, K2Fq) - - ## Points with vertical tg in K(pi_2, Fq) ## - v>0 && println("W-critical points with vertical tangent") - K1WmFq = fbr(computepolar(e+2, K2Fqmins, dim=e+1, dimproj=e+0, v=max(v-1,0)), q) - K1WmFq = real_solutions(K1WmFq, info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) - - ## New base and query points ## - Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] - K1W = vcat(K1Fq, K1WmFq) - # Heuristic to be proven (Reeb's th) - #K1W = K1W[2:end-1] - ########## - K1WRat = MidRationalPoints(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) - newQ = vcat.(Ref(q), K1WRat) - - # Recursively compute roadmap of possible fibers - if !isempty(newQ) - RMFq = roadmap(V, Q=newQ, C=Cq) + ## sing(Fq) ## + if checks + v>0 && println("Check real quasi-smoothness") + singFq = computepolar(e, V, v=max(v-1,0)) |> Ideal + @assert(isempty(real_solutions(fbr(singFq, q), info_level=max(v-1,0), nr_thrds=Threads.nthreads())), + "Non-empty real sing locus!") + end + + ## K(pi_1,Fq) ## + v>0 && println("V-critical points") + K1Fq = computepolar(e+1, V, v=max(v-1,0)) |> Ideal + K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## K(pi_2, Fq) ## + v>0 && println("Polar variety") + K2Fqmins = computepolar(e+2, V, v=max(v-1,0), only_mins=true) + K2Fq = vcat(V.gens, K2Fqmins) |> Ideal + if checks + @assert(fbr(K2Fq, q) |> dimension == 1, "Non-generic polar variety") + else + K2Fq.dim = e + 1 + end + push!(RM, (K2Fqmins, q)) + + ## Points with vertical tg in K(pi_2, Fq) ## + v>0 && println("W-critical points with vertical tangent") + K1WmFq = computepolar(e+2, K2Fq, dimproj=e, v=max(v-1,0)) |> Ideal + K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + + ## New base and query points ## + Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] + K1W = vcat(K1Fq, K1WmFq) + # Heuristic to be proven (Reeb's th) + #K1W = K1W[2:end-1] + ########## + K1WRat = MidRationalPoints(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) + newQ = vcat.(Ref(q), K1WRat) + + # Recursively compute roadmap of possible fibers + if !isempty(newQ) + for newq in newQ + RMFq = roadmap(V, q=newq, C=Cq) append!(RM, RMFq) end end @@ -91,6 +87,8 @@ function roadmap( return RM end + + function roadmap( V::Ideal{P}, # input ideal C::Ideal{P}; # ideal defining query points @@ -102,9 +100,15 @@ function roadmap( return roadmap(V, C=CQ, v=v, checks=checks) end + +function fbr(F::Vector{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) + @assert(!isempty(F), "Empty polynomial vector") + vars = gens(parent(first(F))) + return Ideal(vcat(F, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) +end + function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) - vars = gens(parent(I)) - return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) + return fbr(I.gens, Q) end #= From 4366ed4326d2e67474a911eac96c7ba8507d7e7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Thu, 24 Apr 2025 20:52:15 +0200 Subject: [PATCH 10/33] start Roadmap structure --- Project.toml | 1 + src/connectivity/computeRM.jl | 66 ++++++++++++++++++++++++++++++----- src/types.jl | 4 ++- 3 files changed, 61 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index e01efb9..e8efddd 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["ederc ", "Mohab Safey El Din 0 && println("W-critical points with vertical tangent") @@ -80,15 +78,17 @@ function roadmap( if !isempty(newQ) for newq in newQ RMFq = roadmap(V, q=newq, C=Cq) - append!(RM, RMFq) + push!(RM.children, RMFq) end end - return RM + if e == 0 + return Roadmap(V, RM) + else + return RM + end end - - function roadmap( V::Ideal{P}, # input ideal C::Ideal{P}; # ideal defining query points @@ -111,6 +111,54 @@ function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) return fbr(I.gens, Q) end +mutable struct RMnode + polar_eqs::Vector{QQMPolyRingElem} + base_pt::Vector{QQFieldElem} + children::Vector{RMnode} +end + +mutable struct Roadmap + initial_ideal::Ideal{QQMPolyRingElem} + root::RMnode +end + +function all_eqs(RM::Roadmap) + function all_eqs_rec(RMn::RMnode) + eqs = [fbr(vcat(RM.initial_ideal.gens, RMn.polar_eqs), RMn.base_pt)] + for child in RMn.children + append!(eqs, all_eqs_rec(child)) + end + return eqs + end + return all_eqs_rec(RM.root) +end + +function all_base_pts(RM::Roadmap) + function all_base_pts_rec(RMn::RMnode) + eqs = [RMn.base_pt] + for child in RMn.children + append!(eqs, all_base_pts_rec(child)) + end + return eqs + end + return all_base_pts_rec(RM.root) +end + +function nb_nodes(RM::Roadmap) + function nb_nodes_rec(RMn::RMnode) + nb = 1 + for child in RMn.children + nb += nb_nodes_rec(child) + end + return nb + end + return nb_nodes_rec(RM.root) +end + +Base.show(io::IO, RM::Roadmap) = print(io, all_base_pts(RM)) +Base.getindex(RM::Roadmap, idx::Union{Int, UnitRange}) = all_eqs(RM)[idx] +Base.lastindex(RM::Roadmap) = nb_nodes(RM) + #= ## Test ## R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) diff --git a/src/types.jl b/src/types.jl index a041ef9..381dc63 100644 --- a/src/types.jl +++ b/src/types.jl @@ -75,4 +75,6 @@ Base.parent(I::Ideal) = Nemo.parent(I.gens[1]) Base.show(io::IO, I::Ideal) = print(io, I.gens) -Base.getindex(I::Ideal, idx::Int) = I.gens[idx] +Base.getindex(I::Ideal, idx::Union{Int, UnitRange}) = I.gens[idx] + +Base.lastindex(I::Ideal) = lastindex(I.gens) From 4a9f502e31cd21f8c2f300065e476f1baa9bef70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Mon, 28 Apr 2025 14:46:39 +0200 Subject: [PATCH 11/33] remove abstract trees --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index e8efddd..e01efb9 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ authors = ["ederc ", "Mohab Safey El Din Date: Sat, 28 Jun 2025 23:09:59 +0200 Subject: [PATCH 12/33] maj dim nothing --- src/connectivity/computeRM.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 67192fa..8c63ea4 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -21,7 +21,7 @@ function roadmap( return [V.gens] end # Some preprocessing - V.dim == -1 && dimension(V) + isnothing(V.dim) && dimension(V) # Base points e = length(q) From ee70bc5009b7670287ec91c01e9ff384176f8f0d Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Mon, 13 Oct 2025 10:22:54 +0200 Subject: [PATCH 13/33] move RM struct to types.jl + improve internal functions --- src/connectivity/computeRM.jl | 50 +---------------------------------- src/types.jl | 37 ++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 49 deletions(-) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 8c63ea4..358012f 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts +export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes include("Cannytools.jl") function roadmap( @@ -111,54 +111,6 @@ function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) return fbr(I.gens, Q) end -mutable struct RMnode - polar_eqs::Vector{QQMPolyRingElem} - base_pt::Vector{QQFieldElem} - children::Vector{RMnode} -end - -mutable struct Roadmap - initial_ideal::Ideal{QQMPolyRingElem} - root::RMnode -end - -function all_eqs(RM::Roadmap) - function all_eqs_rec(RMn::RMnode) - eqs = [fbr(vcat(RM.initial_ideal.gens, RMn.polar_eqs), RMn.base_pt)] - for child in RMn.children - append!(eqs, all_eqs_rec(child)) - end - return eqs - end - return all_eqs_rec(RM.root) -end - -function all_base_pts(RM::Roadmap) - function all_base_pts_rec(RMn::RMnode) - eqs = [RMn.base_pt] - for child in RMn.children - append!(eqs, all_base_pts_rec(child)) - end - return eqs - end - return all_base_pts_rec(RM.root) -end - -function nb_nodes(RM::Roadmap) - function nb_nodes_rec(RMn::RMnode) - nb = 1 - for child in RMn.children - nb += nb_nodes_rec(child) - end - return nb - end - return nb_nodes_rec(RM.root) -end - -Base.show(io::IO, RM::Roadmap) = print(io, all_base_pts(RM)) -Base.getindex(RM::Roadmap, idx::Union{Int, UnitRange}) = all_eqs(RM)[idx] -Base.lastindex(RM::Roadmap) = nb_nodes(RM) - #= ## Test ## R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) diff --git a/src/types.jl b/src/types.jl index e51750e..6ef5453 100644 --- a/src/types.jl +++ b/src/types.jl @@ -77,3 +77,40 @@ Base.show(io::IO, I::Ideal) = print(io, I.gens) Base.getindex(I::Ideal, idx::Union{Int, UnitRange}) = I.gens[idx] Base.lastindex(I::Ideal) = lastindex(I.gens) + +mutable struct RMnode + polar_eqs::Vector{QQMPolyRingElem} + base_pt::Vector{QQFieldElem} + children::Vector{RMnode} +end + +mutable struct Roadmap + initial_ideal::Ideal{QQMPolyRingElem} + root::RMnode +end + +function _collect_roadmap(RMn::RMnode, F) + data = [F(RMn)] + for child in RMn.children + append!(data, _collect_roadmap(child, F)) + end + return data +end + +function all_eqs(RM::Roadmap) + func(s) = fbr(vcat(RM.initial_ideal.gens, s.polar_eqs), s.base_pt) + return _collect_roadmap(RM.root, func) +end + +function all_base_pts(RM::Roadmap) + return _collect_roadmap(RM.root, s->s.base_pt) +end + +function nb_nodes(RM::Roadmap) + return length(_collect_roadmap(RM.root, s -> true)) +end + +Base.show(io::IO, RM::Roadmap) = print(io, all_base_pts(RM)) +Base.getindex(RM::Roadmap, idx::Union{Int, UnitRange}) = all_eqs(RM)[idx] +Base.lastindex(RM::Roadmap) = nb_nodes(RM) +Base.length(RM::Roadmap) = nb_nodes(RM) \ No newline at end of file From 4a9393611e4584ed9f20ba2510dcb42240d189dc Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Wed, 15 Oct 2025 00:13:01 +0200 Subject: [PATCH 14/33] start doc + some renamings --- src/connectivity/Cannytools.jl | 7 ++- src/connectivity/computeRM.jl | 112 ++++++++++++++++++++++++--------- 2 files changed, 85 insertions(+), 34 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 456a324..2327ff1 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -4,16 +4,16 @@ function computepolar( V::Ideal{P}; # input ideal phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) dimproj = j-1, # maximum dimension of tangent space of phi - v=0, # verbosity level only_mins = false # return only minors without eqns of V ) where (P <: MPolyRingElem) - V.dim == -1 && dimension(V) + isnothing(V.dim) && dimension(V) R = parent(V) n = nvars(R) c = n - V.dim nphi = length(phi) - JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in vcat(V.gens, phi[1:min(j,nphi)])]) + psi = vcat(phi, V.gens[nphi+1:end]) + JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in psi]) sizeminors = c + min(nphi,j) + min(dimproj,j-1) - (j-1) minors = compute_minors(sizeminors, JW, R) @@ -24,6 +24,7 @@ function computepolar( end end + function compute_minors(p, A, R) #Computes the p-minors of a matrix A n, m = size(A) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 358012f..3719727 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,55 +4,105 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes +export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis include("Cannytools.jl") +@doc Markdown.doc""" + roadmap(I::Ideal{T} where T <: MPolyRingElem, ) + +Given a **radical** ideal `I` with solution set X, that is smooth and +whose real trace XR is bounded, return a roadmap of XR + +The output is given as a Roadmap structure, encoding the recursive structure +of roadmaps. It is encoded as a chained list, whose root is containing the equations defining X +and each node representing a curve component, that is defined by additional polar equation and base point< +Moreover it is linked to fibers, that share the same base point. + +# Arguments +- `I::Ideal{T} where T <: QQMPolyRingElem`: input generators. +- `q::Vector{QQFieldElem}=QQFieldElem[]`: single base point with rational coefficients +- `C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[]`: query points with rational coefficients +- `info_level::Int=0`: verbosity level +- `checks::Bool=false`: whether perform checks (dimension, regularity, etc.) +- 'generic_change=false": whether it performs a prior random linear change of variables (TODO) +) + +# Examples +```jldoctest +julia> using AlgebraicSolving + +julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) +(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4]) + +julia> I = Ideal([x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2)]) +QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64] + +julia> RM = roadmap(I) +Vector{QQFieldElem}[[], [-3], [-3, -2], [-1], [-1, -2], [-1, -1], [-1, 2], [3], [3, 2]] + +julia> nb_nodes(RM) +9 + +julia> all_eqs(RM) +9-element Vector{Ideal{QQMPolyRingElem}}: + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x3 + 4*x2^2*x3 + 4*x3^3 + 4*x3*x4^2 - 40*x3, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 + 3] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 3, x2 + 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 + 1] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 + 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 + 1] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 - 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 - 3] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 - 3, x2 - 2] + +julia> roadmap(V, generic_change=true) +..... +``` +""" function roadmap( - V::Ideal{T} where T <: QQMPolyRingElem; # input ideal + I::Ideal{T} where T <: QQMPolyRingElem; # input ideal q::Vector{QQFieldElem}=QQFieldElem[], # single base point with rational coefficients C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients - v::Int=0, # verbosity level + info_level::Int=0, # verbosity level checks::Bool=false # perform checks (dimension, regularity, etc.) ) # Some base cases - A = parent(V) - varias = gens(A) - if length(varias)<=2 - return [V.gens] + if nvars(parent(I))<=2 + return [I.gens] end # Some preprocessing - isnothing(V.dim) && dimension(V) + isnothing(I.dim) && dimension(I) # Base points e = length(q) ## Fq ## # Genericity assumption (can be checked) if checks - @assert(dimension(fbr(V,q)) == V.dim - e, "Non-generic coordinates") + @assert(dimension(fbr(I,q)) == I.dim - e, "Non-generic coordinates") end # Terminal case (dim <=1) - if V.dim - e <= 1 + if I.dim - e <= 1 return RMnode([], q, RMnode[]) end ## sing(Fq) ## if checks - v>0 && println("Check real quasi-smoothness") - singFq = computepolar(e, V, v=max(v-1,0)) |> Ideal - @assert(isempty(real_solutions(fbr(singFq, q), info_level=max(v-1,0), nr_thrds=Threads.nthreads())), + info_level>0 && println("Check real quasi-smoothness") + singFq = computepolar(e, I) |> Ideal + @assert(isempty(real_solutions(fbr(singFq, q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads())), "Non-empty real sing locus!") end ## K(pi_1,Fq) ## - v>0 && println("V-critical points") - K1Fq = computepolar(e+1, V, v=max(v-1,0)) |> Ideal - K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + info_level>0 && println("I-critical points") + K1Fq = computepolar(e+1, I) |> Ideal + K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## - v>0 && println("Polar variety") - K2Fqmins = computepolar(e+2, V, v=max(v-1,0), only_mins=true) - K2Fq = vcat(V.gens, K2Fqmins) |> Ideal + info_level>0 && println("Polar variety") + K2Fqmins = computepolar(e+2, I, only_mins=true) + K2Fq = vcat(I.gens, K2Fqmins) |> Ideal if checks @assert(fbr(K2Fq, q) |> dimension == 1, "Non-generic polar variety") else @@ -61,9 +111,9 @@ function roadmap( RM = RMnode(K2Fqmins, q, RMnode[]) ## Points with vertical tg in K(pi_2, Fq) ## - v>0 && println("W-critical points with vertical tangent") - K1WmFq = computepolar(e+2, K2Fq, dimproj=e, v=max(v-1,0)) |> Ideal - K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(v-1,0), nr_thrds=Threads.nthreads(), interval=true) + info_level>0 && println("W-critical points with vertical tangent") + K1WmFq = computepolar(e+2, K2Fq, dimproj=e) |> Ideal + K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## New base and query points ## Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] @@ -77,27 +127,27 @@ function roadmap( # Recursively compute roadmap of possible fibers if !isempty(newQ) for newq in newQ - RMFq = roadmap(V, q=newq, C=Cq) + RMFq = roadmap(I, q=newq, C=Cq) push!(RM.children, RMFq) end end if e == 0 - return Roadmap(V, RM) + return Roadmap(I, RM) else return RM end end function roadmap( - V::Ideal{P}, # input ideal + I::Ideal{P}, # input ideal C::Ideal{P}; # ideal defining query points - v::Int=0, # verbosity level + info_level::Int=0, # verbosity level checks::Bool=false # perform checks (dimension, regularity, etc.) ) where (P <: QQMPolyRingElem) - @assert(parent(V)==parent(C), "Equations for variety and query points must live the same ring") - CQ = real_solutions(C, info_level=max(v-1,0), nr_thrds=Threads.nthreads()) - return roadmap(V, C=CQ, v=v, checks=checks) + @assert(parent(I)==parent(C), "Equations for variety and query points must live the same ring") + CQ = real_solutions(C, info_level=max(info_level-1,0), nr_thrds=Threads.nthreads()) + return roadmap(I, C=CQ, info_level=info_level, checks=checks) end @@ -117,9 +167,9 @@ R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) #h = x1^2+x2^2+x3^2+x4^2-1 h = (x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) h = evaluate(h,[x1,x2,x3],[x1+rand(-10:10)*x2+rand(-10:10)*x3+rand(-10:10)*x4,x2+rand(-10:10)*x3+rand(-10:10)*x4,x3+rand(-10:10)*x4]) -V = Ideal([h]) +I = Ideal([h]) -carte = (computeRM(V, 3, [Vector{QQFieldElem}(undef,0)])) +carte = (computeRM(I, 3, [Vector{QQFieldElem}(undef,0)])) #println(carte) =# From da73da700be98c4fad8f0fbe48d4eb6ff2083fe6 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Wed, 15 Oct 2025 19:44:56 +0200 Subject: [PATCH 15/33] add index list for compute polar --- src/connectivity/Cannytools.jl | 35 ++++++++++++++++++++++------------ src/connectivity/computeRM.jl | 10 +++++----- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 2327ff1..8969f00 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -1,20 +1,32 @@ -# Compute the set of points x where phi_j(T_x(V)) has dimension < dimproj +# Compute the set of points x where, if psi in the map whose elts are the ones in [phi_1,..,phi_p, x_{p+1},...x_n] indexed in J, +# then psi(T_x(V)) has dimension < dimproj function computepolar( - j::Int, # j-th first coordinate images of phi - V::Ideal{P}; # input ideal - phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) - dimproj = j-1, # maximum dimension of tangent space of phi - only_mins = false # return only minors without eqns of V + J::Union{Vector{Int},UnitRange{Int}}, # coordinate images of [phi,proj] (as above) + V::Ideal{P}; # input ideal + phi::Vector{P} = P[], # polynomial map in consideration (completed by sufficiently many projections) + dimproj = length(J)-1, # maximum dimension of tangent space of phi + only_mins = false # return only minors without eqns of V ) where (P <: MPolyRingElem) - isnothing(V.dim) && dimension(V) + R = parent(V) n = nvars(R) - c = n - V.dim nphi = length(phi) + @assert all([1<=j<=n for j in J]) - psi = vcat(phi, V.gens[nphi+1:end]) - JW = transpose([ derivative(f, k) for k=max(j+1-nphi,0):n, f in psi]) - sizeminors = c + min(nphi,j) + min(dimproj,j-1) - (j-1) + isnothing(V.dim) && dimension(V) + c = n - V.dim + + ## Is it correct AND useful? + sort!(J) + ## + Jphi = [ j for j in J if j <= nphi ] + Jproj = setdiff(J, Jphi) + + # Construct the truncated Jacobian matrix + psi = vcat(V.gens, phi[Jphi]) + JW = transpose([ derivative(f, k) for k in setdiff(1:n, Jproj), f in psi]) + # Compute the minors + sizeminors = c + length(Jphi) + min(dimproj, length(J)-1) - (length(J)-1) minors = compute_minors(sizeminors, JW, R) if only_mins @@ -24,7 +36,6 @@ function computepolar( end end - function compute_minors(p, A, R) #Computes the p-minors of a matrix A n, m = size(A) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 3719727..08ca373 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis +export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL include("Cannytools.jl") @doc Markdown.doc""" @@ -89,19 +89,19 @@ function roadmap( ## sing(Fq) ## if checks info_level>0 && println("Check real quasi-smoothness") - singFq = computepolar(e, I) |> Ideal + singFq = computepolar(1:e, I) |> Ideal @assert(isempty(real_solutions(fbr(singFq, q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads())), "Non-empty real sing locus!") end ## K(pi_1,Fq) ## info_level>0 && println("I-critical points") - K1Fq = computepolar(e+1, I) |> Ideal + K1Fq = computepolar(1:e+1, I) |> Ideal K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## info_level>0 && println("Polar variety") - K2Fqmins = computepolar(e+2, I, only_mins=true) + K2Fqmins = computepolar(1:e+2, I, only_mins=true) K2Fq = vcat(I.gens, K2Fqmins) |> Ideal if checks @assert(fbr(K2Fq, q) |> dimension == 1, "Non-generic polar variety") @@ -112,7 +112,7 @@ function roadmap( ## Points with vertical tg in K(pi_2, Fq) ## info_level>0 && println("W-critical points with vertical tangent") - K1WmFq = computepolar(e+2, K2Fq, dimproj=e) |> Ideal + K1WmFq = computepolar(1:e+2, K2Fq, dimproj=e) |> Ideal K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## New base and query points ## From 7e67eeb7a8882eb858ee0f357af407406660a901 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Wed, 15 Oct 2025 23:01:04 +0200 Subject: [PATCH 16/33] doc --- src/connectivity/Cannytools.jl | 62 ++++++++++++++++++++++++++++++++-- src/connectivity/computeRM.jl | 3 -- 2 files changed, 60 insertions(+), 5 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 8969f00..8ceb3ce 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -1,5 +1,63 @@ -# Compute the set of points x where, if psi in the map whose elts are the ones in [phi_1,..,phi_p, x_{p+1},...x_n] indexed in J, -# then psi(T_x(V)) has dimension < dimproj +@doc Markdown.doc""" + computepolar( + J::Union{Vector{Int},UnitRange{Int}}, + V::Ideal{P}; + phi::Vector{P}=P[], + dimproj::Int=length(J)-1, + only_mins::Bool=false + ) where {P <: MPolyRingElem} + +Compute the **polar variety** associated with the map whose components are +`[phi_1, …, phi_p, x_{p+1}, …, x_n]` indexed by `J`, for an algebraic variety defined +by the ideal `V`. + +More precisely, this function computes the set of points `x` in `V` such that, +if `psi` denotes the above map, then the image of the tangent space `Tₓ(V)` under +the differential `psi` has dimension strictly less than `dimproj`. + +This is a key geometric construction in the computation of polar varieties, +used in the roadmap algorithm and critical point methods. + +# Arguments +- `J::Union{Vector{Int},UnitRange{Int}}`: + Indices of the selected coordinate functions +- `V::Ideal{P} where P <: MPolyRingElem`: + Input ideal defining the variety V on which critical loci is computed +- `phi::Vector{P}=P[]`: + Polynomial map possibly completed by projection coordinates to have a map of total length `n`. +- `dimproj::Int=length(J)-1`: + Expected maximal dimension of the image of the tangent space. + Typically equals the number of projection coordinates minus one. +- `only_mins::Bool=false`: + If `true`, only the computed minors of the Jacobian are returned; + otherwise, the output includes both the generators of `V` and those minors. + +# Returns +- If `only_mins=false` (default): + A `Vector{P}` containing the union of the generators of `V` and the computed minors. + This defines the ideal of the polar variety. +- If `only_mins=true`: + A `Vector{P}` containing only the minors (without the original equations of `V`). + +# Example +```jldoctest +julia> using AlgebraicSolving + +julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) +(Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4]) + +julia> I = Ideal([(x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2)]) +QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64] + +julia> computepolar(1:3, I, dimproj=1, phi=[x1^2+x2^2+x3^2+x4^2]) +5-element Vector{QQMPolyRingElem}: + x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64 + 4*x1^3 + 4*x1*x2^2 + 4*x1*x3^2 + 4*x1*x4^2 - 40*x1 + 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4 + 2*x1 + 2*x4 +``` +""" function computepolar( J::Union{Vector{Int},UnitRange{Int}}, # coordinate images of [phi,proj] (as above) V::Ideal{P}; # input ideal diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 08ca373..c33b907 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -54,9 +54,6 @@ julia> all_eqs(RM) QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 - 2] QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 - 3] QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 - 3, x2 - 2] - -julia> roadmap(V, generic_change=true) -..... ``` """ function roadmap( From bafd3ce2264f083dad43ab213d8c66afd13b6e45 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Sat, 18 Oct 2025 13:33:37 +0200 Subject: [PATCH 17/33] improve genericity checks + choice for determinant function (after some tests) --- src/connectivity/Cannytools.jl | 25 +++++++-- src/connectivity/computeRM.jl | 93 ++++++++++++++++++++++------------ 2 files changed, 83 insertions(+), 35 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index 8ceb3ce..a730125 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -82,7 +82,7 @@ function computepolar( # Construct the truncated Jacobian matrix psi = vcat(V.gens, phi[Jphi]) - JW = transpose([ derivative(f, k) for k in setdiff(1:n, Jproj), f in psi]) + JW = matrix(R, QQMPolyRingElem[ derivative(f, k) for f in psi, k in setdiff(1:n, Jproj)]) # Compute the minors sizeminors = c + length(Jphi) + min(dimproj, length(J)-1) - (length(J)-1) minors = compute_minors(sizeminors, JW, R) @@ -101,9 +101,24 @@ function compute_minors(p, A, R) colsmins = collect(combinations(1:m, p)) mins = Vector{eltype(A)}(undef, length(rowsmins) * length(colsmins)) k = 1 + # for performance tweaks, check if there are non-linear polynomials or if all are linear or constant + degmax = 0 + for a in A + degmax = min(2, max(degmax, total_degree(a))) + if degmax > 1 + break + end + end + # TODO: use algos from FLINT (in particular berkowitz) + # Naive function performs better for high degree or small matrices + if (degmax == 0 && p <= 2) || (degmax == 1 && p <= 6) || (degmax == 2) + detfct = s->detmpoly(s, R) + else # else use fraction-free LU from AbstractAlgebra.jl + detfct = det + end for rowsmin in rowsmins for colsmin in colsmins - mins[k] = detmpoly(A[rowsmin, colsmin], R) + mins[k] = detfct(A[rowsmin, colsmin]) k += 1 end end @@ -129,7 +144,7 @@ function combinations(a, n) return _combinations(a, n, 1, Vector{Int}([])) end -function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) +function detmpoly(A, R) # Get the size of the matrix n = size(A, 1) if n != size(A, 2) @@ -140,6 +155,10 @@ function detmpoly(A::Matrix{T} where T<:MPolyRingElem, R) return A[1, 1] end + if n == 2 + return A[1, 1] * A[2, 2] - A[1,2] * A[2, 1] + end + # Initialize the determinant polynomial detA = zero(R) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index c33b907..3278497 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL +export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly include("Cannytools.jl") @doc Markdown.doc""" @@ -67,15 +67,23 @@ function roadmap( if nvars(parent(I))<=2 return [I.gens] end - # Some preprocessing - isnothing(I.dim) && dimension(I) - # Base points + + # Some preprocessng + if isnothing(I.dim) + lucky_prime = _generate_lucky_primes(I.gens, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(lucky_prime)), I.gens)) + dimension(INEW) + I.dim = INEW.dim + end e = length(q) ## Fq ## # Genericity assumption (can be checked) if checks - @assert(dimension(fbr(I,q)) == I.dim - e, "Non-generic coordinates") + local Fnew = fbr(I,q).gens + new_lucky_prime = _generate_lucky_primes(Fnew, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), Fnew)) + @assert(dimension(INEW) == I.dim - e, "Non-generic coordinates") end # Terminal case (dim <=1) @@ -85,30 +93,34 @@ function roadmap( ## sing(Fq) ## if checks - info_level>0 && println("Check real quasi-smoothness") - singFq = computepolar(1:e, I) |> Ideal - @assert(isempty(real_solutions(fbr(singFq, q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads())), - "Non-empty real sing locus!") + info_level>0 && println("Check smoothness") + local FNEW = fbr(computepolar(1:e, I)|> Ideal, q).gens + new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) + @assert(dimension(INEW) == -1, "Non-empty sing locus!") end ## K(pi_1,Fq) ## - info_level>0 && println("I-critical points") + info_level>0 && println("Compute x1-critical points: K1") K1Fq = computepolar(1:e+1, I) |> Ideal K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## - info_level>0 && println("Polar variety") + info_level>0 && println("Compute (x1,x2)-polar variety: W") K2Fqmins = computepolar(1:e+2, I, only_mins=true) K2Fq = vcat(I.gens, K2Fqmins) |> Ideal if checks - @assert(fbr(K2Fq, q) |> dimension == 1, "Non-generic polar variety") + local FNEW = fbr(K2Fq, q).gens + new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first + local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) + @assert(dimension(INEW) == 1, "Non-generic polar variety") else K2Fq.dim = e + 1 end RM = RMnode(K2Fqmins, q, RMnode[]) ## Points with vertical tg in K(pi_2, Fq) ## - info_level>0 && println("W-critical points with vertical tangent") + info_level>0 && println("Compute W-critical points with vertical tangent: K1W") K1WmFq = computepolar(1:e+2, K2Fq, dimproj=e) |> Ideal K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) @@ -148,26 +160,43 @@ function roadmap( end -function fbr(F::Vector{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) - @assert(!isempty(F), "Empty polynomial vector") - vars = gens(parent(first(F))) - return Ideal(vcat(F, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) -end - function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) - return fbr(I.gens, Q) + @assert(!isempty(I.gens), "Empty polynomial vector") + vars = gens(parent(first(I.gens))) + return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) end -#= -## Test ## -R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) -#h = x1^2+x2^2+x3^2+x4^2-1 -h = (x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) -h = evaluate(h,[x1,x2,x3],[x1+rand(-10:10)*x2+rand(-10:10)*x3+rand(-10:10)*x4,x2+rand(-10:10)*x3+rand(-10:10)*x4,x3+rand(-10:10)*x4]) -I = Ideal([h]) - -carte = (computeRM(I, 3, [Vector{QQFieldElem}(undef,0)])) -#println(carte) -=# - +# Generate N primes > start that do not divide any numerator/denominator +# of any coefficient in polynomials from LP +function _generate_lucky_primes( + LF::Vector{P} where P<:MPolyRingElem, + low::ZZRingElem, + up::ZZRingElem, + N::Int64 + ) + # Avoid repetitive enumeration and redundant divisibility check + CF = ZZRingElem[] + for f in LF, c in coefficients(f), part in (numerator(c), denominator(c)) + if !isone(part) + push!(CF, part) + end + end + sort!(CF, rev=true) + unique!(CF) + + # Test primes + Lprim = ZZRingElem[] + while length(Lprim) < N + cur_prim = next_prime(rand(low:up)) + is_lucky = !(cur_prim in Lprim) + i = firstindex(CF) + # Exploit decreasing order of CF + while is_lucky && i <= lastindex(CF) && CF[i] > cur_prim + is_lucky = !is_divisible_by(CF[i], cur_prim) + i += 1 + end + is_lucky && push!(Lprim, cur_prim) + end + return Lprim +end From 8486af01bc22c44e50398a53d1ba2185357d70a5 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Sat, 18 Oct 2025 13:39:57 +0200 Subject: [PATCH 18/33] internalize functions that must be --- src/connectivity/Cannytools.jl | 34 ++++++++--------- src/connectivity/computeRM.jl | 69 ++++++++++++++++++++++++++++++---- 2 files changed, 76 insertions(+), 27 deletions(-) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/Cannytools.jl index a730125..897925f 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/Cannytools.jl @@ -85,7 +85,7 @@ function computepolar( JW = matrix(R, QQMPolyRingElem[ derivative(f, k) for f in psi, k in setdiff(1:n, Jproj)]) # Compute the minors sizeminors = c + length(Jphi) + min(dimproj, length(J)-1) - (length(J)-1) - minors = compute_minors(sizeminors, JW, R) + minors = _compute_minors(sizeminors, JW, R) if only_mins return minors @@ -94,11 +94,11 @@ function computepolar( end end -function compute_minors(p, A, R) +function _compute_minors(p, A, R) #Computes the p-minors of a matrix A n, m = size(A) - rowsmins = collect(combinations(1:n, p)) - colsmins = collect(combinations(1:m, p)) + rowsmins = collect(_combinations(1:n, p, 1, Vector{Int}([]))) + colsmins = collect(_combinations(1:m, p, 1, Vector{Int}([]))) mins = Vector{eltype(A)}(undef, length(rowsmins) * length(colsmins)) k = 1 # for performance tweaks, check if there are non-linear polynomials or if all are linear or constant @@ -126,22 +126,18 @@ function compute_minors(p, A, R) return mins end -function combinations(a, n) - # Helper function to recursively generate combinations - function _combinations(a, n, start, chosen) - if length(chosen) == n - return [chosen] - elseif start > length(a) - return Vector{Int}([]) - else - # Include the current element and recurse - include_current = _combinations(a, n, start + 1, [chosen; a[start]]) - # Exclude the current element and recurse - exclude_current = _combinations(a, n, start + 1, chosen) - return vcat(include_current, exclude_current) - end +function _combinations(a, n, start, chosen) + if length(chosen) == n + return [chosen] + elseif start > length(a) + return Vector{Int}([]) + else + # Include the current element and recurse + include_current = _combinations(a, n, start + 1, [chosen; a[start]]) + # Exclude the current element and recurse + exclude_current = _combinations(a, n, start + 1, chosen) + return vcat(include_current, exclude_current) end - return _combinations(a, n, 1, Vector{Int}([])) end function detmpoly(A, R) diff --git a/src/connectivity/computeRM.jl b/src/connectivity/computeRM.jl index 3278497..2031b36 100644 --- a/src/connectivity/computeRM.jl +++ b/src/connectivity/computeRM.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, MidRationalPoints, fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly +export roadmap, computepolar, _mid_rational_points, _fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly include("Cannytools.jl") @doc Markdown.doc""" @@ -80,7 +80,7 @@ function roadmap( ## Fq ## # Genericity assumption (can be checked) if checks - local Fnew = fbr(I,q).gens + local Fnew = _fbr(I,q).gens new_lucky_prime = _generate_lucky_primes(Fnew, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), Fnew)) @assert(dimension(INEW) == I.dim - e, "Non-generic coordinates") @@ -94,7 +94,7 @@ function roadmap( ## sing(Fq) ## if checks info_level>0 && println("Check smoothness") - local FNEW = fbr(computepolar(1:e, I)|> Ideal, q).gens + local FNEW = _fbr(computepolar(1:e, I)|> Ideal, q).gens new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) @assert(dimension(INEW) == -1, "Non-empty sing locus!") @@ -103,14 +103,14 @@ function roadmap( ## K(pi_1,Fq) ## info_level>0 && println("Compute x1-critical points: K1") K1Fq = computepolar(1:e+1, I) |> Ideal - K1Fq = real_solutions(fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) + K1Fq = real_solutions(_fbr(K1Fq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## K(pi_2, Fq) ## info_level>0 && println("Compute (x1,x2)-polar variety: W") K2Fqmins = computepolar(1:e+2, I, only_mins=true) K2Fq = vcat(I.gens, K2Fqmins) |> Ideal if checks - local FNEW = fbr(K2Fq, q).gens + local FNEW = _fbr(K2Fq, q).gens new_lucky_prime = _generate_lucky_primes(FNEW, one(ZZ)<<30, one(ZZ)<<31-1, 1) |> first local INEW = Ideal(change_base_ring.(Ref(GF(new_lucky_prime)), FNEW)) @assert(dimension(INEW) == 1, "Non-generic polar variety") @@ -122,7 +122,7 @@ function roadmap( ## Points with vertical tg in K(pi_2, Fq) ## info_level>0 && println("Compute W-critical points with vertical tangent: K1W") K1WmFq = computepolar(1:e+2, K2Fq, dimproj=e) |> Ideal - K1WmFq = real_solutions(fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) + K1WmFq = real_solutions(_fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## New base and query points ## Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] @@ -130,7 +130,7 @@ function roadmap( # Heuristic to be proven (Reeb's th) #K1W = K1W[2:end-1] ########## - K1WRat = MidRationalPoints(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) + K1WRat = _mid_rational_points(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) newQ = vcat.(Ref(q), K1WRat) # Recursively compute roadmap of possible fibers @@ -160,12 +160,65 @@ function roadmap( end -function fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) +function _fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) @assert(!isempty(I.gens), "Empty polynomial vector") vars = gens(parent(first(I.gens))) return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) end +function _mid_rational_points(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T <: QQFieldElem} + # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] + # such that the [l_i, r_i] are rational and disjoint open intervals. + # * Q is a list of rationals that intersects no [l_i,r_i] + # + # It orders the [l_i,r_i], and compute a list ratioP such that + # strictly between each of these intervals there is: + # - either at least one element of Q + # - or the simplest rational number + isempty(S) && return Q + + S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q) + ratioP = T[] + qidx = 1 + qlen = length(Q1) + + # Handle left gap before first interval + while qidx <= qlen && Q1[qidx] < S1[1][1] + push!(ratioP, Q1[qidx]) + qidx += 1 + end + + # Loop through gaps between sorted disjoint intervals + for i in 1:(length(S1) - 1) + ri, li1 = S1[i][2], S1[i+1][1] + @assert ri < li1 "Intervals are not disjoint." + inserted = false + while qidx <= qlen && Q1[qidx] < li1 + @assert(Q1[qidx] > ri, "A query point is singular") + push!(ratioP, Q1[qidx]) + inserted = true + qidx += 1 + end + if !inserted + eps = (li1 - ri)//1000 # for open interval + # We choose the simplest in absolute value + if -ri > li1 # this means ri is negative and the largest in absolute value + push!(ratioP, -simplest_between(-ri - eps, -li1 + eps)) + else + push!(ratioP, simplest_between(ri + eps, li1 - eps)) + end + end + end + + # Append remaining right-side Q points + while qidx <= qlen + push!(ratioP, Q1[qidx]) + qidx += 1 + end + + return ratioP +end + # Generate N primes > start that do not divide any numerator/denominator # of any coefficient in polynomials from LP From 731b7d4ec4e5c5e157071a8f25c0dc7273004d01 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Sat, 18 Oct 2025 13:42:05 +0200 Subject: [PATCH 19/33] rename files --- src/connectivity/{Cannytools.jl => polar.jl} | 57 +------------------ src/connectivity/{computeRM.jl => roadmap.jl} | 0 2 files changed, 3 insertions(+), 54 deletions(-) rename src/connectivity/{Cannytools.jl => polar.jl} (75%) rename src/connectivity/{computeRM.jl => roadmap.jl} (100%) diff --git a/src/connectivity/Cannytools.jl b/src/connectivity/polar.jl similarity index 75% rename from src/connectivity/Cannytools.jl rename to src/connectivity/polar.jl index 897925f..2937a22 100644 --- a/src/connectivity/Cannytools.jl +++ b/src/connectivity/polar.jl @@ -112,7 +112,7 @@ function _compute_minors(p, A, R) # TODO: use algos from FLINT (in particular berkowitz) # Naive function performs better for high degree or small matrices if (degmax == 0 && p <= 2) || (degmax == 1 && p <= 6) || (degmax == 2) - detfct = s->detmpoly(s, R) + detfct = s->_detmpoly(s, R) else # else use fraction-free LU from AbstractAlgebra.jl detfct = det end @@ -140,7 +140,7 @@ function _combinations(a, n, start, chosen) end end -function detmpoly(A, R) +function _detmpoly(A, R) # Get the size of the matrix n = size(A, 1) if n != size(A, 2) @@ -161,62 +161,11 @@ function detmpoly(A, R) # Compute the determinant polynomial for j = 1:n submatrix = A[2:end, [i for i = 1:n if i != j]] - detA += (-1)^(1+j)*A[1, j] * detmpoly(submatrix, R) + detA += (-1)^(1+j)*A[1, j] * _detmpoly(submatrix, R) end return detA end -function MidRationalPoints(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T <: QQFieldElem} - # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] - # such that the [l_i, r_i] are rational and disjoint open intervals. - # * Q is a list of rationals that intersects no [l_i,r_i] - # - # It orders the [l_i,r_i], and compute a list ratioP such that - # strictly between each of these intervals there is: - # - either at least one element of Q - # - or the simplest rational number - isempty(S) && return Q - - S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q) - ratioP = T[] - qidx = 1 - qlen = length(Q1) - - # Handle left gap before first interval - while qidx <= qlen && Q1[qidx] < S1[1][1] - push!(ratioP, Q1[qidx]) - qidx += 1 - end - - # Loop through gaps between sorted disjoint intervals - for i in 1:(length(S1) - 1) - ri, li1 = S1[i][2], S1[i+1][1] - @assert ri < li1 "Intervals are not disjoint." - inserted = false - while qidx <= qlen && Q1[qidx] < li1 - @assert(Q1[qidx] > ri, "A query point is singular") - push!(ratioP, Q1[qidx]) - inserted = true - qidx += 1 - end - if !inserted - eps = (li1 - ri)//1000 # for open interval - # We choose the simplest in absolute value - if -ri > li1 # this means ri is negative and the largest in absolute value - push!(ratioP, -simplest_between(-ri - eps, -li1 + eps)) - else - push!(ratioP, simplest_between(ri + eps, li1 - eps)) - end - end - end - - # Append remaining right-side Q points - while qidx <= qlen - push!(ratioP, Q1[qidx]) - qidx += 1 - end - return ratioP -end diff --git a/src/connectivity/computeRM.jl b/src/connectivity/roadmap.jl similarity index 100% rename from src/connectivity/computeRM.jl rename to src/connectivity/roadmap.jl From 6dd94b8beb21ce95b847cceb21554c1df613d18c Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Sat, 18 Oct 2025 14:11:57 +0200 Subject: [PATCH 20/33] change name of roadmap fct + remove code duplication --- src/AlgebraicSolving.jl | 2 +- src/connectivity/roadmap.jl | 76 ++++++++++++------------------------- 2 files changed, 26 insertions(+), 52 deletions(-) diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index 612d3b1..e3fa4ff 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -19,7 +19,7 @@ include("algorithms/hilbert.jl") #= siggb =# include("siggb/siggb.jl") #= connectivity =# -include("connectivity/computeRM.jl") +include("connectivity/roadmap.jl") #= examples =# include("examples/katsura.jl") include("examples/cyclic.jl") diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 2031b36..98d2f94 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -5,7 +5,7 @@ #using Nemo export roadmap, computepolar, _mid_rational_points, _fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly -include("Cannytools.jl") +include("polar.jl") @doc Markdown.doc""" roadmap(I::Ideal{T} where T <: MPolyRingElem, ) @@ -20,7 +20,6 @@ Moreover it is linked to fibers, that share the same base point. # Arguments - `I::Ideal{T} where T <: QQMPolyRingElem`: input generators. -- `q::Vector{QQFieldElem}=QQFieldElem[]`: single base point with rational coefficients - `C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[]`: query points with rational coefficients - `info_level::Int=0`: verbosity level - `checks::Bool=false`: whether perform checks (dimension, regularity, etc.) @@ -57,11 +56,31 @@ julia> all_eqs(RM) ``` """ function roadmap( - I::Ideal{T} where T <: QQMPolyRingElem; # input ideal - q::Vector{QQFieldElem}=QQFieldElem[], # single base point with rational coefficients + I::Ideal{P}; # input ideal C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients - info_level::Int=0, # verbosity level + info_level::Int=0, # verbosity level checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + return _roadmap_rec(I, QQFieldElem[], C, info_level, checks) +end + +function roadmap( + I::Ideal{P}, # input ideal + C::Ideal{P}; # ideal defining query points + info_level::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) +) where (P <: QQMPolyRingElem) + @assert(parent(I)==parent(C), "Equations for variety and query points must live the same ring") + CI = real_solutions(C, info_level=max(info_level-1,0), nr_thrds=Threads.nthreads()) + return _roadmap_rec(I, QFieldElem[], CI, info_level, checks) +end + +function _roadmap_rec( + I::Ideal{T} where T <: QQMPolyRingElem, # input ideal + q::Vector{QQFieldElem}, # single base point with rational coefficients + C::Vector{Vector{QQFieldElem}}, # query points with rational coefficients + info_level::Int, # verbosity level + checks::Bool # perform checks (dimension, regularity, etc.) ) # Some base cases if nvars(parent(I))<=2 @@ -136,7 +155,7 @@ function roadmap( # Recursively compute roadmap of possible fibers if !isempty(newQ) for newq in newQ - RMFq = roadmap(I, q=newq, C=Cq) + RMFq = _roadmap_rec(I, newq, Cq, info_level, checks) push!(RM.children, RMFq) end end @@ -148,18 +167,6 @@ function roadmap( end end -function roadmap( - I::Ideal{P}, # input ideal - C::Ideal{P}; # ideal defining query points - info_level::Int=0, # verbosity level - checks::Bool=false # perform checks (dimension, regularity, etc.) -) where (P <: QQMPolyRingElem) - @assert(parent(I)==parent(C), "Equations for variety and query points must live the same ring") - CQ = real_solutions(C, info_level=max(info_level-1,0), nr_thrds=Threads.nthreads()) - return roadmap(I, C=CQ, info_level=info_level, checks=checks) -end - - function _fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) @assert(!isempty(I.gens), "Empty polynomial vector") vars = gens(parent(first(I.gens))) @@ -220,36 +227,3 @@ function _mid_rational_points(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T end -# Generate N primes > start that do not divide any numerator/denominator -# of any coefficient in polynomials from LP -function _generate_lucky_primes( - LF::Vector{P} where P<:MPolyRingElem, - low::ZZRingElem, - up::ZZRingElem, - N::Int64 - ) - # Avoid repetitive enumeration and redundant divisibility check - CF = ZZRingElem[] - for f in LF, c in coefficients(f), part in (numerator(c), denominator(c)) - if !isone(part) - push!(CF, part) - end - end - sort!(CF, rev=true) - unique!(CF) - - # Test primes - Lprim = ZZRingElem[] - while length(Lprim) < N - cur_prim = next_prime(rand(low:up)) - is_lucky = !(cur_prim in Lprim) - i = firstindex(CF) - # Exploit decreasing order of CF - while is_lucky && i <= lastindex(CF) && CF[i] > cur_prim - is_lucky = !is_divisible_by(CF[i], cur_prim) - i += 1 - end - is_lucky && push!(Lprim, cur_prim) - end - return Lprim -end From d2646af9ca6cc882eba877a24ca28bda97bf8fd8 Mon Sep 17 00:00:00 2001 From: Remi Prebet Date: Sat, 18 Oct 2025 15:01:44 +0200 Subject: [PATCH 21/33] new function for interval query points --- src/connectivity/roadmap.jl | 64 ++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 98d2f94..82d72f8 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, _mid_rational_points, _fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly +export roadmap, computepolar, _mid_rational_points, _mid_rational_points_inter, _fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly include("polar.jl") @doc Markdown.doc""" @@ -226,4 +226,66 @@ function _mid_rational_points(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T return ratioP end +function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = Vector{T}[]) where {T <: QQFieldElem} + # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] + # such that the [l_i, r_i] are rational and disjoint open intervals. + # * Same assumptions on Q + # * Intervals in S and Q do not intersect as well + # + # It orders the [l_i,r_i], and compute a list ratioP such that + # strictly between each of these intervals there is: + # - either at least one element inside an interval of Q + # - or the simplest rational number + # TODO: + isempty(S) && return Q + + S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q, lt=(x, y) -> x[2] <= y[1]) + ratioP = T[] + qidx = 1 + qlen = length(Q1) + + # Handle left gap before first interval + while qidx <= qlen && Q1[qidx][2] < S1[1][1] + ql, qr = Q1[qidx] + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + qidx += 1 + end + + # Loop through gaps between sorted disjoint intervals + for i in 1:(length(S1) - 1) + ri, li1 = S1[i][2], S1[i+1][1] + @assert ri < li1 "Intervals are not disjoint." + inserted = false + while qidx <= qlen && Q1[qidx][2] < li1 + ql, qr = Q1[qidx] + @assert(ql > ri, "A query point might be singular") + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + inserted = true + qidx += 1 + end + @assert qidx > qlen || Q1[qidx][1] > S1[i+1][2] "A query point might be singular" + # If there's already rational betwee no need to add new + !inserted && push!(ratioP, _open_simplest_between(ri, li1, abs(li1 - ri)//1000)) + end + + # Append remaining right-side Q points + while qidx <= qlen + ql, qr = Q1[qidx] + push!(ratioP, _open_simplest_between(ql, qr, abs(qr - ql)//1000)) + qidx += 1 + end + + return ratioP +end + + +function _open_simplest_between(a::QQFieldElem, b::QQFieldElem, eps::QQFieldElem) + # We choose the simplest in absolute value + if -a > b # this means a is negative and the largest in absolute value + return -simplest_between(-a - eps, -b + eps) + else + return simplest_between( a + eps, b - eps) + end +end + From e8381ead2def3432c5a59e9173928e0efbcbed8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Tue, 21 Oct 2025 19:09:45 +0200 Subject: [PATCH 22/33] improve and simplify minor computations --- src/connectivity/polar.jl | 93 ++++++++----------------------------- src/connectivity/roadmap.jl | 2 +- 2 files changed, 21 insertions(+), 74 deletions(-) diff --git a/src/connectivity/polar.jl b/src/connectivity/polar.jl index 2937a22..8913967 100644 --- a/src/connectivity/polar.jl +++ b/src/connectivity/polar.jl @@ -74,8 +74,6 @@ function computepolar( isnothing(V.dim) && dimension(V) c = n - V.dim - ## Is it correct AND useful? - sort!(J) ## Jphi = [ j for j in J if j <= nphi ] Jproj = setdiff(J, Jphi) @@ -85,7 +83,7 @@ function computepolar( JW = matrix(R, QQMPolyRingElem[ derivative(f, k) for f in psi, k in setdiff(1:n, Jproj)]) # Compute the minors sizeminors = c + length(Jphi) + min(dimproj, length(J)-1) - (length(J)-1) - minors = _compute_minors(sizeminors, JW, R) + minors = _compute_minors(sizeminors, JW) if only_mins return minors @@ -94,78 +92,27 @@ function computepolar( end end -function _compute_minors(p, A, R) - #Computes the p-minors of a matrix A - n, m = size(A) - rowsmins = collect(_combinations(1:n, p, 1, Vector{Int}([]))) - colsmins = collect(_combinations(1:m, p, 1, Vector{Int}([]))) - mins = Vector{eltype(A)}(undef, length(rowsmins) * length(colsmins)) - k = 1 - # for performance tweaks, check if there are non-linear polynomials or if all are linear or constant - degmax = 0 - for a in A - degmax = min(2, max(degmax, total_degree(a))) - if degmax > 1 - break - end - end - # TODO: use algos from FLINT (in particular berkowitz) - # Naive function performs better for high degree or small matrices - if (degmax == 0 && p <= 2) || (degmax == 1 && p <= 6) || (degmax == 2) - detfct = s->_detmpoly(s, R) - else # else use fraction-free LU from AbstractAlgebra.jl - detfct = det - end - for rowsmin in rowsmins - for colsmin in colsmins - mins[k] = detfct(A[rowsmin, colsmin]) - k += 1 - end - end - - return mins +function _compute_minors(p, A) + # Computes the p-minors of a matrix A + rowsmins = _combinations(1:nrows(A), p) + colsmins = _combinations(1:ncols(A), p) + # We use charpoly for a division-free determinant method + return [ coeff(charpoly(A[rows, cols]), 0) for rows in rowsmins for cols in colsmins ] end -function _combinations(a, n, start, chosen) - if length(chosen) == n - return [chosen] - elseif start > length(a) - return Vector{Int}([]) - else - # Include the current element and recurse - include_current = _combinations(a, n, start + 1, [chosen; a[start]]) - # Exclude the current element and recurse - exclude_current = _combinations(a, n, start + 1, chosen) - return vcat(include_current, exclude_current) - end +function _combinations(v::UnitRange{Int}, k::Int) + # Compute the k-subsets of v + n = length(v) + ans = Vector{Int}[] + k > n && return ans + _combinations_dfs!(ans, Vector{Int}(undef, k), v, n, k) + return ans end -function _detmpoly(A, R) - # Get the size of the matrix - n = size(A, 1) - if n != size(A, 2) - throw(ArgumentError("Matrix must be square")) - end - - if n == 1 - return A[1, 1] - end - - if n == 2 - return A[1, 1] * A[2, 2] - A[1,2] * A[2, 1] +function _combinations_dfs!(ans::Vector{Vector{Int}}, comb::Vector{Int}, v::UnitRange{Int}, n::Int, k::Int) + k < 1 && (pushfirst!(ans, comb[:]); return) + for m in n:-1:k + comb[k] = v[m] + _combinations_dfs!(ans, comb, v, m - 1, k - 1) end - - # Initialize the determinant polynomial - detA = zero(R) - - # Compute the determinant polynomial - for j = 1:n - submatrix = A[2:end, [i for i = 1:n if i != j]] - detA += (-1)^(1+j)*A[1, j] * _detmpoly(submatrix, R) - end - - return detA -end - - - +end \ No newline at end of file diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 82d72f8..21be64e 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -4,7 +4,7 @@ #using AlgebraicSolving #using Nemo -export roadmap, computepolar, _mid_rational_points, _mid_rational_points_inter, _fbr, all_eqs, all_base_pts, nb_nodes, compute_minors, compute_minors_bis, computepolarL, detmpoly +export roadmap, computepolar, _mid_rational_points, _mid_rational_points_inter, _fbr, all_eqs, all_base_pts, nb_nodes, computepolar include("polar.jl") @doc Markdown.doc""" From ef3487822632751931475dc34e2826424c221abb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 22 Oct 2025 15:34:46 +0200 Subject: [PATCH 23/33] complete the doc --- docs/make.jl | 6 +- docs/src/connectivity.md | 48 ++++++++++++++++ docs/src/types.md | 14 ++--- src/algorithms/solvers.jl | 2 +- src/connectivity/polar.jl | 24 +++----- src/connectivity/roadmap.jl | 108 +++++++++--------------------------- src/types.jl | 8 ++- 7 files changed, 104 insertions(+), 106 deletions(-) create mode 100644 docs/src/connectivity.md diff --git a/docs/make.jl b/docs/make.jl index 9c761f2..8a25b1b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,9 @@ +using Pkg +Pkg.activate("../../AlgebraicSolving.jl") using AlgebraicSolving using Documenter + push!(LOAD_PATH, "../src") DocMeta.setdocmeta!(AlgebraicSolving, :DocTestSetup, :(using AlgebraicSolving); recursive=true) @@ -24,7 +27,8 @@ makedocs( "hilbert.md", "dimension.md", "solvers.md", - "decomposition.md"], + "decomposition.md", + "connectivity.md"], "Examples" => "katsura.md" ] ) diff --git a/docs/src/connectivity.md b/docs/src/connectivity.md new file mode 100644 index 0000000..a932210 --- /dev/null +++ b/docs/src/connectivity.md @@ -0,0 +1,48 @@ +```@meta +CurrentModule = AlgebraicSolving +DocTestSetup = quote + using AlgebraicSolving +end +``` + +```@setup algebraicsolving +using AlgebraicSolving +``` + +```@contents +Pages = ["connectivity.md"] +``` + +# Algorithms for connectivity queries on smooth bounded algebraic sets + +## Introduction + +AlgebraicSolving allows to compute a roadmap for the real trace of the zero-set +of the ideal spanned by given input generators over the rationals. + +It assumes that the underlying algebraic set is **smooth**, and its real trace is **bounded**. + +The underlying engine is provided by msolve. + +## Functionality + +```@docs + roadmap( + I::Ideal{P} where P <: QQMPolyRingElem; + C::Vector{Vector{Vector{QQFieldElem}}}=Vector{Vector{QQFieldElem}}[], + info_level::Int=0, + checks::Bool=false + ) +``` + +In addition, AlgebraicSolving can compute equations definition critical loci of polynomial maps over the given algebraic set. + +```@docs + computepolar( + J::Union{Vector{Int},UnitRange{Int}}, + V::Ideal{P}; + phi::Vector{P} = P[], + dimproj = length(J)-1, + only_mins = false + ) where (P <: MPolyRingElem) +``` \ No newline at end of file diff --git a/docs/src/types.md b/docs/src/types.md index b60edce..bc4a2a0 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -17,21 +17,21 @@ Pages = ["types.md"] ## Introduction -AlgebraicSolving handles ideals in multivariate polynomial rings over a prime +AlgebraicSolving handles ideals in multivariate polynomial rings over a prime field of characteristic $0$ or $p$ where $p$ is a prime number $<2^{31}$. ## Polynomial Rings -We use [Nemo](https://www.nemocas.org/)'s multivariate polynomial +We use [Nemo](https://www.nemocas.org/)'s multivariate polynomial ring structures: ```@repl using AlgebraicSolving R, (x,y,z) = polynomial_ring(QQ, ["x", "y", "z"], internal_ordering=:degrevlex) ``` -The above example defines a multivariate polynomial ring in three variables `x`, -`y`, and `z` over the rationals using the dgree reverse lexicographical ordering -for printing polynomials in the following. One can also define polynomial rings +The above example defines a multivariate polynomial ring in three variables `x`, +`y`, and `z` over the rationals using the dgree reverse lexicographical ordering +for printing polynomials in the following. One can also define polynomial rings over finite fields: ```@repl @@ -41,8 +41,8 @@ R, (x,y,z) = polynomial_ring(GF(101), ["x", "y", "z"], internal_ordering=:degrev ## Ideals -Ideals can be constructed by giving an array of generators. Ideals cache varies -data structures connected to ideals in order to make computational algebra more +Ideals can be constructed by giving an array of generators. Ideals cache varies +data structures connected to ideals in order to make computational algebra more effective: ```@repl diff --git a/src/algorithms/solvers.jl b/src/algorithms/solvers.jl index 52c9445..234e84e 100644 --- a/src/algorithms/solvers.jl +++ b/src/algorithms/solvers.jl @@ -295,8 +295,8 @@ QQMPolyRingElem[x1 + 2*x2 + 2*x3 - 1, x1^2 - x1 + 2*x2^2 + 2*x3^2, 2*x1*x2 + 2*x julia> rat_sols = rational_solutions(I) 2-element Vector{Vector{QQFieldElem}}: - [1, 0, 0] [1//3, 0, 1//3] + [1, 0, 0] julia> map(r->map(p->evaluate(p, r), I.gens), rat_sols) 2-element Vector{Vector{QQFieldElem}}: diff --git a/src/connectivity/polar.jl b/src/connectivity/polar.jl index 8913967..ad0805c 100644 --- a/src/connectivity/polar.jl +++ b/src/connectivity/polar.jl @@ -1,11 +1,5 @@ @doc Markdown.doc""" - computepolar( - J::Union{Vector{Int},UnitRange{Int}}, - V::Ideal{P}; - phi::Vector{P}=P[], - dimproj::Int=length(J)-1, - only_mins::Bool=false - ) where {P <: MPolyRingElem} + computepolar(J::Union{Vector{Int},UnitRange{Int}}, V::Ideal{P} where P <: MPolyRingElem; ) Compute the **polar variety** associated with the map whose components are `[phi_1, …, phi_p, x_{p+1}, …, x_n]` indexed by `J`, for an algebraic variety defined @@ -46,16 +40,16 @@ julia> using AlgebraicSolving julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) (Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4]) -julia> I = Ideal([(x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2)]) -QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64] +julia> I = Ideal([(x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) + 1]) +QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65] julia> computepolar(1:3, I, dimproj=1, phi=[x1^2+x2^2+x3^2+x4^2]) 5-element Vector{QQMPolyRingElem}: - x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64 - 4*x1^3 + 4*x1*x2^2 + 4*x1*x3^2 + 4*x1*x4^2 - 40*x1 - 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4 - 2*x1 - 2*x4 + x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65 + -4*x1^3 - 4*x1*x2^2 - 4*x1*x3^2 - 4*x1*x4^2 + 40*x1 + -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4 + -2*x1 + -2*x4 ``` """ function computepolar( @@ -100,7 +94,7 @@ function _compute_minors(p, A) return [ coeff(charpoly(A[rows, cols]), 0) for rows in rowsmins for cols in colsmins ] end -function _combinations(v::UnitRange{Int}, k::Int) +function _combinations(v::UnitRange{Int}, k::Int) # Compute the k-subsets of v n = length(v) ans = Vector{Int}[] diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 21be64e..a49e079 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -31,39 +31,43 @@ Moreover it is linked to fibers, that share the same base point. julia> using AlgebraicSolving julia> R,(x1,x2,x3,x4) = polynomial_ring(QQ, ["x1","x2","x3","x4"]) -(Multivariate polynomial ring in 3 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4]) +(Multivariate polynomial ring in 4 variables over QQ, QQMPolyRingElem[x1, x2, x3, x4]) -julia> I = Ideal([x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2)]) -QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64] +julia> I = Ideal([(x1^2+x2^2+x3^2+x4^2+9-1)^2-4*9*(x1^2+x2^2+x3^2) + 1]) +QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65] -julia> RM = roadmap(I) -Vector{QQFieldElem}[[], [-3], [-3, -2], [-1], [-1, -2], [-1, -1], [-1, 2], [3], [3, 2]] +julia> RM = roadmap(I, checks=true) +Vector{QQFieldElem}[[], [-3], [-3, -2], [-2], [-2, -1], [-2, 0], [-2, 1], [3], [3, 2]] julia> nb_nodes(RM) 9 julia> all_eqs(RM) 9-element Vector{Ideal{QQMPolyRingElem}}: - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x3 + 4*x2^2*x3 + 4*x3^3 + 4*x3*x4^2 - 40*x3, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 + 3] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 3, x2 + 2] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 + 1] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 + 2] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 + 1] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 + 1, x2 - 2] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, 4*x1^2*x4 + 4*x2^2*x4 + 4*x3^2*x4 + 4*x4^3 + 32*x4, x1 - 3] - QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 64, x1 - 3, x2 - 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, -4*x1^2*x3 - 4*x2^2*x3 - 4*x3^3 - 4*x3*x4^2 + 40*x3, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 + 3] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, x1 + 3, x2 + 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 + 2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, x1 + 2, x2 + 1] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, x1 + 2, x2] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, x1 + 2, x2 - 1] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, -4*x1^2*x4 - 4*x2^2*x4 - 4*x3^2*x4 - 4*x4^3 - 32*x4, x1 - 3] + QQMPolyRingElem[x1^4 + 2*x1^2*x2^2 + 2*x1^2*x3^2 + 2*x1^2*x4^2 - 20*x1^2 + x2^4 + 2*x2^2*x3^2 + 2*x2^2*x4^2 - 20*x2^2 + x3^4 + 2*x3^2*x4^2 - 20*x3^2 + x4^4 + 16*x4^2 + 65, x1 - 3, x2 - 2] ``` """ function roadmap( - I::Ideal{P}; # input ideal - C::Vector{Vector{QQFieldElem}}=Vector{QQFieldElem}[], # query points with rational coefficients - info_level::Int=0, # verbosity level - checks::Bool=false # perform checks (dimension, regularity, etc.) + I::Ideal{P}; # input ideal + C::Vector{Vector{Vector{QQFieldElem}}}=Vector{Vector{QQFieldElem}}[], # query points: interval with rational coefficients + info_level::Int=0, # verbosity level + checks::Bool=false # perform checks (dimension, regularity, etc.) ) where (P <: QQMPolyRingElem) return _roadmap_rec(I, QQFieldElem[], C, info_level, checks) end +@doc Markdown.doc""" + roadmap(I::Ideal{T} where T <: MPolyRingElem, I::Ideal{P}, C::Ideal{P}; info_level::Int=0, checks::Bool=false) +``` +""" function roadmap( I::Ideal{P}, # input ideal C::Ideal{P}; # ideal defining query points @@ -78,7 +82,7 @@ end function _roadmap_rec( I::Ideal{T} where T <: QQMPolyRingElem, # input ideal q::Vector{QQFieldElem}, # single base point with rational coefficients - C::Vector{Vector{QQFieldElem}}, # query points with rational coefficients + C::Vector{Vector{Vector{QQFieldElem}}}, # query points with rational coefficients info_level::Int, # verbosity level checks::Bool # perform checks (dimension, regularity, etc.) ) @@ -144,12 +148,12 @@ function _roadmap_rec( K1WmFq = real_solutions(_fbr(K1WmFq,q), info_level=max(info_level-1,0), nr_thrds=Threads.nthreads(), interval=true) ## New base and query points ## - Cq = isempty(q) ? C : [ c for c in C if c[e] == q[e]] + Cq = isempty(q) ? C : [ c for c in C if q[e] in c[e]] K1W = vcat(K1Fq, K1WmFq) # Heuristic to be proven (Reeb's th) #K1W = K1W[2:end-1] ########## - K1WRat = _mid_rational_points(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) + K1WRat = _mid_rational_points_inter(getindex.(K1W,e+1), unique(getindex.(Cq, e+1))) newQ = vcat.(Ref(q), K1WRat) # Recursively compute roadmap of possible fibers @@ -167,72 +171,14 @@ function _roadmap_rec( end end -function _fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) - @assert(!isempty(I.gens), "Empty polynomial vector") - vars = gens(parent(first(I.gens))) - return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) -end - -function _mid_rational_points(S::Vector{Vector{T}}, Q::Vector{T} = T[]) where {T <: QQFieldElem} - # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] - # such that the [l_i, r_i] are rational and disjoint open intervals. - # * Q is a list of rationals that intersects no [l_i,r_i] - # - # It orders the [l_i,r_i], and compute a list ratioP such that - # strictly between each of these intervals there is: - # - either at least one element of Q - # - or the simplest rational number - isempty(S) && return Q - - S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q) - ratioP = T[] - qidx = 1 - qlen = length(Q1) - - # Handle left gap before first interval - while qidx <= qlen && Q1[qidx] < S1[1][1] - push!(ratioP, Q1[qidx]) - qidx += 1 - end - - # Loop through gaps between sorted disjoint intervals - for i in 1:(length(S1) - 1) - ri, li1 = S1[i][2], S1[i+1][1] - @assert ri < li1 "Intervals are not disjoint." - inserted = false - while qidx <= qlen && Q1[qidx] < li1 - @assert(Q1[qidx] > ri, "A query point is singular") - push!(ratioP, Q1[qidx]) - inserted = true - qidx += 1 - end - if !inserted - eps = (li1 - ri)//1000 # for open interval - # We choose the simplest in absolute value - if -ri > li1 # this means ri is negative and the largest in absolute value - push!(ratioP, -simplest_between(-ri - eps, -li1 + eps)) - else - push!(ratioP, simplest_between(ri + eps, li1 - eps)) - end - end - end - - # Append remaining right-side Q points - while qidx <= qlen - push!(ratioP, Q1[qidx]) - qidx += 1 - end - - return ratioP -end - function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = Vector{T}[]) where {T <: QQFieldElem} # * S is a list of [ [l_1,r_1], ..., [l_n, r_n] ] # such that the [l_i, r_i] are rational and disjoint open intervals. # * Same assumptions on Q # * Intervals in S and Q do not intersect as well # - # It orders the [l_i,r_i], and compute a list ratioP such that + # It orders the [l_i,r_i], and compute a list ratioP + # intersecting all intervals of Q and such that # strictly between each of these intervals there is: # - either at least one element inside an interval of Q # - or the simplest rational number diff --git a/src/types.jl b/src/types.jl index 6ef5453..bd4b7f8 100644 --- a/src/types.jl +++ b/src/types.jl @@ -97,8 +97,14 @@ function _collect_roadmap(RMn::RMnode, F) return data end +function _fbr(I::Ideal{P} where P <: QQMPolyRingElem, Q::Vector{QQFieldElem}) + @assert(!isempty(I.gens), "Empty polynomial vector") + vars = gens(parent(first(I.gens))) + return Ideal(vcat(I.gens, [vars[i] - Q[i] for i in 1:min(length(vars),length(Q))])) +end + function all_eqs(RM::Roadmap) - func(s) = fbr(vcat(RM.initial_ideal.gens, s.polar_eqs), s.base_pt) + func(s) = _fbr(vcat(RM.initial_ideal.gens, s.polar_eqs) |> Ideal, s.base_pt) return _collect_roadmap(RM.root, func) end From 8881239e145ec4a0ac484caf505c8b9f551c9565 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 22 Oct 2025 15:35:30 +0200 Subject: [PATCH 24/33] revert a change --- docs/make.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 8a25b1b..722d8f2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,3 @@ -using Pkg -Pkg.activate("../../AlgebraicSolving.jl") using AlgebraicSolving using Documenter From 872c7154eb13743233d20503ee38ae67e10ad36f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Tue, 4 Nov 2025 17:35:34 +0100 Subject: [PATCH 25/33] start documenting ratparams and roadmaps --- docs/src/types.md | 122 +++++++++++++++++++++++++++++++++- src/algorithms/param-curve.jl | 3 +- 2 files changed, 121 insertions(+), 4 deletions(-) diff --git a/docs/src/types.md b/docs/src/types.md index bc4a2a0..fa440d9 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -41,9 +41,7 @@ R, (x,y,z) = polynomial_ring(GF(101), ["x", "y", "z"], internal_ordering=:degrev ## Ideals -Ideals can be constructed by giving an array of generators. Ideals cache varies -data structures connected to ideals in order to make computational algebra more -effective: +Ideals can be constructed by giving an array of generators: ```@repl using AlgebraicSolving @@ -51,3 +49,121 @@ R, (x,y,z) = polynomial_ring(QQ, ["x", "y", "z"], internal_ordering=:degrevlex) I = Ideal([x+y+1, y*z^2-13*y^2]) ``` +Ideals cache various data structures connected to ideals in order to make +computational algebra more effective. Assume that `T <: MPolyRingElem`, then + * `gens::Vector{T}`: generators of the Ideal provded by the user; + * `deg::Int`: degree of the ideal (`Nothing` if unknown); + * `dim::Int`: Krull dimension of the Ideal (`Nothing` if unknown); + * `gb::Dict{Int, Vector{T}}`: Gröbner bases of the Ideal where the key correspond to the number of variables eliminated (starting with the first); + * `inter_sols::Vector{Vector{Vector{QQFieldElem}}}`: intervals with rational endpoints, containing the real solutions of the system (if `dim=0`); + * `real_sols::Vector{Vector{QQFieldElem}}`: midpoints of the above intervals; + * `rat_sols::Vector{Vector{QQFieldElem}}`: rational solutions of the system (if `dim=0`). + * `rat_param::Union{RationalParametrization, RationalCurveParametrization}`: rational parametrization encoding the complex solution set of the Ideal (see below). +## Rational Parametrizations + +Rational Parametrizations are special data structure for representing equidimensionnal +algebraic sets of dimension 0 and 1. They are typical outputs of algorithm for solving +polynomial systems encoding such algebraic sets thanks to their nice behaviour with +subsequent computations. + +### Zero-dimensionnal parametrizations: finite set of points (see also [there](https://msolve.lip6.fr/downloads/msolve-tutorial.pdf#section.6)) + +A *zero-dimensional parametrization* $\mathscr{P}$ with +coefficients in a field $\mathbb{Q}$ consists of: +* polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t]$ where $t$ is + a new variable and such that + * $w$ is a monic square-free polynomial; + * $w'=1$ when $\mathbb{K}$ is a prime field, $w'=\frac{dw}{dt}$ else; + * $\deg(v_i) < \deg(w)$. +* a linear form $l$ in the variables $x_1,\dotsc, x_n$, +such that +$$ +l(\rho_1,\dotsc,\rho_n) = t \mod w. +$$ + +Such a data-structure encodes the following finite set of points +$$ +\left\{\left (\frac{\rho_1(\vartheta)}{w'(\vartheta)}, +\ldots, \frac{\rho_n(\vartheta)}{w'(\vartheta)} \right) \middle|\, w(\vartheta) = 0\right\}. +$$ +According to this definition, the roots of $w$ are exactly the +values taken by $l$ on this set. + +The type `RationalParametrization` therefore caches the following attributes: + * `vars::Vector{Symbol}`: variables used for the parametrization the last one playing the role of the above $t$ + (hence, maybe with one more variable than the ones given as input); + * `cfs_lf::Vector{ZZRingElem}`: coefficients on the linear form $l$ (when variables are added); + * `elim::QQPolyRingElem`: elimination polynomial $w$ ; + * `denom::QQPolyRingElem`: denominator polynomial (usually $w'$); + * `param::Vector{QQPolyRingElem}`: numerators $\rho_i$'s. + +See the documentation of the [rational_parametrization](@ref) function for further details. + +### One-dimensionnal parametrizations: algebraic curves + +A *one-dimensional parametrization* $\mathscr{C}$ with +coefficients in a field $\mathbb{C}$ consists of: +* polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t,s]$ where $t,s$ are + new variables and such that + * $w$ is a monic square-free polynomial and $w'=\frac{\partial w}{\partial t}$; + * $\deg(v_i) < \deg(w)$. +* two linear forms $(l,l')$ in the variables $x_1,\dotsc, x_n$, +such that +$$ +l(\rho_1,\dotsc,\rho_n) = t\, w'\mod w \quad \text{and} \quad +l'(\rho_1,\dotsc,\rho_n) = s\,w' \mod w +$$ + +Such a data-structure encodes the curve defined as the Zariski closure of the following set +$$ +\left\{\left (\frac{\rho_1(\vartheta,\eta)}{w'(\vartheta,\eta)}, +\ldots, \frac{\rho_n(\vartheta,\eta)}{w'(\vartheta,\eta)} \right) +\middle|\, w(\vartheta,\eta) = 0,\, w'(\vartheta, \eta) \neq 0\right\}. +$$ +According to this definition, the roots of $w$ are exactly the +values taken by $l$ on this set. + +The type `RationalCurveParametrization` therefore caches the following attributes: + * `vars::Vector{Symbol}`: variables used for the parametrization the last two ones playing the role of the above $(t,s)$ + (hence, maybe with up to two more variable than the ones given as input);; + * `cfs_lfs::Vector{Vector{ZZRingElem}}`: coefficients on the linear form $l$ (when variables are added); + * `elim::QQMPolyRingElem`: elimination polynomial $w$ ; + * `denom::QQMPolyRingElem`: denominator polynomial (usually $w'$); + * `param::Vector{QQMPolyRingElem}`: numerators $\rho_i$'s. + +See the documentation of the [rational_curve_parametrization](@ref) function for further details. + +## Roadmaps + +Roadmaps are algebraic curves capturing the connectivity properties of an algebraic +set containing it. They are computed as union of several curves, organized in a tree +due to the recursive nature of the roadmap algorithms. + +```julia +mutable struct Roadmap + initial_ideal::Ideal{QQMPolyRingElem} + root::RMnode +end +``` + +Hence, it is composed of a main node, containing the equations of the initial algebraic set and a reference to the root node of the tree. + +```julia +mutable struct RMnode + polar_eqs::Vector{QQMPolyRingElem} + base_pt::Vector{QQFieldElem} + children::Vector{RMnode} +end +``` + +Each roadmap node (including the root) correspond to a curve component of the roadmap. These are defined in fibers of the initial variety, and are arranged w.r.t. to the specialized variables they share. In particular, all child nodes (referenced in `children`) of a node will have the same specialized variable, plus a new distinct one. + +Each of these components is defined as the polar variety (whose equations are in `polar_eqs`) over a fiber of the initial algebraic set, which are defined by specializing the first variables to the `base_pt` of its parents, and of itsef (in this order). + + +See the documentation of the [roadmap](@ref) function for further details. + + +``` +TODO : un dessin ? les fonctions all_eqs, all_base_pts et nb_nodes (utiliser la doc du code?) +``` \ No newline at end of file diff --git a/src/algorithms/param-curve.jl b/src/algorithms/param-curve.jl index 1e553f2..e1757b4 100644 --- a/src/algorithms/param-curve.jl +++ b/src/algorithms/param-curve.jl @@ -232,7 +232,8 @@ function _evalvar( return LFeval end -# Generate N primes > start that do not divide any numerator/denominator +# Generate N random primes between low and up +# that do not divide any numerator/denominator # of any coefficient in polynomials from LP function _generate_lucky_primes( LF::Vector{P} where P<:MPolyRingElem, From def0143992793cb53e408520d873cd545699b5b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 5 Nov 2025 15:57:10 +0100 Subject: [PATCH 26/33] doc for roadmap --- docs/src/types.md | 77 +++++++++++++++++++++++++++++++++-------------- src/types.jl | 2 +- 2 files changed, 56 insertions(+), 23 deletions(-) diff --git a/docs/src/types.md b/docs/src/types.md index fa440d9..216eb88 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -54,13 +54,19 @@ computational algebra more effective. Assume that `T <: MPolyRingElem`, then * `gens::Vector{T}`: generators of the Ideal provded by the user; * `deg::Int`: degree of the ideal (`Nothing` if unknown); * `dim::Int`: Krull dimension of the Ideal (`Nothing` if unknown); - * `gb::Dict{Int, Vector{T}}`: Gröbner bases of the Ideal where the key correspond to the number of variables eliminated (starting with the first); - * `inter_sols::Vector{Vector{Vector{QQFieldElem}}}`: intervals with rational endpoints, containing the real solutions of the system (if `dim=0`); + * `gb::Dict{Int, Vector{T}}`: Gröbner bases of the Ideal where the key + correspond to the number of variables eliminated (starting with the first); + * `inter_sols::Vector{Vector{Vector{QQFieldElem}}}`: intervals with rational + endpoints, containing the real solutions of the system (if `dim=0`); * `real_sols::Vector{Vector{QQFieldElem}}`: midpoints of the above intervals; - * `rat_sols::Vector{Vector{QQFieldElem}}`: rational solutions of the system (if `dim=0`). - * `rat_param::Union{RationalParametrization, RationalCurveParametrization}`: rational parametrization encoding the complex solution set of the Ideal (see below). -## Rational Parametrizations + * `rat_sols::Vector{Vector{QQFieldElem}}`: rational solutions of the system + (if `dim=0`). + * `rat_param::Union{RationalParametrization, RationalCurveParametrization}`: + rational parametrization encoding the complex solution set of the Ideal (see + below). + +## Rational Parametrizations Rational Parametrizations are special data structure for representing equidimensionnal algebraic sets of dimension 0 and 1. They are typical outputs of algorithm for solving polynomial systems encoding such algebraic sets thanks to their nice behaviour with @@ -123,21 +129,39 @@ $$ According to this definition, the roots of $w$ are exactly the values taken by $l$ on this set. -The type `RationalCurveParametrization` therefore caches the following attributes: - * `vars::Vector{Symbol}`: variables used for the parametrization the last two ones playing the role of the above $(t,s)$ - (hence, maybe with up to two more variable than the ones given as input);; - * `cfs_lfs::Vector{Vector{ZZRingElem}}`: coefficients on the linear form $l$ (when variables are added); +The type `RationalCurveParametrization` therefore caches the following +attributes: + * `vars::Vector{Symbol}`: variables used for the parametrization the last two + ones playing the role of the above $(t,s)$ (hence, maybe with up to two more + variable than the ones given as input);; + * `cfs_lfs::Vector{Vector{ZZRingElem}}`: coefficients on the linear form $l$ + (when variables are added); * `elim::QQMPolyRingElem`: elimination polynomial $w$ ; * `denom::QQMPolyRingElem`: denominator polynomial (usually $w'$); * `param::Vector{QQMPolyRingElem}`: numerators $\rho_i$'s. -See the documentation of the [rational_curve_parametrization](@ref) function for further details. +See the documentation of the [rational_curve_parametrization](@ref) function for +further details. ## Roadmaps -Roadmaps are algebraic curves capturing the connectivity properties of an algebraic -set containing it. They are computed as union of several curves, organized in a tree -due to the recursive nature of the roadmap algorithms. +Consider an algebraic set $V\subset \mathbb{C}^n$ and a finite set of query +points $\mathcal{P} \subset V$, both defined by polynomials with coefficients in +$\mathbb{Q}$. A *roadmap* $\mathcal{R}$ associated to $(V,\mathcal{P})$ is an +algebraic curve such that + * $\mathcal{P} \subset \mathcal{R} \subset V$; + * $C \cap \mathcal{R}$ is non-empty and connected, for each connected + component $C$ of $V\cap \mathbb{R}^n$. + +Roadmaps are algebraic curves capturing the connectivity properties of an +algebraic set containing it so that points in $\mathcal{P}$ are equivalently +connected both on $\mathcal{R}$ and $V\cap \mathbb{R}^n$. + +The effective construction of roadmaps relies on connectivity statements which +allow one to construct real algebraic subsets of $V\cap \mathbb{R}^n$, of +smaller dimension, having a connected intersection with the connected components +of $V\cap \mathbb{R}^n$ Therefore, the output is union of several curves, +organized in a tree due to the recursive nature of the roadmap algorithms. ```julia mutable struct Roadmap @@ -146,24 +170,33 @@ mutable struct Roadmap end ``` -Hence, it is composed of a main node, containing the equations of the initial algebraic set and a reference to the root node of the tree. +Hence, it is composed of a main node, containing the equations of the initial +algebraic set $V$ and a reference to the root node of the tree. ```julia mutable struct RMnode - polar_eqs::Vector{QQMPolyRingElem} base_pt::Vector{QQFieldElem} + polar_eqs::Vector{QQMPolyRingElem} children::Vector{RMnode} end ``` -Each roadmap node (including the root) correspond to a curve component of the roadmap. These are defined in fibers of the initial variety, and are arranged w.r.t. to the specialized variables they share. In particular, all child nodes (referenced in `children`) of a node will have the same specialized variable, plus a new distinct one. +Each roadmap node (including the root) correspond to a curve component of the +roadmap. More precisely, they are defined by as an algebraic subset $W$ (called +polar variety) of a fiber $F$ of the initial variety $V$ such that: + * if `base_pt` contains $\mathbf{q}=(q_1,\dotsc,q_e) \in \mathbb{Q}^e$ then $$ + F = V \cap \left\{x_1=q_1,\,\dotsc,\, x_e=q_e\right\}; $$ + * if `polar_eqs` contains the polynomials $g_1,\dotsc,g_s \in + \mathbb{Q}[x_1,\dotsc,x_n]$ then $$ W = F \cap \{g_1=0,\, \ldots,\, g_s = 0\}. + $$ -Each of these components is defined as the polar variety (whose equations are in `polar_eqs`) over a fiber of the initial algebraic set, which are defined by specializing the first variables to the `base_pt` of its parents, and of itsef (in this order). +Moreover, `children` rassembles all the tree nodes that contains curves +component for which their attribute `base_pt` contains a point $\mathbf{q}'$ +that shares the same first $e$ coordinates with $\mathbf{q}$. +The reason of this tree arrangement is mainly because of the following property: +*two roadmap components intersects if, and only if, one's node is a descendant +of the other's*. In many cases, one can even replace "descendant" by "parent". -See the documentation of the [roadmap](@ref) function for further details. - -``` -TODO : un dessin ? les fonctions all_eqs, all_base_pts et nb_nodes (utiliser la doc du code?) -``` \ No newline at end of file +See the documentation of the [roadmap](@ref) function for further details. \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index bd4b7f8..56f1c08 100644 --- a/src/types.jl +++ b/src/types.jl @@ -79,8 +79,8 @@ Base.getindex(I::Ideal, idx::Union{Int, UnitRange}) = I.gens[idx] Base.lastindex(I::Ideal) = lastindex(I.gens) mutable struct RMnode - polar_eqs::Vector{QQMPolyRingElem} base_pt::Vector{QQFieldElem} + polar_eqs::Vector{QQMPolyRingElem} children::Vector{RMnode} end From 81183244f65a855d45f78e8b04a88e30ab881813 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Wed, 10 Dec 2025 18:44:31 +0100 Subject: [PATCH 27/33] fix args RMnode --- src/connectivity/roadmap.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index a49e079..e0c9a96 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -82,7 +82,7 @@ end function _roadmap_rec( I::Ideal{T} where T <: QQMPolyRingElem, # input ideal q::Vector{QQFieldElem}, # single base point with rational coefficients - C::Vector{Vector{Vector{QQFieldElem}}}, # query points with rational coefficients + C::Vector{Vector{Vector{QQFieldElem}}}, # query points with rational coefficients info_level::Int, # verbosity level checks::Bool # perform checks (dimension, regularity, etc.) ) @@ -111,7 +111,7 @@ function _roadmap_rec( # Terminal case (dim <=1) if I.dim - e <= 1 - return RMnode([], q, RMnode[]) + return RMnode(q, [], RMnode[]) end ## sing(Fq) ## @@ -140,7 +140,7 @@ function _roadmap_rec( else K2Fq.dim = e + 1 end - RM = RMnode(K2Fqmins, q, RMnode[]) + RM = RMnode(q, K2Fqmins, RMnode[]) ## Points with vertical tg in K(pi_2, Fq) ## info_level>0 && println("Compute W-critical points with vertical tangent: K1W") From fd72a6d76dcaf1e68d4a9855202c4e4de043578b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Fri, 27 Feb 2026 16:59:00 +0100 Subject: [PATCH 28/33] typo --- src/connectivity/roadmap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index e0c9a96..318501f 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -210,7 +210,7 @@ function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = qidx += 1 end @assert qidx > qlen || Q1[qidx][1] > S1[i+1][2] "A query point might be singular" - # If there's already rational betwee no need to add new + # If there's already rational between no need to add new !inserted && push!(ratioP, _open_simplest_between(ri, li1, abs(li1 - ri)//1000)) end From 2bc7a32a4d589b8719d65facc3b8d12b032036f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Tue, 3 Mar 2026 12:14:42 +0100 Subject: [PATCH 29/33] fix roadmap doc --- docs/make.jl | 2 ++ docs/src/types.md | 49 ++++++++++++++----------------------- src/connectivity/roadmap.jl | 2 +- 3 files changed, 21 insertions(+), 32 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 722d8f2..563d01f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,3 +1,5 @@ +#using Pkg +#Pkg.activate("../../AlgebraicSolving.jl/") using AlgebraicSolving using Documenter diff --git a/docs/src/types.md b/docs/src/types.md index 216eb88..1a18966 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -78,20 +78,15 @@ A *zero-dimensional parametrization* $\mathscr{P}$ with coefficients in a field $\mathbb{Q}$ consists of: * polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t]$ where $t$ is a new variable and such that - * $w$ is a monic square-free polynomial; - * $w'=1$ when $\mathbb{K}$ is a prime field, $w'=\frac{dw}{dt}$ else; - * $\deg(v_i) < \deg(w)$. -* a linear form $l$ in the variables $x_1,\dotsc, x_n$, + * $w$ is a monic square-free polynomial; + * $w'=1$ when $\mathbb{K}$ is a prime field, $w'=\frac{dw}{dt}$ else; + * $\deg(v_i) < \deg(w)$. +* a linear form $l$ in the variables $x_1,\dotsc, x_n$, such that -$$ -l(\rho_1,\dotsc,\rho_n) = t \mod w. -$$ +$l(\rho_1,\dotsc,\rho_n) = t \mod w$. Such a data-structure encodes the following finite set of points -$$ -\left\{\left (\frac{\rho_1(\vartheta)}{w'(\vartheta)}, -\ldots, \frac{\rho_n(\vartheta)}{w'(\vartheta)} \right) \middle|\, w(\vartheta) = 0\right\}. -$$ +* $\left\{\left(\frac{\rho_1(\vartheta)}{w'(\vartheta)}, \ldots, \frac{\rho_n(\vartheta)}{w'(\vartheta)} \right) \middle|\, w(\vartheta) = 0\right\}$. According to this definition, the roots of $w$ are exactly the values taken by $l$ on this set. @@ -103,7 +98,7 @@ The type `RationalParametrization` therefore caches the following attributes: * `denom::QQPolyRingElem`: denominator polynomial (usually $w'$); * `param::Vector{QQPolyRingElem}`: numerators $\rho_i$'s. -See the documentation of the [rational_parametrization](@ref) function for further details. +See the documentation of the [`rational_parametrization`](@ref) function for further details. ### One-dimensionnal parametrizations: algebraic curves @@ -111,21 +106,14 @@ A *one-dimensional parametrization* $\mathscr{C}$ with coefficients in a field $\mathbb{C}$ consists of: * polynomials $(w, w', v_1, \ldots, v_n)$ in $\mathbb{K}[t,s]$ where $t,s$ are new variables and such that - * $w$ is a monic square-free polynomial and $w'=\frac{\partial w}{\partial t}$; - * $\deg(v_i) < \deg(w)$. + * $w$ is a monic square-free polynomial and $w'=\frac{\partial w}{\partial t}$; + * $\deg(v_i) < \deg(w)$. * two linear forms $(l,l')$ in the variables $x_1,\dotsc, x_n$, such that -$$ -l(\rho_1,\dotsc,\rho_n) = t\, w'\mod w \quad \text{and} \quad -l'(\rho_1,\dotsc,\rho_n) = s\,w' \mod w -$$ +$l(\rho_1,\dotsc,\rho_n) = t\, w'\mod w \quad \text{and} \quad l'(\rho_1,\dotsc,\rho_n) = s\,w' \mod w$ Such a data-structure encodes the curve defined as the Zariski closure of the following set -$$ -\left\{\left (\frac{\rho_1(\vartheta,\eta)}{w'(\vartheta,\eta)}, -\ldots, \frac{\rho_n(\vartheta,\eta)}{w'(\vartheta,\eta)} \right) -\middle|\, w(\vartheta,\eta) = 0,\, w'(\vartheta, \eta) \neq 0\right\}. -$$ +* $\left\{\left(\frac{\rho_1(\vartheta,\eta)}{w'(\vartheta,\eta)}, \ldots, \frac{\rho_n(\vartheta,\eta)}{w'(\vartheta,\eta)} \right) \middle|\, w(\vartheta,\eta) = 0,\, w'(\vartheta, \eta) \neq 0\right\}$. According to this definition, the roots of $w$ are exactly the values taken by $l$ on this set. @@ -140,7 +128,7 @@ attributes: * `denom::QQMPolyRingElem`: denominator polynomial (usually $w'$); * `param::Vector{QQMPolyRingElem}`: numerators $\rho_i$'s. -See the documentation of the [rational_curve_parametrization](@ref) function for +See the documentation of the [`rational_curve_parametrization`](@ref) function for further details. ## Roadmaps @@ -149,8 +137,8 @@ Consider an algebraic set $V\subset \mathbb{C}^n$ and a finite set of query points $\mathcal{P} \subset V$, both defined by polynomials with coefficients in $\mathbb{Q}$. A *roadmap* $\mathcal{R}$ associated to $(V,\mathcal{P})$ is an algebraic curve such that - * $\mathcal{P} \subset \mathcal{R} \subset V$; - * $C \cap \mathcal{R}$ is non-empty and connected, for each connected + * $\mathcal{P} \subset \mathcal{R} \subset V$; + * $C \cap \mathcal{R}$ is non-empty and connected, for each connected component $C$ of $V\cap \mathbb{R}^n$. Roadmaps are algebraic curves capturing the connectivity properties of an @@ -184,11 +172,10 @@ end Each roadmap node (including the root) correspond to a curve component of the roadmap. More precisely, they are defined by as an algebraic subset $W$ (called polar variety) of a fiber $F$ of the initial variety $V$ such that: - * if `base_pt` contains $\mathbf{q}=(q_1,\dotsc,q_e) \in \mathbb{Q}^e$ then $$ - F = V \cap \left\{x_1=q_1,\,\dotsc,\, x_e=q_e\right\}; $$ - * if `polar_eqs` contains the polynomials $g_1,\dotsc,g_s \in - \mathbb{Q}[x_1,\dotsc,x_n]$ then $$ W = F \cap \{g_1=0,\, \ldots,\, g_s = 0\}. - $$ + * if `base_pt` contains $\mathbf{q}=(q_1,\dotsc,q_e) \in \mathbb{Q}^e$ then + * $F = V \cap \left\{x_1=q_1,\,\dotsc,\, x_e=q_e\right\}$; + * if `polar_eqs` contains the polynomials $g_1,\dotsc,g_s \in \mathbb{Q}[x_1,\dotsc,x_n]$ then + * $W = F \cap \{g_1=0,\, \ldots,\, g_s = 0\}$. Moreover, `children` rassembles all the tree nodes that contains curves component for which their attribute `base_pt` contains a point $\mathbf{q}'$ diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 318501f..0f1e6e8 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -179,7 +179,7 @@ function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = # # It orders the [l_i,r_i], and compute a list ratioP # intersecting all intervals of Q and such that - # strictly between each of these intervals there is: + # strictly between each of the [l_i,r_i] there is: # - either at least one element inside an interval of Q # - or the simplest rational number # TODO: From 8d4abfebc519a518ef6331cc9778a62e39ca1bcb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20P=2E?= Date: Tue, 3 Mar 2026 14:17:51 +0100 Subject: [PATCH 30/33] exported functions --- src/connectivity/roadmap.jl | 13 +------------ src/exports.jl | 3 ++- 2 files changed, 3 insertions(+), 13 deletions(-) diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index 0f1e6e8..ae3da33 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -1,12 +1,3 @@ -#using Pkg -#Pkg.activate("AlgebraicSolving.jl") -#using Revise -#using AlgebraicSolving -#using Nemo - -export roadmap, computepolar, _mid_rational_points, _mid_rational_points_inter, _fbr, all_eqs, all_base_pts, nb_nodes, computepolar -include("polar.jl") - @doc Markdown.doc""" roadmap(I::Ideal{T} where T <: MPolyRingElem, ) @@ -179,10 +170,9 @@ function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = # # It orders the [l_i,r_i], and compute a list ratioP # intersecting all intervals of Q and such that - # strictly between each of the [l_i,r_i] there is: + # *strictly* between each of the [l_i,r_i] there is: # - either at least one element inside an interval of Q # - or the simplest rational number - # TODO: isempty(S) && return Q S1, Q1 = sort(S, lt=(x, y) -> x[2] <= y[1]), sort(Q, lt=(x, y) -> x[2] <= y[1]) @@ -224,7 +214,6 @@ function _mid_rational_points_inter(S::Vector{Vector{T}}, Q::Vector{Vector{T}} = return ratioP end - function _open_simplest_between(a::QQFieldElem, b::QQFieldElem, eps::QQFieldElem) # We choose the simplest in absolute value if -a > b # this means a is negative and the largest in absolute value diff --git a/src/exports.jl b/src/exports.jl index 5e3d57e..ce4142a 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -5,4 +5,5 @@ export polynomial_ring, MPolyRing, GFElem, MPolyRingElem, prime_field, sig_groebner_basis, cyclic, leading_coefficient, equidimensional_decomposition, homogenize, dimension, FqMPolyRingElem, hilbert_series, hilbert_dimension, hilbert_degree, hilbert_polynomial, - rational_curve_parametrization + rational_curve_parametrization, roadmap, computepolar, all_eqs, + all_base_pts, nb_nodes From bc8d1af1abe90f30e5ceb4ab4598f11c7fe14847 Mon Sep 17 00:00:00 2001 From: Mohab Safey El Din Date: Mon, 23 Mar 2026 14:42:20 +0100 Subject: [PATCH 31/33] minor typo in the doc of roadmap --- src/connectivity/roadmap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl index ae3da33..93c5a60 100644 --- a/src/connectivity/roadmap.jl +++ b/src/connectivity/roadmap.jl @@ -6,7 +6,7 @@ whose real trace XR is bounded, return a roadmap of XR The output is given as a Roadmap structure, encoding the recursive structure of roadmaps. It is encoded as a chained list, whose root is containing the equations defining X -and each node representing a curve component, that is defined by additional polar equation and base point< +and each node representing a curve component, that is defined by additional polar equation and base point. Moreover it is linked to fibers, that share the same base point. # Arguments From e7d0fef523ad43277a2528caee186f511c343858 Mon Sep 17 00:00:00 2001 From: Mohab Safey El Din Date: Mon, 23 Mar 2026 14:57:06 +0100 Subject: [PATCH 32/33] includes polar.jl in AlgebraicSolving.jl --- src/AlgebraicSolving.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index e3fa4ff..053771f 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -19,6 +19,7 @@ include("algorithms/hilbert.jl") #= siggb =# include("siggb/siggb.jl") #= connectivity =# +include("connectivity/polar.jl") include("connectivity/roadmap.jl") #= examples =# include("examples/katsura.jl") From 8b35b5797e23883c9ef5a762c24394a09e9c1269 Mon Sep 17 00:00:00 2001 From: Mohab Safey El Din Date: Tue, 24 Mar 2026 10:39:06 +0100 Subject: [PATCH 33/33] fix typo in documentation --- docs/src/types.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/types.md b/docs/src/types.md index 1a18966..a0d0d08 100644 --- a/docs/src/types.md +++ b/docs/src/types.md @@ -182,8 +182,8 @@ component for which their attribute `base_pt` contains a point $\mathbf{q}'$ that shares the same first $e$ coordinates with $\mathbf{q}$. The reason of this tree arrangement is mainly because of the following property: -*two roadmap components intersects if, and only if, one's node is a descendant +*two roadmap components intersect if, and only if, one's node is a descendant of the other's*. In many cases, one can even replace "descendant" by "parent". -See the documentation of the [roadmap](@ref) function for further details. \ No newline at end of file +See the documentation of the [roadmap](@ref) function for further details.