diff --git a/docs/make.jl b/docs/make.jl index 9c761f2..563d01f 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..a0d0d08 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,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,141 @@ 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 + +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 + initial_ideal::Ideal{QQMPolyRingElem} + root::RMnode +end +``` + +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 + 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. 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\}$. + +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 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. diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index bec57f4..053771f 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -18,6 +18,9 @@ include("algorithms/param-curve.jl") include("algorithms/hilbert.jl") #= siggb =# include("siggb/siggb.jl") +#= connectivity =# +include("connectivity/polar.jl") +include("connectivity/roadmap.jl") #= examples =# include("examples/katsura.jl") include("examples/cyclic.jl") 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, diff --git a/src/algorithms/solvers.jl b/src/algorithms/solvers.jl index 05805a2..7646ee5 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/polar.jl b/src/connectivity/polar.jl new file mode 100644 index 0000000..ad0805c --- /dev/null +++ b/src/connectivity/polar.jl @@ -0,0 +1,112 @@ +@doc Markdown.doc""" + 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 +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) + 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 + 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( + 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) + + R = parent(V) + n = nvars(R) + nphi = length(phi) + @assert all([1<=j<=n for j in J]) + + isnothing(V.dim) && dimension(V) + c = n - V.dim + + ## + 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 = 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) + + if only_mins + return minors + else + return vcat(V.gens, minors) + end +end + +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(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 _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 +end \ No newline at end of file diff --git a/src/connectivity/roadmap.jl b/src/connectivity/roadmap.jl new file mode 100644 index 0000000..93c5a60 --- /dev/null +++ b/src/connectivity/roadmap.jl @@ -0,0 +1,226 @@ +@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. +- `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 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) + 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, 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 + 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{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 + 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{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 + return [I.gens] + end + + # 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 + 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) + if I.dim - e <= 1 + return RMnode(q, [], RMnode[]) + end + + ## sing(Fq) ## + if checks + 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("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("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 + 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(q, K2Fqmins, RMnode[]) + + ## 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) + + ## New base and query points ## + 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_inter(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_rec(I, newq, Cq, info_level, checks) + push!(RM.children, RMFq) + end + end + + if e == 0 + return Roadmap(I, RM) + else + return RM + end +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 + # intersecting all intervals of Q and such that + # *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 + 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 between 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 + + 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 diff --git a/src/types.jl b/src/types.jl index 20896d6..56f1c08 100644 --- a/src/types.jl +++ b/src/types.jl @@ -74,4 +74,49 @@ 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) + +mutable struct RMnode + base_pt::Vector{QQFieldElem} + polar_eqs::Vector{QQMPolyRingElem} + 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 _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) |> Ideal, 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