Symbolica 2.2: symbolic integration

Technical articles and release notes about Symbolica, symbolic computation, numerical methods, and high-performance computer algebra.
symbolica
release
rust
python
integration
Author

Ben Ruijl

Published

July 24, 2026

Introduction

Symbolica 2.2, the latest release of our high-performance symbolic computation framework for Python and Rust, adds fast, step-by-step symbolic integration. Under the hood is a native Rust port of more than 7,000 integration rules, tested against a 72,944-problem corpus. The integrator can return both an antiderivative and the rule-by-rule path used to find it.

Here is a simple integration example:

from symbolica import *
x = E('x')
(1/(x**2+1)).integrate(x) # ➜ atan(x)
use symbolica::prelude::*;
use symbolica_integrate::Integrate;

fn main() {
    let x = symbol!("x");
    parse!("1/(x^2+1)").integrate(x).unwrap(); // ➜ atan(x)
}

Let us first explore why symbolic integration is difficult, and how we can implement it in Symbolica, taking inspiration from symbolic differentiation.

Symbolic differentiation

Differentiation is a common operation and is straightforward to implement. Finding an antiderivative — the inverse operation — is much more difficult. In this section we will explore why this is the case.

A derivative operation on an expression can be cast into a pattern matching rule, by virtue of the chain rule. The chain rule says that:

\[ \frac{d}{dx} f(g(x)) = f'(g(x)) \cdot g'(x) \]

which separates the derivative of the argument of the function from the derivative of the function itself. For example, the derivative of sin(f(x)) is cos(f(x))*f'(x). This means that we can implement differentiation recursively by pattern matching on the expression and applying the appropriate rules.

Below the derivative operation is implemented for a simple expression tree in Rust:

#[derive(Clone)]
enum Expression {
    Num(isize),
    Var(usize),
    Add(Box<(Expression, Expression)>),
    Mul(Box<(Expression, Expression)>),
    Sin(Box<Expression>),
    Cos(Box<Expression>),
}

impl Expression {
    fn derivative(&self, var: usize) -> Expression {
        use Expression::*;
        match self {
            Num(_) => Num(0),
            Var(v) => if *v == var { Num(1) } else { Num(0) },
            Add(arg) => Add(Box::new((arg.0.derivative(var), arg.1.derivative(var)))),
            Mul(arg) => Add(
                Box::new((Mul(Box::new((arg.0.derivative(var), arg.1.clone()))),
                          Mul(Box::new((arg.0.clone(), arg.1.derivative(var)))),
                ),
            )),
            Sin(arg) => {
                Mul(Box::new((Cos(arg.clone()), arg.derivative(var))))
            }
            Cos(arg) => Mul(Box::new((
                Num(-1),
                Mul(Box::new((
                    Sin(arg.clone()),
                    arg.derivative(var),
                ))),
            ))),
        }
    }
}

Running these rules backward is not straightforward. When differentiating, the composition \(f(g(x))\) is visible, so the chain rule tells us exactly what to do. When integrating, we see only the resulting expression and must infer both \(f\) and \(g\), if such a decomposition exists at all. For example,

\[ \frac{d}{dx}e^{-x^2}=-2x e^{-x^2}. \]

The factor \(-2x\) reveals the derivative of the exponent, making \(-2x e^{-x^2}\) easy to integrate. Without that factor, however, \(e^{-x^2}\) has no elementary antiderivative and requires introducing the error function.

Symbolic integration

