Real Algebraic Geometry

points_per_components(eqs::Vector{QQMPolyRingElem}, pos::Vector{QQMPolyRingElem}, ineqs::Vector{QQMPolyRingElem}; nr_thrds::Int=1, worker_pool::AbstractWorkerPool=default_worker_pool(), info_level::Int=0)::Vector{Vector{Vector{QQFieldElem}}}

Given a semi-algebraic set defined by equations eqs, positivity constraints pos, and non-vanishing constraints ineqs, returns a vector of points meeting all connected components of the set. Each point is represented by an isolating box, i.e., a vector of intervals represented by pairs of rational numbers.

NOTE: Sampling points per components is currently only implemented for systems without equations, i.e., eqs must be empty.

Arguments

• eqs::Vector{QQMPolyRingElem}: equations defining the semi-algebraic set.
• pos::Vector{QQMPolyRingElem}: positivity constraints defining the semi-algebraic set.
• ineqs::Vector{QQMPolyRingElem}: non-vanishing constraints defining the semi-algebraic set.
• nr_thrds::Int=1: the number of threads to use for parallel computations.
• worker_pool::AbstractWorkerPool=default_worker_pool(): the worker pool to use for parallel computations.
• info_level::Int=0: info level printout: off (0, default), summary (1), detailed (2).

Examples

julia> using AlgebraicSolving

julia> R, (x, y) = polynomial_ring(QQ, ["x", "y"])
(Multivariate polynomial ring in 2 variables over QQ, QQMPolyRingElem[x, y])

julia> length(points_per_components(QQMPolyRingElem[], [x, y], [x + y - 1])) >= 2
true

real_root_classification(I::ParametricIdeal{K}, g::Vector{MPoly{K}}; <keyword arguments>) -> SemialgebraicSet

Computes the real root classification of the polynomials g with respect to the parametric ideal I. The result is a semi-algebraic set that describes the regions in the parameter space where the number of real roots such that g is positive is constant. Each region is accompanied by the constant number of real roots in that region, and a witness point in that region.

Note: The regions described by the semi-algebraic set may differ from the actual regions due to the fact that the Hermite matrices may not specialize well at certain parameter values. However, their symmetric difference is guaranteed to be contained in a proper algebraic set.

Arguments

• I::ParametricIdeal{K}: input parametric ideal.
• g::Vector{MPoly{K}}: input list of polynomials.
• nr_thrds::Int=1: the number of threads to use for parallel computations.
• worker_pool::AbstractWorkerPool=default_worker_pool(): the worker pool to use for parallel computations.
• info_level::Int=0: info level printout: off (0, default), summary (1), detailed (2).
• show_progress::Bool=info_level >= 1: whether to show progress bars during computations.
• output_form::Symbol=:sign: the output form of the real root classification, can be either polynomials accompanied by their signs (:sign), Hermite matrices accompanied by their signatures (:signature), or Hermite matrices and multiplication matrices (:matrix).
• ignore_no_real_roots::Bool=false: whether to ignore sign conditions with no real roots.

hermite_matrix(gb::Vector{T}, b::Vector{T}, g::T) -> MatElem

Computes the Hermite matrix associated with the polynomial g with respect to the Groebner basis gb and the monomial basis b.

hermite_matrix(I::ParametricIdeal{K}, g::MPoly{K}; <keyword arguments>) -> MatSpaceElem{K}

Computes the parametric Hermite matrix associated with the polynomial g with respect to the parametric ideal I. The Hermite matrix is computed by evaluating the ideal at a number of parameter values and interpolating the coefficients of the resulting Hermite matrices.

Arguments

• retry::Int=10: the maximum number of consecutive failures allowed when computing the Hermite matrix or interpolating the coefficients.
• nr_thrds::Int=1: the number of threads to use for parallel computations.
• worker_pool::AbstractWorkerPool=default_worker_pool(): the worker pool to use for parallel computations.
• show_progress::Bool=false: whether to show a progress bar while computing the Hermite matrix.

hermite_matrix(g::Vector{MPoly{K}}, α::Vector{Int}, I::ParametricIdeal{K}; <keyword arguments>) -> MatSpaceElem{K}

Computes the parametric Hermite matrix associated with the polynomials g raised to the powers specified in α with respect to the parametric ideal I. The Hermite matrix is computed by evaluating the ideal at a number of parameter values and interpolating the coefficients of the resulting Hermite matrices.

Arguments

