Return exact roots from the rule-based backend

Project.toml

@@ -32,7 +32,7 @@ Nemo = "0.51, 0.52, 0.53, 0.54, 0.55, 0.56"
 PolyLog = "2.6.0"
 PrecompileTools = "1.2.1"
 SpecialFunctions = "2"
-SymbolicUtils = "4.27"
+SymbolicUtils = "4.49"
 Symbolics = "7"
 TermInterface = "2"
 Logging = "1.10"

docs/src/methods/risch.md

@@ -158,6 +158,48 @@ integrate(cos(x), x)            # Transformed to exponential form
 integrate(tan(x), x)            # Uses differential field extension
 ```
 
+## Exact algebraic coefficients
+
+A rational function whose denominator does not factor over the rationals needs
+an algebraic extension for its partial fraction decomposition. The roots
+involved are computed exactly, in the field of algebraic numbers, and
+coefficients of degree at most two are returned as exact radicals rather than
+as floating point approximations:
+
+```julia
+integrate(1/(x^2 - 2), x, RischMethod())
+# (1//4)*log(sqrt(2) - x)*sqrt(2) - (1//4)*log(sqrt(2) + x)*sqrt(2)
+
+integrate(1/(x^3 - 1), x, RischMethod())
+# (1//3)*log(1 - x) - (1//6)*log(1 + x + x^2) +
+#     (1//3)*atan(-(1//3)*sqrt(3) - (2//3)*sqrt(3)*x)*sqrt(3)
+
+integrate(1/(x^4 + 1), x, RischMethod())
+# (1//4)*atan(-1 + sqrt(2)*x)*sqrt(2) + (1//4)*atan(1 + sqrt(2)*x)*sqrt(2) -
+#     (1//8)*log(1 - sqrt(2)*x + x^2)*sqrt(2) +
+#     (1//8)*log(1 + sqrt(2)*x + x^2)*sqrt(2)
+```
+
+Radicals are reported in lowest terms, so the largest perfect square factor is
+pulled out of the radicand: a coefficient of `sqrt(12//49)` is returned as
+`(2//7)*sqrt(3)`.
+
+Roots of degree three or higher have no such closed form and are kept as a
+`Root` placeholder holding the exact algebraic number.
+
+!!! note
+    Simplifying an antiderivative keeps it exact. This requires
+    SymbolicUtils ≥ 4.49, whose `simplify` leaves a radical such as `sqrt(2)`
+    exact, which is why the package bounds it there.
+
+    `simplify` still does not always cancel the difference between an
+    antiderivative's derivative and the integrand down to zero — for
+    `1/(x^2 - 2)` and `1/(x^2 - 3)` it stops at an uncancelled product of
+    conjugate factors, which is exact but not visibly zero. To check such a
+    result, evaluate the residual numerically at a few points, or compare
+    against an expected expression with `isequal`, rather than relying on the
+    difference simplifying to zero.
+
 ## Limitations
 
 The Risch method, following Bronstein's book, does not handle:

src/methods/risch/frontend.jl

@@ -98,16 +98,18 @@ function to_symb(t::QQBarFieldElem)
     if degree(f)==2 && iszero(coeff(f,1))
         y = to_symb(-coeff(f,0)//coeff(f, 2))
         if y>=0
+            term_sqrt = SymbolicUtils.simplify(SymbolicUtils.term(sqrt, y))
             if t==maximum(conjugates(t))
-                return sqrt(y)
+                return term_sqrt
             else
-                return -sqrt(y)
+                return -term_sqrt
             end
         else
+            term_sqrt = SymbolicUtils.simplify(SymbolicUtils.term(sqrt, -y))
             if imag(t)==maximum(imag.(conjugates(t)))
-                return sqrt(-y)*1im
+                return term_sqrt*im
             else
-                return -sqrt(-y)*1im
+                return -term_sqrt*im
             end
         end
     elseif degree(f)==2 # coeff(f,1)!=0

src/methods/rule_based/rule2.jl

@@ -253,7 +253,17 @@ function check_expr_r(data::SymsType, rule::Expr, matches::MatchDict)::MatchDict
     end
 
     ((Symbol(operation(data)) !== rule.args[1]) && !neim_pass) && return FAIL_DICT::MatchDict # :):):)
-    (length(arg_data) != length(arg_rule)) && return FAIL_DICT::MatchDict # :):):)
+    if length(arg_data) != length(arg_rule)
+        # ((6.5)) A numeric coefficient folds into a single factor, an exact
+        # constant does not: `-1.4142*x` has two arguments while `-sqrt(2)*x`
+        # has three, [-1, x, sqrt(2)]. A pattern such as `(~!b)*(~x)`, whose
+        # `b` is meant to bind the whole coefficient, then cannot match on
+        # length alone. Retry once with the constant factors collapsed into a
+        # single one. Only reached after the equal-length match has failed, so
+        # this can create a match but never break one.
+        (rule.args[1]===:*) && return collapsed_constants_match(arg_data, arg_rule, matches)
+        return FAIL_DICT::MatchDict # :):):)
+    end
 
     # ((7))
     if (rule.args[1]===:+) || (rule.args[1]===:*)
@@ -271,6 +281,32 @@ function check_expr_r(data::SymsType, rule::Expr, matches::MatchDict)::MatchDict
     return ceoaa(arg_data, arg_rule, matches)::MatchDict
 end
 
+"""
+    collapsed_constants_match(arg_data, arg_rule, matches)
+
+Retry matching a product after collapsing its constant factors into one.
+
+`arg_data` has more factors than `arg_rule` when a coefficient is an exact
+constant that does not fold into a single term, such as the `[-1, x, sqrt(2)]`
+of `-sqrt(2)*x`. Multiplying the constant factors together restores the shape
+the pattern expects, letting a slot bind the whole coefficient.
+"""
+function collapsed_constants_match(arg_data, arg_rule, matches::MatchDict)::MatchDict
+    consts = SymsType[]
+    others = SymsType[]
+    for a in arg_data
+        push!(constant_value(a) === nothing ? others : consts, a)
+    end
+    (length(consts) < 2 || length(others) + 1 != length(arg_rule)) &&
+        return FAIL_DICT::MatchDict
+    collapsed = vcat(SymsType[reduce(*, consts)], others)
+    for perm in permutations(collapsed)
+        m = ceoaa(perm, arg_rule, matches)
+        m !== FAIL_DICT && return m::MatchDict
+    end
+    return FAIL_DICT::MatchDict
+end
+
 """
 check expression of all arguments.
 This function iterates over arguments of data and rule (both at the same time)

src/methods/rule_based/rules_loader.jl

@@ -29,7 +29,9 @@ function load_rules(rules_paths)
         # add rules
         Base.include_string(@__MODULE__, read(file, String), file)
         local_file_rules = Base.invokelatest(() -> file_rules) # Use Base.invokelatest to handle world age issues in Julia 1.12+
-        append!(RULES, [x[2] for x in local_file_rules])
+        # Rewrite literal square roots on the replacement side only; the match
+        # side must keep the `sqrt` the matcher looks for by name.
+        append!(RULES, [exact_roots_in_rule(x[2]) for x in local_file_rules])
         append!(IDENTIFIERS, [x[1] for x in local_file_rules])
     end
     print("\e[1A\e[2K\e[1A\e[2K")
@@ -90,24 +92,27 @@ function reload_rules(path; verbose = true)
 
     for r in local_file_rules
         idx = findfirst(i->identifier_ge(i, r[1]), IDENTIFIERS)
+        # Same rewrite as in `load_rules`: a replacement side installed raw would
+        # build its square roots with Julia's `sqrt` and return floats.
+        rule = exact_roots_in_rule(r[2])
 
         # if there is a identifier >= of r[1]
         if idx !== nothing
             # if r[1] is already in the identifiers
             if IDENTIFIERS[idx]==r[1]
                 # replace rule
-                RULES[idx] = r[2]
+                RULES[idx] = rule
                 verbose && printstyled("replaced rule $(r[1]) at index $idx\n";color = :yellow)
             # else add it at idx
             else
                 insert!(IDENTIFIERS, idx, r[1])
-                insert!(RULES, idx, r[2])
+                insert!(RULES, idx, rule)
                 verbose && printstyled("Inserted rule $(r[1]) at index $idx\n";color=:green)
             end
         # else add it at the end
         else
             push!(IDENTIFIERS, r[1])
-            push!(RULES, r[2])
+            push!(RULES, rule)
             verbose && printstyled("Appended rule $(r[1]) at the end of RULES (index $(length(RULES)))\n";color = :magenta)
         end
     end