diff --git a/src/AlgebraicSolving.jl b/src/AlgebraicSolving.jl index 846f3767..94f40465 100644 --- a/src/AlgebraicSolving.jl +++ b/src/AlgebraicSolving.jl @@ -17,6 +17,11 @@ include("algorithms/decomposition.jl") include("algorithms/param-curve.jl") include("algorithms/hilbert.jl") include("algorithms/sampling.jl") +#= param-ideal =# +include("algorithms/param-ideal/groebner-bases.jl") +include("algorithms/param-ideal/multiplication-matrices.jl") +include("algorithms/param-ideal/hermite-matrices.jl") +include("algorithms/param-ideal/parametrizations.jl") #= siggb =# include("siggb/siggb.jl") #= progress =# diff --git a/src/algorithms/param-ideal/groebner-bases.jl b/src/algorithms/param-ideal/groebner-bases.jl new file mode 100644 index 00000000..38d6f3fb --- /dev/null +++ b/src/algorithms/param-ideal/groebner-bases.jl @@ -0,0 +1,229 @@ +function _map_exponent_vectors(f, p::MPolyRingElem, R::MPolyRing)::MPolyRingElem + cvzip = zip(coefficients(p), exponent_vectors(p)) + M = MPolyBuildCtx(R) + for (c, v) in cvzip + push_term!(M, coefficient_ring(R)(c), f(v)) + end + finish(M) +end + +function _change_ring(p::MPolyRingElem, R::MPolyRing)::MPolyRingElem + n = number_of_variables(parent(p)) + n′ = number_of_variables(R) + if n == 0 || n′ == 0 + return R(constant_coefficient(p)) + end + indices = Vector{Int}() + for s in symbols(R) + i = findfirst(==(s), symbols(parent(p))) + if isnothing(i) + push!(indices, 0) + else + push!(indices, i) + end + end + return _map_exponent_vectors(v -> [i == 0 ? 0 : v[i] for i in indices], p, R) +end + +function _change_ring(v::Vector{<:MPolyRingElem}, R::MPolyRing)::Vector{<:MPolyRingElem} + map(p -> _change_ring(p, R), v) +end + +function _change_ring(A::MatSpaceElem{<:MPolyRingElem}, R::MPolyRing)::MatSpaceElem{<:MPolyRingElem} + map(f -> _change_ring(f, R), A) +end + +function _saturate( + f::Vector{QQMPolyRingElem}, + g::QQMPolyRingElem; + worker_pool::AbstractWorkerPool=default_worker_pool() +)::Vector{QQMPolyRingElem} + @assert !isempty(f) + R = parent(f[1]) + n = number_of_variables(R) + R′, vars = polynomial_ring(QQ, [:sat; symbols(R)]) + f′ = _change_ring(f, R′) + g′ = _change_ring(g, R′) + gb = AlgebraicSolving.groebner_basis( + AlgebraicSolving.Ideal([f′; vars[1] * g′ - 1]); + eliminate=1, + worker_pool=worker_pool + ) + gb = _change_ring(gb, R) + gb +end + +function _radical( + f::Vector{QQMPolyRingElem}; + worker_pool::AbstractWorkerPool=default_worker_pool() +)::Vector{QQMPolyRingElem} + @assert !isempty(f) + R = parent(f[1]) + n = number_of_variables(R) + s = symbols(R) + rads = Vector{QQMPolyRingElem}() + for _ in 1:n + R′, _ = polynomial_ring(QQ, s, internal_ordering=internal_ordering(R)) + disc = AlgebraicSolving.groebner_basis( + AlgebraicSolving.Ideal(_change_ring(f, R′)); + eliminate=n - 1, + worker_pool=worker_pool + ) + r = R′(1) + for (p, _) in factor(prod(disc)) + r *= p + end + push!(rads, _change_ring(r, R)) + s = [s[2:end]; s[1]] + end + gb = AlgebraicSolving.groebner_basis( + AlgebraicSolving.Ideal([f; rads]); + worker_pool=worker_pool + ) + gb = _change_ring(gb, R) + gb +end + +function groebner_basis( + I::ParametricIdeal{K}, + vals::Vector{QQFieldElem}; + worker_pool::AbstractWorkerPool=default_worker_pool() +)::Vector{QQMPolyRingElem} where {K<:FracFieldElem} + gb_spec_hash = hash((I.gens, vals)) + if haskey(I.gb_spec, gb_spec_hash) + return I.gb_spec[gb_spec_hash] + end + n = I.num_vars + elim = number_of_variables(parent(I.gens[1])) - n + R = parent(I.gens[1]) + R₀, _ = polynomial_ring(QQ, symbols(R)[1:end], internal_ordering=:degrevlex) + R₁, _ = polynomial_ring(QQ, I.symbols, internal_ordering=:degrevlex) + fₓ = _change_ring(map(p -> map_coefficients(k -> evaluate(k, vals), p), I.gens), R₀) + gb = AlgebraicSolving.groebner_basis( + AlgebraicSolving.Ideal(fₓ); + eliminate=elim, + worker_pool=worker_pool + ) + gb = _change_ring(gb, R₁) + for sat in I.sats + gb = _saturate(gb, _change_ring(map_coefficients(k -> evaluate(k, vals), sat), R₁); worker_pool=worker_pool) + end + if I.radicalize + gb = _radical(gb; worker_pool=worker_pool) + end + gb = map(p -> p / leading_coefficient(p), gb) + lock(I.gb_spec_lock) do + I.gb_spec[gb_spec_hash] = gb + end + gb +end + +function _is_zero_dimensional(gb::Vector{T})::Bool where {T<:MPolyRingElem} + if isempty(gb) + return false + end + R = parent(gb[1]) + n = number_of_variables(R) + has_power = falses(n) + for f in gb + indices = var_indices(Nemo.leading_monomial(f)) + if isempty(indices) + return true + end + if length(indices) == 1 + has_power[indices[1]] = true + end + end + all(has_power) +end + +function _monomial_basis(gb::Vector{T})::Vector{T} where {T<:MPolyRingElem} + if !_is_zero_dimensional(gb) + error("The ideal is not zero-dimensional, so a finite monomial basis does not exist.") + end + R = parent(gb[1]) + n = number_of_variables(R) + lms = [Nemo.leading_monomial(f) for f in gb] + basis = typeof(gb[1])[] + queue = [one(R)] + while !isempty(queue) + u = popfirst!(queue) + if all(!divides(u, lm)[1] for lm in lms) && !(u in basis) + push!(basis, u) + for i in n:-1:1 + push!(queue, u * gens(R)[i]) + end + end + end + basis +end + +function _expected_degree( + I::ParametricIdeal{K}; + worker_pool::AbstractWorkerPool=default_worker_pool() +)::Int where {K<:FracFieldElem} + gb = groebner_basis(I, QQ.(I.shift); worker_pool=worker_pool) + b = _monomial_basis(gb) + length(b) +end + +@doc Markdown.doc""" + groebner_basis(I::ParametricIdeal{K}; ) -> Vector{Generic.MPoly{FracFieldElem{QQMPolyRingElem}}} + +Computes a parametric Gröbner basis of the parametric ideal `I` w.r.t. the degree reverse lexicographic ordering. The Gröbner basis is computed by evaluating the ideal at a number of parameter values and interpolating the coefficients of the resulting Gröbner bases. + +# Arguments +- `retry::Int=10`: the maximum number of consecutive failures allowed when computing the Gröbner basis 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 Gröbner basis. +""" +function groebner_basis( + I::ParametricIdeal{K}; + retry::Int=10, + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::Vector{Generic.MPoly{FracFieldElem{QQMPolyRingElem}}} where {K<:FracFieldElem} + input_hash = hash(I.gens) + if haskey(I.gb, input_hash) + return I.gb[input_hash] + end + R = parent(I.gens[1]) + R₂ = base_ring(I.base_frac_field) + R₄, _ = polynomial_ring(I.base_frac_field, I.symbols, internal_ordering=internal_ordering(R)) + gb_shift_mons = nothing + function bb(v) + gb = groebner_basis(I, v; worker_pool=worker_pool) + if !isnothing(gb_shift_mons) + gb_mons = [collect(monomials(p)) for p in gb] + if gb_mons != gb_shift_mons + error("Error evaluating black box function: the Gröbner basis has unexpected monomial structure at parameter values $(v).") + end + end + gb + end + gb_shift = bb(QQ.(I.shift)) + gb_shift_mons = [collect(monomials(p)) for p in gb_shift] + pgb = fill(R₄(0), length(gb_shift)) + len = length(gb_shift) + len_mons = length.(gb_shift) + for i in 1:len + for j in 1:len_mons[i] + denom = R₂(1) + c = Interpolation.cuyt_lee( + R₂, v -> evaluate(denom, v) * coeff(bb(v)[i], j); + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric Gröbner basis, size=($len,$(maximum(len_mons))), index=($i,$j)" + ) + pgb[i] += R₄(c) / denom * _change_ring(Nemo.monomial(gb_shift[i], j), R₄) + end + end + lock(I.gb_lock) do + I.gb[input_hash] = pgb + end + pgb +end diff --git a/src/algorithms/param-ideal/hermite-matrices.jl b/src/algorithms/param-ideal/hermite-matrices.jl new file mode 100644 index 00000000..dd98eb77 --- /dev/null +++ b/src/algorithms/param-ideal/hermite-matrices.jl @@ -0,0 +1,175 @@ +@doc Markdown.doc""" + 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`. +""" +function hermite_matrix(gb::Vector{T}, b::Vector{T}, g::T)::MatElem where {T<:MPolyRingElem} + δ = length(b) + H = zero_matrix(coefficient_ring(gb[1]), δ, δ) + M_g = multiplication_matrix(gb, b, g) + M = [multiplication_matrix(gb, b, b[i]) for i in 1:δ] + for i in 1:δ + for j in i:δ + H[i, j] = tr(M[i] * M[j] * M_g) + end + end + for i in 2:δ + for j in 1:(i-1) + H[i, j] = H[j, i] + end + end + H +end + +@doc Markdown.doc""" + hermite_matrix(I::ParametricIdeal{K}, g::MPoly{K}; ) -> 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. +""" +function hermite_matrix( + I::ParametricIdeal{K}, + g::MPoly{K}; + retry::Int=10, + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::MatSpaceElem{K} where {K<:FracFieldElem} + input_hash = hash((I.gens, g)) + if haskey(I.hm, input_hash) + return I.hm[input_hash] + end + if show_progress + @info "Computing parametric Hermite matrix" + end + if I.prefer_gb + gb = groebner_basis(I; retry=retry, nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + b = _monomial_basis(gb) + g′ = _change_ring(g, parent(gb[1])) + H = hermite_matrix(gb, b, g′) + H = map(_canonicalize, H) + lock(I.hm_lock) do + I.hm[input_hash] = H + end + return H + end + δ = _expected_degree(I; worker_pool=worker_pool) + R₁, _ = polynomial_ring(QQ, I.symbols, internal_ordering=:degrevlex) + R₂ = base_ring(I.base_frac_field) + function bb(v) + hm_spec_hash = hash((I.gens, g, v)) + if haskey(I.hm_spec, hm_spec_hash) + return I.hm_spec[hm_spec_hash] + end + gb = groebner_basis(I, v; worker_pool=worker_pool) + b = _monomial_basis(gb) + if length(b) != δ + error("Error evaluating black box function: the ideal has unexpected degree at parameter values $(v).") + end + gₓ = _change_ring(map_coefficients(k -> evaluate(k, v), g), R₁) + Hₓ = hermite_matrix(gb, b, gₓ) + lock(I.hm_spec_lock) do + I.hm_spec[hm_spec_hash] = Hₓ + end + Hₓ + end + H = zero_matrix(fraction_field(R₂), δ, δ) + for i in 1:δ + for j in i:δ + H[i, j] = Interpolation.cuyt_lee( + R₂, v -> bb(v)[i, j]; + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric Hermite matrix, size=($δ,$δ), index=($i,$j)" + ) + end + end + for i in 2:δ + for j in 1:(i-1) + H[i, j] = H[j, i] + end + end + lock(I.hm_lock) do + I.hm[input_hash] = H + end + H +end + +# To compute Her(g^α), computing as above is a bit wasteful +# We proceed by first computing Her(1) and then multiplying by Mul(g)^α as needed +@doc Markdown.doc""" + hermite_matrix(g::Vector{MPoly{K}}, α::Vector{Int}, I::ParametricIdeal{K}; ) -> 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. +""" +function hermite_matrix( + g::Vector{MPoly{K}}, + α::Vector{Int}, + I::ParametricIdeal{K}; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::MatSpaceElem{K} where {K<:FracFieldElem} + R = parent(I.gens[1]) + H = hermite_matrix(I, one(R); nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + for i in 1:length(g) + if α[i] == 0 + continue + end + M = multiplication_matrix(I, g[i]; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + if α[i] == 1 + H = H * M + else + @assert α[i] == 2 + H = H * M * M + end + end + H = map(_canonicalize, H) + H +end + +@doc Markdown.doc""" + hermite_matrix(g::Vector{MPoly{K}}, α::Vector{Int}, I::ParametricIdeal{K}, vals::Vector{QQFieldElem}; ) -> 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. +""" +function hermite_matrix( + g::Vector{MPoly{K}}, + α::Vector{Int}, + I::ParametricIdeal{K}, + vals::Vector{QQFieldElem}; + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::QQMatrix where {K<:FracFieldElem} + R = parent(I.gens[1]) + H = map(a -> evaluate(numerator(a), vals) // evaluate(denominator(a), vals), hermite_matrix(I, one(R); nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)) + for i in 1:length(g) + if α[i] == 0 + continue + end + M = map(a -> evaluate(numerator(a), vals) // evaluate(denominator(a), vals), multiplication_matrix(I, g[i]; nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress)) + if α[i] == 1 + H = H * M + else + @assert α[i] == 2 + H = H * M * M + end + end + H +end diff --git a/src/algorithms/param-ideal/multiplication-matrices.jl b/src/algorithms/param-ideal/multiplication-matrices.jl new file mode 100644 index 00000000..6427fdf2 --- /dev/null +++ b/src/algorithms/param-ideal/multiplication-matrices.jl @@ -0,0 +1,104 @@ + +@doc Markdown.doc""" + multiplication_matrix(gb::Vector{T}, b::Vector{T}, g::T) -> MatElem + +Computes the multiplication matrix associated with the polynomial `g` with respect to the Groebner basis `gb` and the monomial basis `b`. +""" +function multiplication_matrix(gb::Vector{T}, b::Vector{T}, g::T)::MatElem where {T<:MPolyRingElem} + δ = length(b) + M = zero_matrix(coefficient_ring(gb[1]), δ, δ) + for i in 1:δ + prod = divrem(b[i] * g, gb)[2] + for j in 1:δ + M[j, i] = coeff(prod, b[j]) + end + end + M +end + +# For the multiplication matrices, we implement canonicalize functions +# to make computations more efficient + +function _canonicalize!(p::QQMPolyRingElem)::QQMPolyRingElem + combine_like_terms!(sort_terms!(p)) +end + +function _canonicalize(p::FracFieldElem{QQMPolyRingElem})::FracFieldElem{QQMPolyRingElem} + _canonicalize!(numerator(p)) // _canonicalize!(denominator(p)) +end + +@doc Markdown.doc""" + multiplication_matrix(I::ParametricIdeal{K}, g::MPoly{K}; ) -> MatSpaceElem{K} + +Computes the parametric multiplication matrix associated with the polynomial `g` with respect to the parametric ideal `I`. The multiplication matrix is computed by evaluating the ideal at a number of parameter values and interpolating the coefficients of the resulting multiplication matrices. + +# Arguments +- `retry::Int=10`: the maximum number of consecutive failures allowed when computing the multiplication 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 multiplication matrix. +""" +function multiplication_matrix( + I::ParametricIdeal{K}, + g::MPoly{K}; + retry::Int=10, + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=false +)::MatSpaceElem{K} where {K<:FracFieldElem} + input_hash = hash((I.gens, g)) + if haskey(I.mm, input_hash) + return I.mm[input_hash] + end + if show_progress + @info "Computing parametric multiplication matrix" + end + if I.prefer_gb + gb = groebner_basis(I; retry=retry, nr_thrds=nr_thrds, worker_pool=worker_pool, show_progress=show_progress) + b = _monomial_basis(gb) + g′ = _change_ring(g, parent(gb[1])) + M = multiplication_matrix(gb, b, g′) + M = map(_canonicalize, M) + lock(I.mm_lock) do + I.mm[input_hash] = M + end + return M + end + δ = _expected_degree(I; worker_pool=worker_pool) + R₁, _ = polynomial_ring(QQ, I.symbols, internal_ordering=:degrevlex) + R₂ = base_ring(I.base_frac_field) + function bb(v) + mm_spec_hash = hash((I.gens, g, v)) + if haskey(I.mm_spec, mm_spec_hash) + return I.mm_spec[mm_spec_hash] + end + gb = groebner_basis(I, v; worker_pool=worker_pool) + b = _monomial_basis(gb) + if length(b) != δ + error("Error evaluating black box function: the ideal has unexpected degree at parameter values $(v).") + end + gₓ = _change_ring(map_coefficients(k -> evaluate(k, v), g), R₁) + Mₓ = multiplication_matrix(gb, b, gₓ) + lock(I.mm_spec_lock) do + I.mm_spec[mm_spec_hash] = Mₓ + end + Mₓ + end + M = zero_matrix(fraction_field(R₂), δ, δ) + for i in 1:δ + for j in 1:δ + M[i, j] = Interpolation.cuyt_lee( + R₂, v -> bb(v)[i, j]; + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric multiplication matrix, size=($δ,$δ), index=($i,$j)" + ) + end + end + lock(I.mm_lock) do + I.mm[input_hash] = M + end + M +end diff --git a/src/algorithms/param-ideal/parametrizations.jl b/src/algorithms/param-ideal/parametrizations.jl new file mode 100644 index 00000000..f03a29be --- /dev/null +++ b/src/algorithms/param-ideal/parametrizations.jl @@ -0,0 +1,93 @@ +@doc Markdown.doc""" + rational_parametrization(I::ParametricIdeal{K}; ) -> (Vector{Symbol}, Vector{MPoly{K}}, MPoly{K}, MPoly{K}, Vector{MPoly{K}}) + +Computes the parametric rational parametrization associated with the parametric ideal `I`. The rational parametrization is computed by evaluating the ideal at a number of parameter values and interpolating the coefficients of the resulting rational parametrizations. + +# Arguments +- `retry::Int=10`: the maximum number of consecutive failures allowed when computing the rational parametrization 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=true`: whether to show a progress bar while computing the rational parametrization. +""" +function rational_parametrization( + I::ParametricIdeal{K} where K<:FracFieldElem; + retry::Int=10, + nr_thrds::Int=1, + worker_pool::AbstractWorkerPool=default_worker_pool(), + show_progress::Bool=true +) + n = I.num_vars + R = parent(I.gens[1]) + R₂ = base_ring(I.base_frac_field) + R₄, _ = polynomial_ring(fraction_field(R₂), symbols(R)[1:n], internal_ordering=:degrevlex) + param_shift = nothing + function bb(v) + gb = groebner_basis(I, v; worker_pool=worker_pool) + param = AlgebraicSolving.rational_parametrization(AlgebraicSolving.Ideal(gb)) + vars = param.vars + cfs_lf = param.cfs_lf + elim = param.elim / leading_coefficient(param.elim) + denom = param.denom / leading_coefficient(param.denom) + params = map(p -> p / leading_coefficient(param.denom), param.param) + param = (vars, cfs_lf, elim, denom, params) + if !isnothing(param_shift) + if param[1] != param_shift[1] || + param[2] != param_shift[2] || + degree(param[3]) != degree(param_shift[3]) || + degree(param[4]) != degree(param_shift[4]) || + degree.(param[5]) != degree.(param_shift[5]) + error("Error evaluating black box function: the rational parametrization has different monomial structure at different parameter values.") + end + end + param + end + param_shift = bb(QQ.(I.shift)) + param_vars = param_shift[1] + param_cfs_lf = param_shift[2] + if length(param_cfs_lf) > n + error("Error: the last variable in the rational parametrization is not a known variable: is the ideal radical?") + end + param_var_index = findfirst(i -> symbols(R₄)[i] == param_vars[end], 1:n) + param_var = gens(R₄)[param_var_index] + param_elim = zero(R₄) + for i in 0:degree(param_shift[3]) + c = Interpolation.cuyt_lee( + R₂, v -> coeff(bb(v)[3], i); + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric rational parametrization (eliminating polynomial), deg=($(degree(param_shift[3]))), index=($i)" + ) + param_elim += R₄(c) * param_var^i + end + param_denom = zero(R₄) + for i in 0:degree(param_shift[4]) + c = Interpolation.cuyt_lee( + R₂, v -> coeff(bb(v)[4], i); + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric rational parametrization (denominator), deg=($(degree(param_shift[4]))), index=($i)" + ) + param_denom += R₄(c) * param_var^i + end + param_params = [zero(R₄) for _ in 1:(n-1)] + len = n - 1 + degs = degree.(param_shift[5]) + for i in 1:len + for j in 0:degs[i] + c = Interpolation.cuyt_lee( + R₂, v -> coeff(bb(v)[5][i], j); + initial_shift=I.shift, + retry=retry, + nr_thrds=nr_thrds, + show_progress=show_progress, + desc="Parametric rational parametrization (parametrizations), size=($len,$(maximum(degs))), index=($i,$j)" + ) + param_params[i] += R₄(c) * param_var^j + end + end + (param_vars, param_cfs_lf, param_elim, param_denom, param_params) +end diff --git a/src/imports.jl b/src/imports.jl index 7a869c3c..07f20434 100644 --- a/src/imports.jl +++ b/src/imports.jl @@ -59,3 +59,8 @@ import Nemo: ZZMatrix, ZZRing, ZZRingElem + +import Nemo.Generic: + FracFieldElem, + MatSpaceElem, + MPoly diff --git a/src/interp/main.jl b/src/interp/main.jl index 6e09b9c6..c6035867 100644 --- a/src/interp/main.jl +++ b/src/interp/main.jl @@ -4,6 +4,7 @@ using Markdown using Nemo import Nemo.Generic: FracFieldElem using ..Progress +import ..AlgebraicSolving: _map_exponent_vectors, _change_ring export thiele, newton, cuyt_lee diff --git a/src/interp/resultant.jl b/src/interp/resultant.jl index 885012d4..26c8ba21 100644 --- a/src/interp/resultant.jl +++ b/src/interp/resultant.jl @@ -1,30 +1,3 @@ -function _map_exponent_vectors(f, p::MPolyRingElem, R::MPolyRing)::MPolyRingElem - cvzip = zip(coefficients(p), exponent_vectors(p)) - M = MPolyBuildCtx(R) - for (c, v) in cvzip - push_term!(M, coefficient_ring(R)(c), f(v)) - end - finish(M) -end - -function _change_ring(p::MPolyRingElem, R::MPolyRing)::MPolyRingElem - n = length(gens(parent(p))) - n′ = length(gens(R)) - if n == 0 || n′ == 0 - return R(constant_coefficient(p)) - end - indices = Vector{Int}() - for s in symbols(R) - i = findfirst(==(s), symbols(parent(p))) - if isnothing(i) - push!(indices, 0) - else - push!(indices, i) - end - end - return _map_exponent_vectors(v -> [i == 0 ? 0 : v[i] for i in indices], p, R) -end - # TODO: remove this once the next version of Nemo is released function _flint_discriminant(f::QQMPolyRingElem, i::Int)::QQMPolyRingElem R = parent(f) diff --git a/src/types.jl b/src/types.jl index 20896d6e..a0e9e054 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,4 +1,4 @@ -export Ideal +export Ideal, ParametricIdeal mutable struct RationalParametrization vars::Vector{Symbol} @@ -75,3 +75,75 @@ 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] + +mutable struct ParametricIdeal{K<:FracFieldElem} + base_frac_field::FracField + num_params::Int + num_vars::Int + symbols::Vector{Symbol} + gens::Vector{MPoly{K}} + sats::Vector{MPoly{K}} + dim::Union{Int,Nothing} + shift::Vector{Int} + gb::Dict{UInt64,Vector{MPoly{K}}} + gb_lock::ReentrantLock + gb_spec::Dict{UInt64,Vector{QQMPolyRingElem}} + gb_spec_lock::ReentrantLock + mm::Dict{UInt64,MatSpaceElem{K}} + mm_lock::ReentrantLock + mm_spec::Dict{UInt64,QQMatrix} + mm_spec_lock::ReentrantLock + hm::Dict{UInt64,MatSpaceElem{K}} + hm_lock::ReentrantLock + hm_spec::Dict{UInt64,QQMatrix} + hm_spec_lock::ReentrantLock + prefer_gb::Bool + radicalize::Bool + gens_alt::Vector{MPoly{K}} + + # Generate magic numbers for shifting the interpolation points to avoid singularities + # This works surprisingly well in practice, but it would be nice to have a more principled way to choose them + function _magic_shift(n::Int)::Vector{Int} + _magic_shift_cache = [11, 29, 47, 71, 97, 113, 149, 173, 197, 229] + for i in length(_magic_shift_cache)+1:n + push!(_magic_shift_cache, 2^i + rand(Int) % 2^i) + end + return [_magic_shift_cache[i] for i in 1:n] + end + + function ParametricIdeal(gens::Vector{MPoly{K}}; + sats::Vector{MPoly{K}}=Vector{MPoly{K}}(), + num_vars::Int=-1, + prefer_gb::Bool=true, + radicalize::Bool=false, + gens_alt::Vector{MPoly{K}}=Vector{MPoly{K}}() + ) where {K<:FracFieldElem} + @assert !isempty(gens) + I = new{K}() + R = parent(gens[1]) + I.base_frac_field = base_ring(R) + I.num_params = number_of_variables(base_ring(I.base_frac_field)) + I.num_vars = num_vars != -1 ? num_vars : number_of_variables(R) + I.symbols = symbols(R)[end-I.num_vars+1:end] + I.gens = gens + I.sats = sats + I.dim = nothing + I.shift = _magic_shift(I.num_params) + I.gb = Dict() + I.gb_lock = ReentrantLock() + I.gb_spec = Dict() + I.gb_spec_lock = ReentrantLock() + I.mm = Dict() + I.mm_lock = ReentrantLock() + I.mm_spec = Dict() + I.mm_spec_lock = ReentrantLock() + I.hm = Dict() + I.hm_lock = ReentrantLock() + I.hm_spec = Dict() + I.hm_spec_lock = ReentrantLock() + I.prefer_gb = prefer_gb + I.radicalize = radicalize + I.gens_alt = gens_alt + return I + end +end diff --git a/test/algorithms/param-ideal/groebner-bases.jl b/test/algorithms/param-ideal/groebner-bases.jl new file mode 100644 index 00000000..26aa2559 --- /dev/null +++ b/test/algorithms/param-ideal/groebner-bases.jl @@ -0,0 +1,17 @@ +import Nemo: fraction_field + +@testset "Algorithms -> Parametric Gröbner basis" begin + R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([a * x^2 + b * x + c]) + gb = groebner_basis(f) + @test gb == [(a * x^2 + b * x + c) / a] +end + +@testset "Algorithms -> Monomial basis" begin + R, (x,) = polynomial_ring(QQ, [:x]) + f = [x^2 - 2 * x + 1] + g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) + b = AlgebraicSolving._monomial_basis(g) + @test length(b) == 2 +end diff --git a/test/algorithms/param-ideal/hermite-matrices.jl b/test/algorithms/param-ideal/hermite-matrices.jl new file mode 100644 index 00000000..65987338 --- /dev/null +++ b/test/algorithms/param-ideal/hermite-matrices.jl @@ -0,0 +1,40 @@ +import Nemo: zero_matrix + +@testset "Algorithms -> Hermite matrix 1" begin + R, (x,) = polynomial_ring(QQ, [:x]) + f = [R(1)] + g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) + b = AlgebraicSolving._monomial_basis(g) + H₁ = AlgebraicSolving.hermite_matrix(g, b, one(R)) + @test H₁ == zero_matrix(QQ, 0, 0) +end + +@testset "Algorithms -> Hermite matrix 2" begin + R, (x,) = polynomial_ring(QQ, [:x]) + f = [x^2 - 2 * x + 1] + g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) + b = AlgebraicSolving._monomial_basis(g) + H₁ = AlgebraicSolving.hermite_matrix(g, b, one(R)) + @test H₁ == matrix(QQ, [2 2; 2 2]) + Hₓ = AlgebraicSolving.hermite_matrix(g, b, x) + Mₓ = AlgebraicSolving.multiplication_matrix(g, b, x) + @test Hₓ == H₁ * Mₓ +end + +@testset "Algorithms -> Parametric Hermite matrix" begin + R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([a * x^2 + b * x + c]) + g = one(R′) + H = AlgebraicSolving.hermite_matrix(f, g) + @test H == matrix(fraction_field(R), [2 -b//a; -b//a (b^2-2*a*c)//a^2]) +end + +@testset "Algorithms -> Parametric Hermite matrix with zero parameters" begin + R, () = polynomial_ring(QQ, Symbol[]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([x^2 - 2 * x + 1]) + g = one(R′) + H₁ = AlgebraicSolving.hermite_matrix(f, g) + @test H₁ == matrix(QQ, [2 2; 2 2]) +end diff --git a/test/algorithms/param-ideal/multiplication-matrices.jl b/test/algorithms/param-ideal/multiplication-matrices.jl new file mode 100644 index 00000000..1c2bb10f --- /dev/null +++ b/test/algorithms/param-ideal/multiplication-matrices.jl @@ -0,0 +1,25 @@ +import Nemo: identity_matrix, matrix + +@testset "Algorithms -> Multiplication matrix" begin + R, (x,) = polynomial_ring(QQ, [:x]) + f = [x^2 - 2 * x + 1] + g = AlgebraicSolving.groebner_basis(AlgebraicSolving.Ideal(f)) + b = AlgebraicSolving._monomial_basis(g) + M₁ = AlgebraicSolving.multiplication_matrix(g, b, one(R)) + @test M₁ == identity_matrix(QQ, 2) + Mₓ = AlgebraicSolving.multiplication_matrix(g, b, x) + @test Mₓ == matrix(QQ, [0 -1; 1 2]) +end + +@testset "Algorithms -> Parametric multiplication matrix" begin + M = AlgebraicSolving.multiplication_matrix + R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([a * x^2 + b * x + c]) + @test M(f, one(R′)) == matrix(fraction_field(R), + [1 0; 0 1]) + @test M(f, x) == matrix(fraction_field(R), + [0 -c//a; 1 -b//a]) + @test M(f, x^2) == matrix(fraction_field(R), + [-c//a (b*c)//a^2; -b//a (b^2-a*c)//a^2]) +end diff --git a/test/algorithms/param-ideal/parametrizations.jl b/test/algorithms/param-ideal/parametrizations.jl new file mode 100644 index 00000000..2ff04397 --- /dev/null +++ b/test/algorithms/param-ideal/parametrizations.jl @@ -0,0 +1,14 @@ +import Nemo.Generic: FracFieldElem, MPoly + +@testset "Algorithms -> Parametric rational parametrization" begin + R, (a, b, c) = polynomial_ring(QQ, [:a, :b, :c]) + R′, (x,) = polynomial_ring(fraction_field(R), [:x], internal_ordering=:degrevlex) + f = ParametricIdeal([a * x^2 + b * x + c]) + @test rational_parametrization(f) == ( + [:x], + ZZRingElem[], + x^2 + b // a * x + c // a, + x + (1 // 2 * b) // a, + MPoly{FracFieldElem{QQMPolyRingElem}}[] + ) +end diff --git a/test/runtests.jl b/test/runtests.jl index 052157b7..878a3aad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,10 @@ include("algorithms/hilbert.jl") include("algorithms/decomposition.jl") include("algorithms/param-curves.jl") include("algorithms/sampling.jl") +include("algorithms/param-ideal/groebner-bases.jl") +include("algorithms/param-ideal/multiplication-matrices.jl") +include("algorithms/param-ideal/hermite-matrices.jl") +include("algorithms/param-ideal/parametrizations.jl") include("examples/katsura.jl") include("interp/thiele.jl") include("interp/newton.jl")