As you well know from calculus, the integral of \(f'(x)\) with respect to \(x\) is \(f(x) + C\), where \(C\) is the constant of integration. As a result there are infinitely many antiderivatives and some are shorter or more useful than others.

For example

\[ \int \frac{1}{e^x + 1} dx \] has two completely different-looking antiderivatives \(x - \log(e^x + 1)\) and \(-2 \text{arctanh}(1 + 2e^x)\). These results differ by a constant, but that is far from obvious because they use different functions. What answer should the integrator return? Ideally, it should return a compact expression written in functions familiar to the user.

Over the years, various methods for symbolic integration have been designed. For example, Trager’s algorithm can be used to integrate any rational function, and it has been part of Symbolica for a while. The downside is that the output of the algorithm is written in terms of (complex) logarithms, which is not always the most useful form for the end-user. For example, an antiderivative of \(\frac{1}{x^2 + 1}\) is \(\arctan(x)\), but Trager’s algorithm will return \(\frac{i}{2}\log(1 - i x) - \frac{i}{2}\log(1 + i x)\).

Trigonometric functions and special functions look nice for the user, but they are dangerous for a computer algebra system, as they offer more ways to write 0 in a complicated way (e.g., \(\cos(x)^2 + \sin(x)^2 -1\)), which can lead to accidental divisions by 0. In fact, detecting if an expression is zero is undecidable in general!

There is the famous Risch decision process, which can be used to integrate elementary functions1 whose antiderivative can be expressed in terms of elementary functions. However, it has notoriously complicated corner cases and no complete implementation is known. The Risch algorithm does not support special functions, which means it cannot (out of the box) handle error functions or polylogarithms. See the Math for Machines lecture slides for a construction of the Risch algorithm.

1 Elementary functions are built from constants and variables using finitely many arithmetic operations, exponentials, logarithms, algebraic functions, and composition.

Pattern matching for integration

We return to the idea of pattern matching, which worked so well for symbolic differentiation. We will try to implement a set of rules that can be applied to an expression to find an antiderivative. This is not a complete algorithm, but it works for many expressions that are encountered in practice.

Instead of using Rust’s structural pattern matcher (match) like we did for the derivative, we will use Symbolica’s pattern matcher, which is more suited to match mathematical expressions. It supports commutative and associative matching, and has wildcards (variables ending in underscores) that can match to any subexpression.

Let’s use the following rules:

  • a constant in \(x\) integrates to that constant times \(x\).
  • a product of a constant in \(x\) and an expression integrates to the constant times the integral of the expression.
  • a sum integrates to the sum of the integrals.
  • a power \(x^n\) integrates to \(\frac{1}{n+1} x^{n+1}\), except when \(n=-1\), which integrates to \(\log(x)\).

Written in Python, with the Symbolica pattern matching library, this looks like the following2:

2 The for-loop is a terse way to get the next match and the outcome of the next match. In Rust we would have written if let Some(m) = e.matches(...).next().

from symbolica import *

integral, a_, a__, b__, c_, x_, x, y, z = S(
    'integral', 'a_', 'a__', 'b__', 'c_', 'x_', 'x', 'y', 'z')

# Integrate wrt fixed variable x
def integrate(e: Expression) -> Expression:
    for m in e.match(a__, ~a__.req_contains(x), partial=False):
        return m[a__] * x

    for m in e.match(a__ * b__, ~a__.req_contains(x) & b__.req_contains(x), partial=False):
        return m[a__] * integrate(m[b__])

    for m in e.match(a_ + b__, partial=False):
        return integrate(m[a_]) + integrate(m[b__])

    for m in e.match(x**a_.opt(), a_ != -1, partial=False):
        return x**(m[a_]+1) / (m[a_]+1)

    for m in e.match(x**-1, partial=False):
        return x.log()

    return integral(e)
e = integrate(1/x**2 + y/x + x**2 + 3*x + 2)
print(e) # 2*x + y*log(x) + 3/2*x^2 + 1/3*x^3 - 1/x

These few rules nicely integrate polynomials in \(x\) and even include negative powers of \(x\). We can make the patterns broader. For any powers of an affine function in \(x\), we have the following formula:

\[ \begin{align*} \int (a+ b x)^{-1} dx &= \frac{1}{b} \log(a + b x)\\ \int (a+ b x)^n dx &= \frac{(a + b x)^{n+1}}{b (n+1)}; \quad n \neq -1 \end{align*} \]

This relation holds for any \(a\) and \(b \neq0\) that are independent of \(x\), for example \((3+2x)^3\), \((y x)^2\) and \(x\) itself. We can implement this as a pattern, using optional wildcards a__ and b__ (see Section 6). The optionality means that the pattern will also match simpler expressions that are not structurally the same (when \(a=0\) there would be no sum).

Note that the patterns for powers of \(x\) are now redundant, as they are a special case of the more general affine patterns, so we can remove them!

Below is the updated integration function:

# Integrate wrt fixed variable x
def integrate(e: Expression) -> Expression:
    for m in e.match(a__, ~a__.req_contains(x), partial=False):
        return m[a__] * x

    for m in e.match(a__ * b__, ~a__.req_contains(x) & b__.req_contains(x), partial=False):
        return m[a__] * integrate(m[b__])

    for m in e.match(a_ + b__, partial=False):
        return integrate(m[a_]) + integrate(m[b__])

    for m in e.match((a_.opt() + b__.opt()*x)**-1, ~a_.req_contains(x) & ~b__.req_contains(x), partial=False):
        return (m[a_] + m[b__]*x).log() / m[b__]

    for m in e.match((a_.opt() + b__.opt()*x)**c_.opt(), ~a_.req_contains(x) & ~b__.req_contains(x) & ~c_.req_contains(x), partial=False):
        return (m[a_] + m[b__]*x)**(m[c_]+1) / (m[b__]*(m[c_]+1))

    return integral(e)
e = integrate((2+x)**2 + 1/(y+x) + (2+y*x)**-3)
print(e) # log(x+y)-1/2/(y*(2+x*y)^2)+1/3*(2+x)^3

This new function can now integrate a broader class of expressions, including powers of affine functions in \(x\). The question now becomes, how far can we go? Can we implement a set of rules that can integrate any expression?

Rubi

In 2008, Albert Rich started work on Rubi (Rule-based Integrator), a system for symbolic integration based on pattern matching. Over the years it has grown to a set of over 7,000 rules that can integrate a very broad class of expressions. The rules themselves are written in Mathematica, but Rubi is open-source and can be used in other languages.

Here is an example of a more complicated Rubi rule:

Int[x_^m_.*(a_.+b_.*ArcCot[c_.*x_^n_])^p_,x_Symbol] :=
  With[{k=Denominator[m]},
  k \[Star] Subst[Int[x^(k*(m+1)-1)*(a+b*ArcCot[c*x^(k*n)])^p,x],x,x^(1/k)]] /;
FreeQ[{a,b,c},x] && IGtQ[p,1] && IGtQ[n,0] && FractionQ[m]

In a nutshell, it rewrites the integral of x^m*(a+b*acot(c*x^n))^p where m is a fraction into a case where m is an integer by substituting x with x^(1/k), where k is the denominator of m.

The Rubi rules come with a corpus of 72,944 integrals and it can integrate 99.9% of them. It is the most complete symbolic integration system available (outperforming Mathematica’s own integration suite) and it also manages to produce more compact output than most other systems.

Naturally, having the Rubi system natively in Symbolica has been on the back of our minds for quite a while now.

Porting 7,000+ rules

Porting the Rubi rules to a new language is a daunting task due to its sheer volume and intricate control flow. However, there is also a large mechanical component to it, which made it possible to port the rules to Symbolica with the help of present-day AI assistance. We used Codex 5.5 High, and later Codex 5.6 Sol High. The port took about two weeks.

The first 48 hours of (wasted) work were essentially unsupervised. We gave instructions to port Mathematica .nb files into Symbolica and gave it Symbolica documentation. This led to extremely verbose code that did not use Symbolica’s pattern matching capabilities at all. The AI also generated rules that were not part of Rubi. We also noticed that it would preserve whatever (poorly designed) output format it invented for the first few rules rather than refactor it later.

Degradation of the output was also observed in \goal mode, where simple instructions to follow Rubi rules as closely as possible were ignored. We figured this was because of the required context window being too large. One obvious way to reduce context pressure was to preprocess the Mathematica notebooks into compact rule files and prevent verbose code output.

It was time to try again.

Making the port work

The second attempt was much more successful. We changed the Rubi input rules to simple auto-generated .m text files (rather than verbose .nb files), gave concise instructions on how the pattern matcher worked, and gave a very strict output format, with some manually ported rules as examples. The output format is macro-heavy, so that it is as short as possible, and as close as possible to the original Rubi rules. Here is an example of the port of the complicated acot rule above:

fn push_rules_rule_5368(rules: &mut Vec<RubiRule>) {
    rubi_symb!(a__, b__, c__, m_, n_, p_, x_);
    rules.push(rubi_rule!(
        order: 5368,
        source: "Int[x_^m_.*(a_.+b_.*ArcCot[c_.*x_^n_])^p_,x_Symbol] :=
          With[{k=Denominator[m]},
          k \\[Star] Subst[Int[x^(k*(m+1)-1)*(a+b*ArcCot[c*x^(k*n)])^p,x],x,x^(1/k)]] /;
        FreeQ[{a,b,c},x] && IGtQ[p,1] && IGtQ[n,0] && FractionQ[m]",
        pattern: x_.pow(m_) * (a__ + b__ * (c__ * x_.pow(n_)).acot()).pow(p_),
        with: [m_, a__, b__, c__, n_, p_, x_],
        optional: [m_, a__, b__, c__],
        when: {
            freeq!([a__, b__, c__], x_)
                && igtq!(p_, 1)
                && igtq!(n_, 0)
                && fractionq!(m_)
        },
        rhs: {
            let k = Atom::num(rubi_denominator(&m_));
            let substitution_guard = fresh_x().unwrap();
            let x = substitution_guard.symbol();
            let payload = x.pow(&k * (&m_ + 1) - 1)
                * (&a__ + &b__ * (&c__ * x.pow(&k * &n_)).acot()).pow(&p_);
            let primitive = rubi_rhs_int(&payload, x);
            rubi_star(k, rubi_subst(&primitive, x, x_.pow(1 / &k)))
        },
    ));
}

Note that the rule contains the original Mathematica source code, so that for future debugging, the AI has the ground truth in its context. It also contains an order entry. This is an important detail, as the way Mathematica rules are defined in the Rubi source files is not the order in which they are applied. The actual execution order has to be queried in a live Mathematica instance.

During testing and debugging, the AI was still keen on inventing new manipulations on the right-hand-side (rhs) of the substitution, so we had to double down in the instructions to follow the original Rubi rules as closely as possible. What ultimately helped the most was that we wrote a testing setup that can trace the path of the integration process in Rust and in Mathematica, and gave clear instructions that every integration step should follow the original path of Rubi running in Mathematica.

The AI was also keen on working around (perceived and real) bugs in Symbolica or unexpected normalization behavior. We gave instructions to not pave over any bugs but write separate MRE files for us to inspect and fix. Sometimes, the AI came to a wrong conclusion about why a rule wasn’t matching, and would discover while writing an MRE that the problem was somewhere else.

72,944 integrals later

In the end, we have translated all rules to Symbolica’s pattern matching language, and they are now available in Symbolica 2.2. The rules are implemented in a separate MIT-licensed crate symbolica-integrate. For the Python build, this crate is baked into the main Symbolica package.

To the best of our knowledge, Symbolica provides the only complete port of the latest Rubi 4.17 integration system outside the Wolfram Language, preserving its 7,000+ ordered rules and passing the complete 72,944-problem Rubi corpus. SymJa is the only other complete port based on Rubi 4.16, while other public ports are partial, experimental or discontinued.

Benchmarks

symbolica-integrate processes the complete 72,944-problem corpus in 18 minutes of wall time, on a Ryzen 9 5900X parallelized over 8 cores. The distribution of per-problem times had a 57 ms median, 118 ms average, and 10.8 s maximum.

These timings are competitive, as an integration run on the independent test suite consisting of 1,892 problems shows:

Integrator Time
symbolica-integrate 111.24 s
Rubi 4.17.3.0 in Mathematica 13.2 155.78 s

In this run, symbolica-integrate was 1.40× faster.

The benchmarks can be reproduced by running the rubi_corpus test in the repository. All timings were run without derivative verification, which is a separate step that can be run after the integration is complete.

Showing every integration step

Since Rubi is based on pattern matching, it is explainable: one can trace which rule was applied to obtain a particular integral. Using this, you can generate a step-by-step explanation of the integration process, which is useful for educational purposes or for debugging complex integrals.

For example:

from symbolica import *
x = E('x')
result, overview, steps = (x/(x+1)).integrate_with_steps(x)
print(overview)
use symbolica::prelude::*;
use symbolica_integrate::Integrate; // requires `steps` feature

fn main() {
    let x = symbol!("x");
    let explanation = parse!("x/(x+1)").integrate_with_steps(x);
    println!("{}\n", explanation);
}

gives

∫ x/(1+x) dx = ∫ 1-1/(1+x) dx
  ∫ 1-1/(1+x) dx = ∫ 1 dx+∫ 1/(-1-x) dx
    ∫ 1 dx = x
    ∫ -1/(1+x) dx = -log(1+x)
= x-log(1+x)

Each returned IntegrationStep object contains the original Rubi rule, and an explanation of the rule. My hope is that the community will contribute more explanations to the rules, so that the integration process can be explained in more detail.

Pattern matching upgrades

To support the Rubi rules, Symbolica’s pattern matching system has been upgraded. The following features have been added: optional wildcards and alternatives.

Optional wildcards

Symbolica matches structurally, which means that x^n_ will not match the variable x, as there is no literal power in this expression. Sometimes this is undesirable, for example in the earlier case where we wanted to integrate a power of x using

\(\int x^n dx = \frac{1}{n+1}x^{n+1}; \quad n \neq -1\)

which we also wanted to apply to the case \(n=1\). For this purpose, we can wrap the wildcard n_ in the function opt to produce the optional pattern x^opt(n_). This pattern will first try to match x^n_, and when it fails, it will try to match a substructure with n set to a particular value dependent on the context.

Optional pattern Matches Value assigned to n_
x^opt(n_) x 1
x*opt(n_) x 1
x+opt(n_) x 0
f(opt(n_)) f Empty

Optional wildcards can be used in combination with each other to match vastly different structures. For example:

a__, b__, c__, p_, x = S('a__', 'b__', 'c__', 'p_', 'x')

pat = (a__.opt() + b__.opt() + x)**p_.opt() * c__.opt()

((1 + 3 * x)**3 + 4).replace(pat, 1) # ➜ 1
((1 + 3 * x)**3).replace(pat, 1) # ➜ 1
(1 + 2*x).replace(pat, 1) # ➜ 1
x.replace(pat, 1) # ➜ 1

Alternative patterns

Patterns with alternatives can be created using the built-in alt function. For example alt(x,y) will match either x or y:

x, y, z = S('x', 'y', 'z')
pat = E('alt(x, y)') * z

(x*z).replace(pat, 1) # ➜ 1
(y*z).replace(pat, 1) # ➜ 1

Each alternative must have the same wildcards, so that the bindings are the same in every branch. Therefore, alt(x_,y_) is invalid.

Alternative patterns were not needed for the Rubi port, but they came as a byproduct of the work on optional wildcards.

Further updates

For more updates and a migration guide for Symbolica 2.2, see the release notes.

Here are some more changes:

  • Option to factor over the complex numbers
  • Algebraic dependency extraction when converting to polynomials (improves functions such as together with exp and roots)
  • Add specification of a minimum probability density to continuous grids
  • Term iterator for Python API
  • Improved n-th root rendering for LaTeX and Typst
  • Performance improvements for pattern matching
TipSupport Symbolica

Symbolica is free for hobbyists and a single-core instance is free for non-commercial use. If you want new features, more documentation, or support for your workflow, your organization can purchase a license and directly support the project.

Closing

Symbolica 2.2 brings fast, traceable symbolic integration to Rust and Python. Its 7,000+ rules are available in the MIT-licensed symbolica-integrate crate, and the complete Rubi corpus provides a foundation for keeping the implementation correct as it evolves. Next, we want to investigate whether AI can construct useful new integration rules.

Let me know what you think in the comments!