Merge pull request #1192 from n0rbed/RootFinding

RootFinding GSOC project

.JuliaFormatter.toml

@@ -0,0 +1,2 @@
+style = "sciml"
+format_markdown = true

Project.toml

@@ -26,6 +26,7 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
@@ -43,13 +44,15 @@ TermInterface = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 
 [extensions]
 SymbolicsForwardDiffExt = "ForwardDiff"
 SymbolicsGroebnerExt = "Groebner"
 SymbolicsLuxCoreExt = "LuxCore"
+SymbolicsNemoExt = "Nemo"
 SymbolicsPreallocationToolsExt = ["PreallocationTools", "ForwardDiff"]
 SymbolicsSymPyExt = "SymPy"
 
@@ -75,6 +78,7 @@ LogExpFunctions = "0.3"
 LuxCore = "0.1.11"
 MacroTools = "0.5"
 NaNMath = "1"
+Nemo = "0.45, 0.46"
 PreallocationTools = "0.4"
 PrecompileTools = "1"
 RecipesBase = "1.1"
@@ -107,4 +111,4 @@ SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "Lux", "ComponentArrays"]
+test = ["Test", "SafeTestsets", "Pkg", "PkgBenchmark", "PreallocationTools", "ForwardDiff", "Groebner", "BenchmarkTools", "ReferenceTests", "SymPy", "Random", "Lux", "ComponentArrays", "Nemo"]

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -11,7 +13,9 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 [compat]
 BenchmarkTools = "1.3"
 Documenter = "1"
+Groebner = "0.7"
 Latexify = "0.15, 0.16"
+Nemo = "0.45, 0.46"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StaticArrays = "1.5"

docs/make.jl

