From 6a7bb743dd28b64beb2792871eae64233690214a Mon Sep 17 00:00:00 2001 From: n0rbed Date: Sun, 25 Aug 2024 07:36:03 +0300 Subject: [PATCH 01/11] changed output format of zero dimension sols --- ext/SymbolicsGroebnerExt.jl | 5 ----- src/solver/main.jl | 5 ++++- test/solver.jl | 4 ++-- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index 0d9c3f8fb..0baf87103 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -295,11 +295,6 @@ function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplici isempty(tr_basis) && return nothing vars_gen = setdiff(vars, tr_basis) sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns) - for roots in sol - for x in tr_basis - roots[x] = x - end - end sol end diff --git a/src/solver/main.jl b/src/solver/main.jl index be9fd73a7..14b664b9a 100644 --- a/src/solver/main.jl +++ b/src/solver/main.jl @@ -187,7 +187,10 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where isequal(sols, nothing) && return nothing for sol in sols for var in x - sol[var] = postprocess_root(sol[var]) + root = get(sol, var, missing) + if !isequal(wrap(root), missing) + sol[var] = postprocess_root(sol[var]) + end end end diff --git a/test/solver.jl b/test/solver.jl index 6c9e4a905..8a91f0ba0 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -182,7 +182,7 @@ end @test isequal(symbolic_solve([x^4 - 1, x - 2], [x]), []) # TODO: test this properly - sol = symbolic_solve([x^3 + 1, x*y^3 - 1], [x, y]) +# sol = symbolic_solve([x^3 + 1, x*y^3 - 1], [x, y]) eqs = [x*y + 2x^2, y^2 -1] arr_calcd_roots = sort_arr(symbolic_solve(eqs, [x,y]), [x,y]) @@ -278,7 +278,7 @@ end @variables t w u v sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w,u,v]) - @test isequal(sol, [Dict(u => u, t => -5 / (-1 + u + v), v => v, w => 1 - u - v)]) + @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)]) end @testset "Factorisation" begin From 048d7f0a4822a158723a1e77965ae1851cd98e7b Mon Sep 17 00:00:00 2001 From: n0rbed Date: Sun, 25 Aug 2024 07:46:59 +0300 Subject: [PATCH 02/11] simpler fix --- src/solver/main.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solver/main.jl b/src/solver/main.jl index 14b664b9a..2e4b91c2a 100644 --- a/src/solver/main.jl +++ b/src/solver/main.jl @@ -187,8 +187,7 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where isequal(sols, nothing) && return nothing for sol in sols for var in x - root = get(sol, var, missing) - if !isequal(wrap(root), missing) + if haskey(sol, var) sol[var] = postprocess_root(sol[var]) end end From a2705410f83c881dbc05469449502236f025b6d9 Mon Sep 17 00:00:00 2001 From: n0rbed Date: Sun, 25 Aug 2024 07:49:52 +0300 Subject: [PATCH 03/11] uncommented test --- test/solver.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/solver.jl b/test/solver.jl index 8a91f0ba0..c931b8567 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -182,7 +182,7 @@ end @test isequal(symbolic_solve([x^4 - 1, x - 2], [x]), []) # TODO: test this properly -# sol = symbolic_solve([x^3 + 1, x*y^3 - 1], [x, y]) + sol = symbolic_solve([x^3 + 1, x*y^3 - 1], [x, y]) eqs = [x*y + 2x^2, y^2 -1] arr_calcd_roots = sort_arr(symbolic_solve(eqs, [x,y]), [x,y]) From 8aac5980dab1150a558f0a5bd66b7256cd905fad Mon Sep 17 00:00:00 2001 From: n0rbed Date: Mon, 26 Aug 2024 03:35:14 +0300 Subject: [PATCH 04/11] warn for non-cyclic case --- ext/SymbolicsGroebnerExt.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index 0baf87103..8c0abf07d 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -233,7 +233,10 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa end # non-cyclic case - n_iterations > 10 && return [] + if n_iterations > 10 + warns && @warn("symbolic_solve can not currently solve this system of polynomials.") + return nothing + end n_iterations += 1 end From fd3f498e9a5d42471b25b8f520eb89e87276c884 Mon Sep 17 00:00:00 2001 From: n0rbed Date: Tue, 27 Aug 2024 12:35:37 +0300 Subject: [PATCH 05/11] added test --- src/solver/main.jl | 8 ++++---- test/solver.jl | 8 +++++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/solver/main.jl b/src/solver/main.jl index 2e4b91c2a..414972eca 100644 --- a/src/solver/main.jl +++ b/src/solver/main.jl @@ -154,8 +154,8 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where sols = [] if expr_univar sols = check_poly_inunivar(expr, x) ? - solve_univar(expr, x, dropmultiplicity = dropmultiplicity) : - ia_solve(expr, x, warns = warns) + solve_univar(expr, x, dropmultiplicity=dropmultiplicity, warns=warns) : + ia_solve(expr, x, warns=warns) isequal(sols, nothing) && return nothing else for i in eachindex(expr) @@ -165,7 +165,7 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where end end sols = solve_multipoly( - expr, x, dropmultiplicity = dropmultiplicity, warns = warns) + expr, x, dropmultiplicity=dropmultiplicity, warns=warns) isequal(sols, nothing) && return nothing end @@ -183,7 +183,7 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where end end - sols = solve_multivar(expr, x, dropmultiplicity = dropmultiplicity) + sols = solve_multivar(expr, x, dropmultiplicity=dropmultiplicity, warns=warns) isequal(sols, nothing) && return nothing for sol in sols for var in x diff --git a/test/solver.jl b/test/solver.jl index c931b8567..58ac51588 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -236,13 +236,13 @@ end # cyclic 3 @variables z1 z2 z3 eqs = [z1 + z2 + z3, z1*z2 + z1*z3 + z2*z3, z1*z2*z3 - 1] - sol = Symbolics.symbolic_solve(eqs, [z1,z2,z3]) + sol = symbolic_solve(eqs, [z1,z2,z3]) backward = [Symbolics.substitute(eqs, s) for s in sol] @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward) @variables x y eqs = [2332//232*x + 2131232*y - 1//343434, x + y + 1] - sol = Symbolics.symbolic_solve(eqs, [x,y]) + sol = symbolic_solve(eqs, [x,y]) backward = [Symbolics.substitute(eqs, s) for s in sol] @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward) @@ -259,10 +259,12 @@ end # at most 4 roots by Bézout's theorem rand_eq(xs, d) = rand(-10:10) + rand(-10:10)*x + rand(-10:10)*y + rand(-10:10)*x*y + rand(-10:10)*x^2 + rand(-10:10)*y^2 eqs = [rand_eq([x,y],2), rand_eq([x,y],2)] - sol = Symbolics.symbolic_solve(eqs, [x,y]) + sol = symbolic_solve(eqs, [x,y]) backward = [Symbolics.substitute(eqs, s) for s in sol] @test all(x -> all(isapprox.(eval(Symbolics.toexpr(x)), 0; atol=1e-6)), backward) end + + @test isnothing(symbolic_solve([x^2, x*y, y^2], [x,y], warns=false)) end @testset "Multivar parametric" begin From 58a045a5d503d95bfc532e9bf06725b80cbbd626 Mon Sep 17 00:00:00 2001 From: n0rbed Date: Tue, 27 Aug 2024 12:58:54 +0300 Subject: [PATCH 06/11] woops --- src/solver/main.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solver/main.jl b/src/solver/main.jl index 414972eca..e4152e0ec 100644 --- a/src/solver/main.jl +++ b/src/solver/main.jl @@ -154,7 +154,7 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where sols = [] if expr_univar sols = check_poly_inunivar(expr, x) ? - solve_univar(expr, x, dropmultiplicity=dropmultiplicity, warns=warns) : + solve_univar(expr, x, dropmultiplicity=dropmultiplicity) : ia_solve(expr, x, warns=warns) isequal(sols, nothing) && return nothing else @@ -238,7 +238,7 @@ implemented in the function `get_roots` and its children. # Examples """ -function solve_univar(expression, x; dropmultiplicity = true) +function solve_univar(expression, x; dropmultiplicity=true) args = [] mult_n = 1 expression = unwrap(expression) From 3cc5fd0bcd182d05adc595590d477a3166b628b9 Mon Sep 17 00:00:00 2001 From: n0rbed Date: Wed, 28 Aug 2024 02:30:20 +0300 Subject: [PATCH 07/11] some changes --- ext/SymbolicsGroebnerExt.jl | 11 +++++++++-- src/solver/main.jl | 5 +---- test/solver.jl | 2 +- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index 8c0abf07d..438ddf460 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -298,6 +298,13 @@ function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplici isempty(tr_basis) && return nothing vars_gen = setdiff(vars, tr_basis) sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns) + + for roots in sol + for x in tr_basis + roots[x] = x + end + end + sol end @@ -311,8 +318,8 @@ PrecompileTools.@setup_workload begin PrecompileTools.@compile_workload begin symbolic_solve(equation1, x) symbolic_solve(equation_actually_polynomial) - symbolic_solve(simple_linear_equations, [x, y]) - symbolic_solve(equations_intersect_sphere_line, [x, y, z]) + symbolic_solve(simple_linear_equations, [x, y], warns=false) + symbolic_solve(equations_intersect_sphere_line, [x, y, z], warns=false) end end diff --git a/src/solver/main.jl b/src/solver/main.jl index e4152e0ec..237fee425 100644 --- a/src/solver/main.jl +++ b/src/solver/main.jl @@ -121,10 +121,6 @@ function symbolic_solve(expr, x::T; dropmultiplicity = true, warns = true) where for var in x check_x(var) end - if length(x) == 1 - x = x[1] - x_univar = true - end end if !(expr isa Vector) @@ -213,6 +209,7 @@ function symbolic_solve(expr; x...) vars = wrap.(vars) @assert all(v isa Num for v in vars) "All variables should be Nums or BasicSymbolics" + vars = isone(length(vars)) ? vars[1] : vars return symbolic_solve(expr, vars; x...) end diff --git a/test/solver.jl b/test/solver.jl index 58ac51588..3dff04a47 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -279,7 +279,7 @@ end @test isnothing(symbolic_solve([x*y - a, sin(x)], [x, y])) @variables t w u v - sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w,u,v]) + sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w]) @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)]) end From c16933bbb2f10a20da89d539c91671842a48387b Mon Sep 17 00:00:00 2001 From: n0rbed Date: Fri, 30 Aug 2024 13:08:11 +0300 Subject: [PATCH 08/11] symbolic_solve([...], [x]) added --- ext/SymbolicsGroebnerExt.jl | 11 ++++++++++- test/solver.jl | 3 +++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index 438ddf460..ebf5174a5 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -108,6 +108,7 @@ end # Given a GB in k[params][vars] produces a GB in k(params)[vars] function demote(gb, vars::Vector{Num}, params::Vector{Num}) isequal(gb, [1]) && return gb + gb = Symbolics.wrap.(SymbolicUtils.toterm.(gb)) Symbolics.check_polynomial.(gb) @@ -126,7 +127,7 @@ function demote(gb, vars::Vector{Num}, params::Vector{Num}) ring_param, params_demoted = Nemo.polynomial_ring(Nemo.base_ring(ring_flat), map(string, nemo_params)) ring_demoted, vars_demoted = Nemo.polynomial_ring(Nemo.fraction_field(ring_param), map(string, nemo_vars), internal_ordering=:lex) varmap = Dict((nemo_vars .=> vars_demoted)..., (nemo_params .=> params_demoted)...) - gb_demoted = map(f -> nemo_crude_evaluate(f, varmap), nemo_gb) + gb_demoted = map(f -> ring_demoted(nemo_crude_evaluate(f, varmap)), nemo_gb) result = empty(gb_demoted) while true gb_demoted = map(f -> Nemo.map_coefficients(c -> c // Nemo.leading_coefficient(f), f), gb_demoted) @@ -176,6 +177,7 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa # Use a new variable to separate the input polynomials (Reference above) new_var = gen_separating_var(vars) old_len = length(vars) + old_vars = deepcopy(vars) vars = vcat(vars, new_var) new_eqs = [] @@ -204,6 +206,13 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa return [] end + for i in reverse(eachindex(new_eqs)) + all_present = Symbolics.get_variables(new_eqs[i]) + if length(intersect(all_present, vars)) < 1 + deleteat!(new_eqs, i) + end + end + new_eqs = demote(new_eqs, vars, params) new_eqs = map(Symbolics.unwrap, new_eqs) diff --git a/test/solver.jl b/test/solver.jl index 3dff04a47..def1507bf 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -281,6 +281,9 @@ end @variables t w u v sol = symbolic_solve([t*w - 1 ~ 4, u + v + w ~ 1], [t,w]) @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)]) + + sol = symbolic_solve([x-y, y-z], [x]) + @test isequal(sol, [Dict(x=>z)]) end @testset "Factorisation" begin From b0e130cee5587b82fca38bd406bc83d94748f950 Mon Sep 17 00:00:00 2001 From: n0rbed Date: Fri, 30 Aug 2024 16:50:09 +0300 Subject: [PATCH 09/11] changed filter location --- ext/SymbolicsGroebnerExt.jl | 16 +++++++++------- test/solver.jl | 5 ++++- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index ebf5174a5..3cd142dda 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -206,13 +206,6 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa return [] end - for i in reverse(eachindex(new_eqs)) - all_present = Symbolics.get_variables(new_eqs[i]) - if length(intersect(all_present, vars)) < 1 - deleteat!(new_eqs, i) - end - end - new_eqs = demote(new_eqs, vars, params) new_eqs = map(Symbolics.unwrap, new_eqs) @@ -301,10 +294,19 @@ function transendence_basis(sys, vars) end function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, warns=true) + for i in reverse(eachindex(eqs)) + present_vars = Symbolics.get_variables(eqs[i]) + if length(intersect(present_vars, vars)) < 1 && length(present_vars) != 0 + deleteat!(eqs, i) + end + end + sol = solve_zerodim(eqs, vars; dropmultiplicity=dropmultiplicity, warns=warns) !isnothing(sol) && return sol + tr_basis = transendence_basis(eqs, vars) isempty(tr_basis) && return nothing + vars_gen = setdiff(vars, tr_basis) sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns) diff --git a/test/solver.jl b/test/solver.jl index def1507bf..44b9a8901 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -283,7 +283,10 @@ end @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)]) sol = symbolic_solve([x-y, y-z], [x]) - @test isequal(sol, [Dict(x=>z)]) + @test isequal(sol, [Dict(x=>y)]) + + sol = symbolic_solve([x-y, y-z], [x, y]) + @test isequal(sol, [Dict(x=>z, y=>z)]) end @testset "Factorisation" begin From 60de1a148de739f0259f973dcf2c4d4ed84fb55c Mon Sep 17 00:00:00 2001 From: n0rbed Date: Fri, 30 Aug 2024 19:35:35 +0300 Subject: [PATCH 10/11] reverted to original filter position --- ext/SymbolicsGroebnerExt.jl | 16 +++++++--------- test/solver.jl | 5 ++++- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/ext/SymbolicsGroebnerExt.jl b/ext/SymbolicsGroebnerExt.jl index 3cd142dda..ebf5174a5 100644 --- a/ext/SymbolicsGroebnerExt.jl +++ b/ext/SymbolicsGroebnerExt.jl @@ -206,6 +206,13 @@ function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, wa return [] end + for i in reverse(eachindex(new_eqs)) + all_present = Symbolics.get_variables(new_eqs[i]) + if length(intersect(all_present, vars)) < 1 + deleteat!(new_eqs, i) + end + end + new_eqs = demote(new_eqs, vars, params) new_eqs = map(Symbolics.unwrap, new_eqs) @@ -294,19 +301,10 @@ function transendence_basis(sys, vars) end function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, warns=true) - for i in reverse(eachindex(eqs)) - present_vars = Symbolics.get_variables(eqs[i]) - if length(intersect(present_vars, vars)) < 1 && length(present_vars) != 0 - deleteat!(eqs, i) - end - end - sol = solve_zerodim(eqs, vars; dropmultiplicity=dropmultiplicity, warns=warns) !isnothing(sol) && return sol - tr_basis = transendence_basis(eqs, vars) isempty(tr_basis) && return nothing - vars_gen = setdiff(vars, tr_basis) sol = solve_zerodim(eqs, vars_gen; dropmultiplicity=dropmultiplicity, warns=warns) diff --git a/test/solver.jl b/test/solver.jl index 44b9a8901..5cc678fc2 100644 --- a/test/solver.jl +++ b/test/solver.jl @@ -283,10 +283,13 @@ end @test isequal(sol, [Dict(t => -5 / (-1 + u + v), w => 1 - u - v)]) sol = symbolic_solve([x-y, y-z], [x]) - @test isequal(sol, [Dict(x=>y)]) + @test isequal(sol, [Dict(x=>z)]) sol = symbolic_solve([x-y, y-z], [x, y]) @test isequal(sol, [Dict(x=>z, y=>z)]) + + sol = symbolic_solve([x + y - z, y - z], [x]) + @test isequal(sol, [Dict(x=>0)]) end @testset "Factorisation" begin From 1ebb2c94d12024ccd775e88936633905a827c1de Mon Sep 17 00:00:00 2001 From: n0rbed Date: Fri, 30 Aug 2024 21:54:49 +0300 Subject: [PATCH 11/11] small --- docs/src/manual/solver.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/src/manual/solver.md b/docs/src/manual/solver.md index 7b4a33013..3380d6632 100644 --- a/docs/src/manual/solver.md +++ b/docs/src/manual/solver.md @@ -63,6 +63,14 @@ Symbolics.symbolic_solve(eqs, [x,y,z]) - [ ] Systems of polynomial equations with parameters and positive dimensional systems - [ ] Inequalities +### Expressions we can not solve (but aim to) +``` +# Mathematica + +In[1]:= Reduce[x^2 - x - 6 > 0, x] +Out[1]= x < -2 || x > 3 +``` + # References [^1]: [Rouillier, F. Solving Zero-Dimensional Systems Through the Rational Univariate Representation. AAECC 9, 433–461 (1999).](https://doi.org/10.1007/s002000050114)