• nr_thrds::Int=1: the number of threads to use for parallel computations.
• worker_pool::AbstractWorkerPool=default_worker_pool(): the worker pool to use for parallel computations.
• show_progress::Bool=false: whether to show a progress bar while computing the Hermite matrix.

hermite_matrix(g::Vector{MPoly{K}}, α::Vector{Int}, I::ParametricIdeal{K}, vals::Vector{QQFieldElem}; <keyword arguments>) -> QQMatrix

Computes the Hermite matrix associated with the polynomials g raised to the powers specified in α with respect to the parametric ideal I evaluated at the parameter values specified in vals. The Hermite matrix is computed by evaluating the ideal at the specified parameter values and interpolating the coefficients of the resulting Hermite matrices.

Arguments

• nr_thrds::Int=1: the number of threads to use for parallel computations.
• worker_pool::AbstractWorkerPool=default_worker_pool(): the worker pool to use for parallel computations.
• show_progress::Bool=false: whether to show a progress bar while computing the Hermite matrix.

roadmap(I::Ideal{T} where T <: MPolyRingElem, <keyword arguments>)

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

)

Examples

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)
Vector{Tuple{QQMPolyRingElem, QQFieldElem}}[[], [(x1, -3)], [(x1, -3), (x2, -2)], [(x1, -2)], [(x1, -2), (x2, -1)], [(x1, -2), (x2, 0)], [(x1, -2), (x2, 1)], [(x1, 3)], [(x1, 3), (x2, 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, 0, 0]
 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]

computepolar(J::Union{Vector{Int},UnitRange{Int}}, V::Ideal{P} where P <: MPolyRingElem; <keyword arguments>)

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

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

curve_graph(I::Ideal; generic=true, outf=true, kwargs...)
curve_graph(P::CurveRationalParametrization; kwargs...)

# Both functions accept optional control points C in various formats:
curve_graph(I, C::Vector{P}; kwargs...)
curve_graph(I, C::Vector{Vector{P}}; kwargs...)
curve_graph(I, C::Dict{Int, Vector{P}}; kwargs...)

Computes a planar straight-line graph that is homeomorphic to a real algebraic (space) curve.

Core Workflow & Pre-processing

1. Parametrization: If given an Ideal, it first computes a CurveRationalParametrization. If generic!=0, it applies a random linear transformation to place the curve into generic position. The absolute size of the used coefficients is generic if positive, 100 else.
2. Coefficient Extraction: Pulls the planar projection f(x,y) = 0 and the vertical lift z = g(x,y) / df/dy(x,y) from the parametrization.
3. Graph Construction: Computes bounding boxes for critical points, routes connections, and identifies "apparent singularities" (2D crossings that do not intersect in 3D).

Output Data Structure

Returns a CurveGraph{T} object (where T is determined by the outf flag), containing:

• Vert::Vector{Tuple{T, T}}: 2D coordinates of the graph vertices (critical points, routing nodes, control points).
• Edg::Vector{Tuple{Int, Int}}: Index pairs defining undirected edges between Vert indices.
• Vcon::Dict{Int, Vector{Int}}: A mapping from the original index of a control point keys its vertex index in Vert.

Arguments

• I (Ideal): The algebraic ideal defining the curve.
• P (CurveRationalParametrization): A pre-computed rational parametrization.
• C (Optional): User-defined plane points on the curve (control points). Given either as Ideal or RationalParametrization.
• generic (Union{ZZRingElem, Int64}, default 0): If !=0, applies a random shear transformation.
• precx (Int, default 150): Base numerical precision for real root isolation.
• v (Int, default 0): Verbosity level.
• force_app (Bool, default false): Skips 3D intersection checks, treating all 2D nodes as apparent singularities.
• outf (Bool, default true): Output coordinates as Float64. If false, outputs exact QQFieldElem.

Example

julia> using AlgebraicSolving

julia> R, (x, y, z) = polynomial_ring(QQ, ["x", "y", "z"]);

julia> f = -42*x^2 + 101*x*y - 8*x*z - 88*x + 53*y^2 + 71*y*z + 154*y + 53*z^2 + 2*z - 32;
       g = -7*x^2 - 144*x*y - 126*x*z - 7*x + 30*y^2 - 34*y*z + 72*y + 71*z^2 + 46*z + 5;

julia> I = Ideal([f, g]);

julia> G = curve_graph(I);

julia> number_of_connected_components(G)
3

curve_arrangement_graph(curves::Vector{Ideal}; generic::Union{ZZRingElem, Int64}=0, outf=true, v=0, kwargs...)

Computes the combined planar graph of an arrangement of multiple space curves in complete intersection. Automatically computes the mutual intersections between all curves and guarantees they are projected using a shared, unified linear form.

connected_components(G::CurveGraph{T}) where T

Decomposes a graph into its topologically connected components using a Depth-First Search (DFS).

This function is particularly useful when a single algebraic curve resolves into multiple disconnected real branches in space. It guarantees that isolated vertices are preserved as their own zero-edge components.

Crucially, this function preserves the control_nodes mappings. If a global control node belonged to a specific vertex, it is correctly remapped to that vertex's new local index within its specific disconnected subgraph.

Inputs

• G::CurveGraph{T}: The unified curve graph to be decomposed.

Outputs

• Vector{CurveGraph{T}}: A list of disjoint subgraphs. Each subgraph retains the mapped coordinates, edges, and control nodes belonging to its component.

group_by_component(G::CurveGraph)

Extracts the control_nodes dictionaries partitioned by the graph's connected components. This is primarily used to determine which external connections (or user-defined points) belong to which isolated branch of the curve.

Inputs

• G::CurveGraph: The parent graph structure.

Outputs

• Vector{Dict{Int, Vector{Int}}}: A list of dictionaries, one for each connected component,

detailing the local vertex indices mapped to specific control IDs.

number_of_connected_components(G::CurveGraph)

Computes the total number of topologically disconnected branches/components in the graph. This is a convenience wrapper for topological validation.

Inputs

• G::CurveGraph: The parent graph structure.

Outputs

• Int: The count of isolated subgraphs.

merge_graphs(graphs::Vector{CurveGraph{T}}) where T

Merges a collection of disjoint graphs into a single graph, mathematically fusing vertices based on their control_nodes dictionnaries that describe a one-to-one correspondance between subset of vertices of each graph.

Fusion Mechanism:

This work as follows. Given V1 = [u1,u2] and V2 = [w1,w2,w3] the vertices of two graphs G1 and G2. Let C1= {2 => [1]} and C2 = {1 => [2]} be their respective control_nodes disctionnary. Then it merges the first vertex w1 of G2 with the second vertex u2 of G1. The merged graph has vertices [u1,u2,w2,w3].

Note that control nodes that refer to external bounds (e.g., indices k <= 0 or k > length(graphs)) are preserved in the final merged graph to allow for subsequent higher-level merging or tracking.

Inputs

• graphs::Vector{CurveGraph{T}}: A list of graphs. graphs[i].control_nodes[k] should contain the local vertex indices in graph i that connect to graph k. This vertex correspondance must be one-to-one.

Outputs

• CurveGraph{T}: A unified graph struct.

build_graph_data(G::CurveGraph{T}; width=3.0, vert=true, color="#FF0000") where T <: Union{Float64,QQFieldElem}

Construct plotting data from a CurveGraph.

This function converts the piecewise linear structure encoded in G into a GraphPlotData object, separating edges and point groups (vertices and control points), ready for visualization with plotting libraries such as Plots.jl (see documentation of GraphPlotData)

Arguments

• G::CurveGraph{T}: Input graph containing vertices, control nodes, and edges.
• width::Real=3.0: Line width used for rendering edges.
• vert::Bool=true: If true, include vertices in the graph.
• color::AbstractString="#FF0000": Color used for both edges and points.

Behavior

• Edges are converted into line segments between vertex coordinates.
• Vertices (if enabled) are displayed with marker :x.
• Control nodes are grouped and displayed with marker :+.

Returns

• GraphPlotData: A structure containing:
  • one EdgeGroup for all edges
  • multiple PointGroups for vertices and control nodes

Notes

• All coordinates are converted to Float64 for compatibility with plotting backends.

build_graph_data(CG::Vector{CurveGraph{T}}; width=3.0, vert=true) where T <: Union{Float64,QQFieldElem}

Construct plotting data for multiple CurveGraphs with automatically assigned colors.

Arguments

• CG::Vector{CurveGraph{T}}: Collection of curve graphs to visualize.
• width::Real=3.0: Line width used for rendering edges.
• vert::Bool=true: If true, include vertices in each graph.

Behavior

• Colors are generated using distinguishable_colors to ensure visual contrast.
• Each graph is assigned a unique color consistently across its edges and points.

Returns

• Vector{GraphPlotData}: One plotting data object per input graph.

Notes

• The number of colors scales with length(CG).