@@ -42,6 +42,7 @@ makedocs(
             "manual/expression_manipulation.md",
             "manual/derivatives.md",
             "manual/groebner.md",
+            "manual/solver.md",
             "manual/arrays.md",
             "manual/build_function.md",
             "manual/functions.md",

docs/src/manual/solver.md

@@ -0,0 +1,69 @@
+# Solver
+
+The main symbolic solver for Symbolics.jl is `symbolic_solve`. Symbolic solving
+means that it only uses symbolic (algebraic) methods and outputs exact solutions.
+
+```@docs
+Symbolics.symbolic_solve
+```
+
+One other symbolic solver is `symbolic_linear_solve` which is limited compared to 
+`symbolic_solve` as it only solves linear equations.
+```@docs
+Symbolics.symbolic_linear_solve
+```
+
+`symbolic_solve` only supports symbolic, i.e. non-floating point computations, and thus prefers equations
+where the coefficients are integer, rational, or symbolic. Floating point coefficients are transformed into
+rational values and BigInt values are used internally with a potential performance loss, and thus it is recommended
+that this functionality is only used with floating point values if necessary. In contrast, `symbolic_linear_solve`
+directly handles floating point values using standard factorizations.
+
+### More technical details and examples
+
+#### Technical details
+
+The `symbolic_solve` function uses 4 hidden solvers in order to solve the user's input. Its base,
+`solve_univar`, uses analytic solutions up to polynomials of degree 4 and factoring as its method
+for solving univariate polynomials. The function's `solve_multipoly` uses GCD on the input polynomials then throws passes the result
+to `solve_univar`. The function's `solve_multivar` uses Groebner basis and a separating form in order to create linear equations in the
+input variables and a single high degree equation in the separating variable [^1]. Each equation resulting from the basis is then passed
+to `solve_univar`. We can see that essentially, `solve_univar` is the building block of `symbolic_solve`. If the input is not a valid polynomial and can not be solved by the algorithm above, `symbolic_solve` passes
+it to `ia_solve`, which attempts solving by attraction and isolation [^2]. This only works when the input is a single expression
+and the user wants the answer in terms of a single variable. Say `log(x) - a == 0` gives us `[e^a]`.
+
+#### Nice examples
+
+```@example solver
+using Symbolics, Nemo;
+@variables x;
+Symbolics.symbolic_solve(9^x + 3^x ~ 8, x)
+```
+
+```@example solver
+@variables x y z;
+Symbolics.symbolic_linear_solve(2//1*x + y - 2//1*z ~ 9//1*x, 1//1*x)
+```
+
+```@example solver
+using Groebner;
+@variables x y z;
+
+eqs = [x^2 + y + z - 1, x + y^2 + z - 1, x + y + z^2 - 1]
+Symbolics.symbolic_solve(eqs, [x,y,z])
+```
+
+### Feature completeness
+
+- [x] Linear and polynomial equations
+- [x] Systems of linear and polynomial equations
+- [x] Some transcendental functions
+- [x] Systems of linear equations with parameters (via `symbolic_linear_solve`)
+- [ ] Equations with radicals
+- [ ] Systems of polynomial equations with parameters and positive dimensional systems
+- [ ] Inequalities
+
+# References
+
+[^1]: [Rouillier, F. Solving Zero-Dimensional Systems Through the Rational Univariate Representation. AAECC 9, 433–461 (1999).](https://doi.org/10.1007/s002000050114)
+[^2]: [R. W. Hamming, Coding and Information Theory, ScienceDirect, 1980](https://www.sciencedirect.com/science/article/pii/S0747717189800070).

ext/SymbolicsGroebnerExt.jl

@@ -1,6 +1,7 @@
 module SymbolicsGroebnerExt
 
 using Groebner
+const Nemo = Groebner.Nemo
 
 if isdefined(Base, :get_extension)
     using Symbolics
@@ -70,4 +71,241 @@ function Symbolics.is_groebner_basis(polynomials::Vector{Num}; kwargs...)
     Groebner.isgroebner(polynoms; kwargs...)
 end
 
+### Solver ###
+
+# Map each variable of the given poly.
+# Can be used to transform Nemo polynomial to expression.
+function nemo_crude_evaluate(poly::Nemo.MPolyRingElem, varmap)
+    new_poly = 0
+    for (i, term) in enumerate(Nemo.terms(poly))
+        new_term = nemo_crude_evaluate(Nemo.coeff(poly, i), varmap)
+        for var in Nemo.vars(term)
+            exp = Nemo.degree(term, var)
+            exp == 0 && continue
+            new_var = varmap[var]
+            new_term *= new_var^exp
+        end
+        new_poly += new_term
+    end
+    new_poly
+end
+
+function nemo_crude_evaluate(poly::Nemo.FracElem, varmap)
+    nemo_crude_evaluate(numerator(poly), varmap) // nemo_crude_evaluate(denominator(poly), varmap)
+end
+
+function nemo_crude_evaluate(poly::Nemo.ZZRingElem, varmap)
+    BigInt(poly)
+end
+
+function gen_separating_var(vars)
+    n = 1
+    new_var = (Symbolics.@variables _T)[1]
+    present = any(isequal(new_var, var) for var in vars)
+    while present
+        new_var = Symbolics.variables(repeat("_", n) * "_T")[1]
+        present = any(isequal(new_var, var) for var in vars)
+        n += 1
+    end
+    return new_var
+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)
+
+    all_vars = [vars..., params...]
+    nemo_ring, nemo_all_vars = Nemo.polynomial_ring(Nemo.QQ, map(string, all_vars))
+
+    sym_to_nemo = Dict(all_vars .=> nemo_all_vars)
+    nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
+    nemo_gb = Symbolics.substitute(gb, sym_to_nemo)
+    nemo_gb = Symbolics.substitute(nemo_gb, sym_to_nemo)
+
+    nemo_vars = filter(v -> string(v) in string.(vars), nemo_all_vars)
+    nemo_params = filter(v -> string(v) in string.(params), nemo_all_vars)
+
+    ring_flat = parent(nemo_vars[1])
+    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)
+    result = empty(gb_demoted)
+    while true
+        gb_demoted = map(f -> Nemo.map_coefficients(c -> c // Nemo.leading_coefficient(f), f), gb_demoted)
+        for i in 1:length(gb_demoted)
+            f = gb_demoted[i]
+            f_nf = Nemo.normal_form(f, result)
+            if !iszero(f_nf)
+                push!(result, f_nf)
+            end
+        end
+        isequal(gb_demoted, result) && break
+        gb_demoted = result
+        result = empty(result)
+    end
+    @assert all(f -> isone(Nemo.leading_coefficient(f)), result)
+
+    sym_to_nemo = Dict(sym => nem for sym in all_vars for nem in [vars_demoted..., params_demoted...] if isequal(string(sym),string(nem)))
+    nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
+
+    final_result = Num[]
+
+    for i in eachindex(result)
+
+        monoms = collect(Nemo.monomials(result[i]))
+        coeffs = collect(Nemo.coefficients(result[i]))
+
+        poly = 0
+        for j in eachindex(monoms)
+            poly += nemo_crude_evaluate(coeffs[j], nemo_to_sym) * nemo_crude_evaluate(monoms[j], nemo_to_sym)
+        end
+        push!(final_result, poly)
+    end
+
+    final_result
+end
+
+function solve_zerodim(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, warns=true)
+    # Reference: Rouillier, F. Solving Zero-Dimensional Systems
+    # Through the Rational Univariate Representation.
+    # AAECC 9, 433–461 (1999). https://doi.org/10.1007/s002000050114
+
+    rng = Groebner.Random.Xoshiro(42)
+
+    all_indeterminates = reduce(union, map(Symbolics.get_variables, eqs))
+    params = map(Symbolics.Num ∘ Symbolics.wrap, setdiff(all_indeterminates, vars))
+
+    # Use a new variable to separate the input polynomials (Reference above)
+    new_var = gen_separating_var(vars)
+    old_len = length(vars)
+    vars = vcat(vars, new_var)
+
+    new_eqs = []
+    generating = true
+    n_iterations = 1
+    separating_form = new_var
+
+    while generating
+        new_eqs = copy(eqs)
+        separating_form = new_var
+        for i = 1:(old_len)
+            separating_form += BigInt(rand(rng, -n_iterations:n_iterations))*vars[i]
+        end
+
+        if isequal(separating_form, new_var)
+            continue
+        end
+
+        push!(new_eqs, separating_form)
+
+        new_eqs = Symbolics.groebner_basis(new_eqs, ordering=Lex(vcat(vars, params)))
+
+        # handle "unsolvable" case
+        if isequal(1, new_eqs[1])
+            return []
+        end
+
+        new_eqs = demote(new_eqs, vars, params)
+        new_eqs = map(Symbolics.unwrap, new_eqs)
+
+        # condition for positive dimensionality, i.e. infinite solutions
+        if length(new_eqs) < length(vars)
+            warns && @warn("Infinite number of solutions")
+            return nothing
+        end
+
+        # Exit in the Shape Lemma case:
+        # g(T, params) = 0
+        # x1 - f1(T, params) = 0
+        # ...
+        # xn - fn(T, params) = 0
+        generating = !(length(new_eqs) == length(vars))
+        if length(new_eqs) == length(vars)
+            generating |= !(isequal(setdiff(Symbolics.get_variables(new_eqs[1]), params), [new_var]))
+            for i in eachindex(new_eqs)[2:end]
+                present_vars = setdiff(Symbolics.get_variables(new_eqs[i]), new_var)
+                present_vars = setdiff(present_vars, params)
+                isempty(present_vars) && (generating = false; break;)
+                var_i = present_vars[1]
+                condition1 = isequal(present_vars, [var_i])
+                condition2 = Symbolics.degree(new_eqs[i], var_i) == 1
+                generating |= !(condition1 && condition2)
+            end
+        end
+
+        # non-cyclic case
+        n_iterations > 10 && return []
+
+        n_iterations += 1
+    end
+
+    solutions = []
+
+    # first, solve the first minimal polynomial
+    @assert length(new_eqs) == length(vars)
+    @assert isequal(setdiff(Symbolics.get_variables(new_eqs[1]), params), [new_var])
+    minpoly_sols = Symbolics.symbolic_solve(Symbolics.wrap(new_eqs[1]), new_var, dropmultiplicity=dropmultiplicity)
+    solutions = [Dict{Num, Any}(new_var => sol) for sol in minpoly_sols]
+
+    new_eqs = new_eqs[2:end]
+
+    # second, iterate over eqs and sub each found solution
+    # then add the roots of the remaining unknown variables 
+    for (i, eq) in enumerate(new_eqs)
+        present_vars = setdiff(Symbolics.get_variables(eq), params)
+        present_vars = setdiff(present_vars, new_var)
+        @assert length(present_vars) == 1
+        var_tosolve = present_vars[1]
+        @assert Symbolics.degree(eq, var_tosolve) == 1
+        @assert !isempty(solutions)
+        for roots in solutions
+            subbded_eq = Symbolics.substitute(eq, Dict([new_var => roots[new_var]]); fold=false)
+            subbded_eq = Symbolics.substitute(subbded_eq, Dict([var_tosolve => 0]); fold=false)
+            new_var_sols = [-subbded_eq]
+            @assert length(new_var_sols) == 1
+            root = new_var_sols[1]
+            roots[var_tosolve] = root
+        end
+    end
+
+    vars = vars[1:end-1]
+    for roots in solutions
+        delete!(roots, new_var)
+    end
+
+    return solutions
+end
+
+function transendence_basis(sys, vars)
+    J = Symbolics.jacobian(sys, vars)
+    x0 = Dict(v => rand(-10:10) for v in vars)
+    J_x0 = substitute(J, x0)
+    rk, rref = Nemo.rref(Nemo.matrix(Nemo.QQ, J_x0))
+    pivots = Int[]
+    for i in 1:length(sys)
+        col = findfirst(!iszero, rref[i, :])
+        !isnothing(col) && push!(pivots, col)
+    end
+    vars[setdiff(collect(1:length(vars)), pivots)]
+end
+
+function Symbolics.solve_multivar(eqs::Vector, vars::Vector{Num}; dropmultiplicity=true, warns=true)
+    sol = solve_zerodim(eqs, vars; dropmultiplicity=dropmultiplicity, warns=warns)
+    !isnothing(sol) && return sol
+    tr_basis = transendence_basis(eqs, vars)
+    @info "" tr_basis
+    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
+
 end # module