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,68 @@
+# 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`)
+- [ ] Systems of polynomial equations with parameters
+- [ ] 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

@@ -70,4 +70,115 @@ function Symbolics.is_groebner_basis(polynomials::Vector{Num}; kwargs...)
     Groebner.isgroebner(polynoms; kwargs...)
 end
 
+function Symbolics.solve_multivar(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
+
+    # Use a new variable to separate the input polynomials (Reference above)
+    new_var = Symbolics.gen_separating_var(vars)
+    old_len = length(vars)
+    push!(vars, new_var)
+
+    new_eqs = []
+    generating = true
+    n_iterations = 1
+
+    while generating
+        new_eqs = copy(eqs)
+        eq = new_var
+        for i = 1:(old_len)
+            eq += BigInt(rand(-n_iterations:n_iterations))*vars[i]
+        end
+
+        if isequal(eq, new_var)
+            continue
+        end
+
+        push!(new_eqs, eq)
+
+        new_eqs = Symbolics.groebner_basis(new_eqs, ordering=Lex(vars))
+
+        if length(new_eqs) <= length(vars) 
+            generating &= false
+        end
+
+        for i  in eachindex(new_eqs)[2:end]
+            generating |= all(Symbolics.degree(var) > 1 for var in Symbolics.get_variables(new_eqs[i]))
+        end
+
+        if isequal(new_eqs[1], new_var)
+            generating = true
+        end
+
+        n_iterations += 1
+    end
+
+    solutions = []
+
+    # handle "unsolvable" cases
+    if isequal(1, new_eqs[1])
+        return solutions
+    end
+    if length(new_eqs) < length(vars)
+        warns && (@warn("Infinite number of solutions"); return nothing) || return nothing
+    end
+
+    new_eqs = SymbolicUtils.toterm.(new_eqs)
+
+    # first, solve any single variable equations
+    i = 1
+    while !(i > length(new_eqs))
+            present_vars = Symbolics.get_variables(new_eqs[i])
+        for var in vars
+            if size(present_vars, 1) == 1 && isequal(var, present_vars[1])
+                new_sols = Symbolics.solve_univar(Symbolics.wrap(new_eqs[i]), var, dropmultiplicity=dropmultiplicity)
+
+                if length(solutions) == 0
+                    append!(solutions, [Dict{Num, Any}(var => sol) for sol in new_sols])
+                else
+                    solutions = Symbolics.add_sol_to_all(solutions, new_sols, var)
+                end
+
+                deleteat!(new_eqs, i)
+                i = i - 1
+                break
+            end
+        end
+        i = i + 1
+    end
+
+
+    # second, iterate over eqs and sub each found solution
+    # then add the roots of the remaining unknown variables 
+    for eq in new_eqs
+        solved = false
+        present_vars = Symbolics.get_variables(eq)
+        size_of_sub = length(solutions[1])
+
+        if size(present_vars, 1) <= (size_of_sub + 1)
+            while !solved 
+                subbed_eq = eq
+                for (var, root) in solutions[1]
+                    subbed_eq = Symbolics.substitute(subbed_eq, Dict([var => root]), fold=false)
+                end
+
+                var_tosolve = Symbolics.get_variables(subbed_eq)[1]
+                new_var_sols = Symbolics.solve_univar(subbed_eq, var_tosolve, dropmultiplicity=dropmultiplicity)
+                Symbolics.add_sol!(solutions, new_var_sols, var_tosolve, 1)
+
+                solved = all(x -> length(x) == size_of_sub+1, solutions)
+            end
+        end
+    end
+
+    pop!(vars)
+    for roots in solutions
+        delete!(roots, new_var)
+    end
+
+    return solutions
+end
+
 end # module

ext/SymbolicsNemoExt.jl

@@ -0,0 +1,127 @@
+module SymbolicsNemoExt
+using Nemo
+
+if isdefined(Base, :get_extension)
+    using Symbolics
+    using Symbolics: Num, symtype
+else
+    using ..Symbolics
+    using ..Symbolics: Num, symtype
+end
+
+# 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)
+    Rational(poly)
+end
+
+# factor(x^2*y + b*x*y - a*x - a*b)  ->  (x*y - a)*(x + b)
+function Symbolics.factor_use_nemo(poly::Num)
+    Symbolics.check_polynomial(poly)
+    Symbolics.degree(poly) == 0 && return poly, Num[]
+    vars = Symbolics.get_variables(poly)
+    nemo_ring, nemo_vars = Nemo.polynomial_ring(Nemo.QQ, map(string, vars))
+    sym_to_nemo = Dict(vars .=> nemo_vars)
+    nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
+    nemo_poly = Symbolics.substitute(poly, sym_to_nemo)
+    nemo_fac = Nemo.factor(nemo_poly)
+    nemo_unit = Nemo.unit(nemo_fac)
+    nemo_factors = collect(keys(nemo_fac.fac)) 
+    sym_unit = Rational(Nemo.coeff(nemo_unit, 1))
+    sym_factors = map(f -> Symbolics.wrap(nemo_crude_evaluate(f, nemo_to_sym)), nemo_factors)
+
+    for (i, fac) in enumerate(sym_factors)
+        sym_factors[i] = fac^(collect(values(nemo_fac.fac))[i])
+    end
+
+    return sym_unit, sym_factors
+end
+
+# gcd(x^2 - y^2, x^3 - y^3) -> x - y
+function Symbolics.gcd_use_nemo(poly1::Num, poly2::Num)
+    Symbolics.check_polynomial(poly1)
+    Symbolics.check_polynomial(poly2)
+    vars1 = Symbolics.get_variables(poly1)
+    vars2 = Symbolics.get_variables(poly2)
+    vars = vcat(vars1, vars2)
+    nemo_ring, nemo_vars = Nemo.polynomial_ring(Nemo.QQ, map(string, vars))
+    sym_to_nemo = Dict(vars .=> nemo_vars)
+    nemo_to_sym = Dict(v => k for (k, v) in sym_to_nemo)
+    nemo_poly1 = Symbolics.substitute(poly1, sym_to_nemo)
+    nemo_poly2 = Symbolics.substitute(poly2, sym_to_nemo)
+    nemo_gcd = Nemo.gcd(nemo_poly1, nemo_poly2)
+    sym_gcd = Symbolics.wrap(nemo_crude_evaluate(nemo_gcd, nemo_to_sym))
+    return sym_gcd
+end
+
+
+function Symbolics.demote(gb, vars::Vector{Num}, params::Vector{Num})
+    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 = [v for (k, v) in sym_to_nemo if any(isequal(k, var) for var in vars)]
+    nemo_params = [v for (k, v) in sym_to_nemo if any(isequal(k, param) for param in params)]
+
+    ring_flat = parent(nemo_vars[1])
+    ring_param, params_demoted = Nemo.polynomial_ring(base_ring(ring_flat), map(string, nemo_params))
+    ring_demoted, vars_demoted = Nemo.polynomial_ring(fraction_field(ring_param), map(string, nemo_vars), internal_ordering=Nemo.internal_ordering(ring_flat))
+    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)
+    for i in 1:length(gb_demoted)
+        gb_demoted = map(f -> map_coefficients(c -> c // leading_coefficient(f), f), gb_demoted)
+        f = gb_demoted[i]
+        f_nf = Nemo.normal_form(f, result)
+        if !iszero(f_nf)
+            push!(result, f_nf)
+        end
+    end
+
+    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 = []
+
+    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
+
+end # module