Skip to content

Instantly share code, notes, and snippets.

@asmeurer

asmeurer/RISCH_ALGEBRAIC_PLAN.md

Created August 15, 2026 18:20
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save asmeurer/5ad380f4ed1cdd573e6328ffd0d8a947 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save asmeurer/5ad380f4ed1cdd573e6328ffd0d8a947 to your computer and use it in GitHub Desktop.
Download ZIP
Plan: branch-correct algebraic integration in SymPy's Risch code, and correctness testing against the Rubi corpus
Raw
RISCH_ALGEBRAIC_PLAN.md

Plan: branch-correct algebraic integration, and correctness testing

against the Rubi corpus

Companion to RISCH_PLAN.md (the transcendental Risch plan) and BRONSTEIN_ERRATA.md. This file covers the experimental algebraic track on branch risch-algebraic and is written to be picked up by a fresh session.

Status: research and prototyping done, nothing implemented. Two work items are queued, in this order:

  1. Correctness testing against the Rubi corpus' expected antiderivatives (Aaron: "Correctness is the most important thing, so let's be thorough").
  2. Branch-correct radicand normalization via the signum algorithm.

Item 1 first: it is the instrument that measures item 2, and it will also re-audit everything already landed.

1. Where things stand

Branches (SHAs shift under rebases -- match on commit subjects)

  • risch-gaps -- Phase 0, PR #30180. Tip at time of writing d01ca185c7.
  • risch-rde-cancellation -- Phases 1-4 plus fixes; stacked on risch-gaps. Tip bf88623d8e ("Rewrite all computable residue terms via log_to_real() in the Risch code"). The last three commits are Aaron's Rioboo LogToReal work (real arc-tangents for complex-conjugate residue pairs) -- relevant here, since output form and branch behavior interact.
  • risch-algebraic -- the experimental exp-log-tower representation of radicals; stacked on risch-rde-cancellation. Tip a1c2226052. Nine commits: Aaron's two proof-of-concept commits, repairs, the algebraic=False gate, the acceptance filter, the radicand normalization, and regression tests.

No PRs opened for the latter two branches. Aaron's rule: keep the branches sharing identical commits so merged work disappears from the others; rebase only when necessary.

What the algebraic mode is

risch_integrate(f, x, algebraic=True) represents radicals as exp-log towers (sqrt(x) as exp(log(x)/2)), runs the transcendental machinery over the (non-transcendental) tower, and degrades every nonelementary conclusion to a plain Integral, since those proofs assume transcendence. Off by default: ungated tower building sent the integrals test suite from 90 s to a 10-minute timeout, because integrate() attempts Risch automatically.

Results over the Rubi corpus (~16,800 attemptable radical cases): about 1,357 integrals solved that non-Risch integrate() cannot do, 2,771 solved by both, zero false nonelementary claims. Timing: solved cases have a 0.09-0.26 s median; roughly 2% of cases exceed 300 s, all in bounded-arithmetic sites (no loops). Full details, per-case tables and the run log-book: https://gist.github.com/asmeurer/b4b8ceb7c364566f5e7a3d07ce133300

The acceptance filter (already implemented)

_nontrans_accept() in sympy/integrals/risch.py, called from integrate_primitive() and integrate_hyperexponential() when not DE.transcendental. It accepts a candidate only if

  1. f == D(elem) + D(residue part) + i holds as a formal rational-function identity in the tower generators (decided by cancel(), on the polynomial representation -- never on the final Expr, per Aaron: sympy's Expr-level simplification cannot be trusted or afforded for algebraic expressions), and
  2. no denominator (of elem, i, or the residue terms' root polynomials and logarithm arguments) lies in the kernel of the evaluation map, decided by reduction modulo the tower's power relations t**q == u**p.

Check 2 was added after a backward-constructed test exposed a real hole; the corpus never triggers either check (400-case instrumented sample), so both are covered by constructed tests only.

Known limitation, and the reason for item 2 below: the filter certifies the candidate against the tower image of the integrand as the tower builder rewrote it. Any rewriting applied before the tower is built is outside its guarantee.

2. The open correctness problem

5779866b43 ("Factor radicands before erecting algebraic tower generators") factors each radicand and distributes the fractional power over the factors, so that content like a perfect square does not become a spurious generator. Measured benefit: +387 SOLVED-NEW cases corpus-wide (23.2% -> 25.6% overall solve rate).

This is an algebraic equivalence only where the factors are positive. sqrt((x+1)**2*(x+2)) -> (x+1)*sqrt(x+2) drops an absolute value and has the wrong sign for x < -1. Confirmed numerically on the current branch:

risch_integrate(sqrt(x**2 + 2*x + 1)/x, x, algebraic=True)
    -> sqrt(x**2 + 2*x + 1) + log(x) - 1
at x = -3:  D(answer) = -1.3333  but  f(-3) = -0.6667

and Timofeev's (x**3 - 5*x**2 + 3*x + 9)**(-2/3) (radicand (x-3)**2*(x+1), one of the nine chapter-0 headline solves) is wrong on all of x < 3, returning complex values where the integrand is real.

Scope of the damage:

  • The core machinery is not affected. exp(log(u)/2) is identically the principal sqrt(u) on the cut plane, so the tower representation substitutes nothing; core-path answers spot-checked at negative and complex points are faithful. Only the pre-tower rewriting is at fault.
  • The regression test that pins the normalization uses Symbol('y', positive=True), where the collapse is sound. The corpus failures are all on assumption-less symbols.
  • Continuity is a separate, milder matter: the core path inherits sympy's principal-branch conventions (log jumps), the same as stock integrate(). Aaron's new Rioboo work improves output form here.

Interim options considered and rejected: gating the split on provable nonnegativity is sound but loses essentially all corpus gains (Rubi's symbols carry no assumptions); keeping generic-only validity contradicts sympy's assumption discipline. The signum algorithm below is strictly better than both.

3. Item 1: correctness testing against the corpus' expected answers

The corpus is Upabjojr/rubi-integration-test-suite, cloned at ~/Documents/Python/sympy/rubi-integration-test-suite (a sibling of the sympy checkout; not to be checked into sympy). Each RubiTestSuiteCase carries integrand, variable, num_steps and integral -- the expected antiderivative, which we have never used. The runner is risch_test_suite_runner.py in that clone, also proposed upstream as PR #1 there (opened as Claude, per Aaron).

Goal: use integral as an oracle to find answers that are potentially wrong, then investigate each. Our answers legitimately differ in form, so a mismatch is a signal to investigate, not a verdict.

3.1 Comparison ladder

Apply in order; stop at the first conclusive outcome.

  1. Difference is constant. ours - expected should be constant. Test by differentiating and simplifying to zero, on the polynomial representation where possible.
  2. Derivative check. D(ours) - f == 0. Note this is what the internal filter already does at tower level; at Expr level it is unreliable for radicals (during this session a correct answer was flagged "WRONG" purely because Expr-level zero-testing of a radical identity failed -- exactly the reason Aaron insisted the internal check stay on the polynomial representation).
  3. Numerical evaluation at many points -- the decisive test for this class of bug, and the one that catches branch errors that every symbolic check misses. Requirements:
    • sample both signs of every factor that appears under a radical, and points on either side of each real root of those factors (the branch bug is invisible if you only sample where everything is positive);
    • sample complex points as well as real ones;
    • compare D(ours) against f pointwise, and separately check ours - expected is constant across points within a connected region;
    • freeze sign() factors locally before differentiating (they are locally constant; naive differentiation produces DiracDelta).

3.2 The symbolic-constants caution (Aaron's point)

Rubi's a, b, c, ... are pattern variables with no assumptions, whereas real sympy usage substitutes actual numbers. Both directions matter and both must be tested:

  • Instantiation testing. Replace the symbolic constants with random concrete values -- rational and irrational, positive and negative, and some complex -- then re-run the comparison ladder. This is the strongest available oracle: it turns a generically-valid formula into a checkable one, and it is precisely how a domain-dependent error (a dropped absolute value, a division by a constant that vanishes) becomes visible. Skip instantiations that make the integrand degenerate (a denominator identically zero).
  • Assumption sensitivity. Run each case with the constants assumption-free and again with positive=True, and diff the outcomes. Divergence is informative in both directions: it flags results that silently depend on assumptions (our normalization is sound under positive=True and wrong without), and it flags cases where sympy refuses to proceed without assumptions.
  • Special values. Related to the Piecewise/conds work already documented in RISCH_PLAN.md item 6: a result obtained by dividing by a symbolic constant is valid only generically. Instantiating at the vanishing values of those divisors is the test for it.

3.3 Mismatch taxonomy

Every mismatch should be classified, not just counted:

  • our answer genuinely wrong (branch, sign, dropped condition);
  • our answer correct, different form (constant of integration, algebraically equal, or log vs atan forms -- especially now that Rioboo rewriting is in);
  • both correct on different domains (Rubi's conventions are real- oriented; sympy's are principal-branch);
  • the corpus' own answer questionable (Rubi answers are not immune to branch issues, and Rubi's sqrt conventions differ from sympy's);
  • integrand degenerate under the chosen instantiation.

3.4 Deliverables for item 1

  • Extend the runner with an expected-comparison mode (keep it environment-variable driven, as with RISCH_RESULTS, RISCH_MODE, RISCH_HANDLE_FIRST, so the CLI stays compatible with the upstream runner) and push it to the same PR branch.
  • Run it over every case the algebraic mode currently solves (~4,100 SOLVED-NEW + SOLVED-both), and over the transcendental chapters as a control -- this re-audits everything already landed, including the pre-normalization solves.
  • Publish the mismatch table (SymPy expression, our answer, expected answer, classification) as a new page in the run gist.
  • File a regression test in sympy for every confirmed wrong answer.

4. Item 2: branch-correct normalization via the signum algorithm

4.1 Sources (both in ~/Dropbox/papers/symbolic-computation/)

  • Jeffrey 1993, Integration to obtain expressions valid on domains of maximum extent, ISSAC 93. §2 has the log-combining theory, including the rule for fractional coefficients a*ln f1 + b*ln f2 -> (m/n)*ln(f1**p * f2**q) (needed because combining under a fractional power reintroduces discontinuities). §5 works ∫sqrt(x**(2/3) + x**(4/3)), whose correct antiderivative carries sgn**(5/3) factors -- our failure shape.
  • Jeffrey, Labahn, von Mohrenschildt & Rich, Integration of the signum, piecewise and related functions. §5 gives the complete algorithm; Theorem 5 gives the jump-correction formula; Example 1 is the case reproduced by the prototype below.
  • Background: numerical-analysis/Kahan 1986 - Branch Cuts for Complex Elementary Functions; theses/von Mohrenschildt 1994 - Symbolic Solutions of Discontinuous Differential Equations.

4.2 The algorithm, as it applies here

Rewrite sqrt(w**2 * v) as s*w*sqrt(v) where s = sgn(w) is introduced as a symbolic constant, integrate as usual, then

  • substitute s -> sgn(w) in the result, and
  • add a jump correction -J_k * sgn(x - x_k) at each breakpoint x_k (each real root or pole of w), where J_k = (G(x_k, s=+1) - G(x_k, s=-1))/2.

The corrections restore continuity across the breakpoints, giving answers on Jeffrey's "domain of maximum extent" -- better than merely correct. The approach fits this machinery unusually well because the signs ride through as ordinary symbolic constants, which the corpus runs already showed the towers handle at rates comparable to concrete coefficients.

4.3 Prototype evidence (verified this session)

signum_proto2.py, attached to the run gist. Results:

  • Jeffrey's Example 1 reproduced exactly: 3*x**2*sqrt(1 + 1/x**2) -> sgn(x)*((1 + x**2)**(3/2) - 1), with J = 1.
  • sqrt((x+1)**2*(x+2)): J = -4/15; derivative correct at x = -1.9, -1.5, -0.5, 3 (the first two are where the current code is wrong), and the correction takes the discontinuity at x = -1 from 0.533 to 0.
  • sqrt(x**2 + 2*x + 1)/x: correct at x = -3 and -1.5, the exact failing points.

4.4 Implementation sketch

  1. Replace the current unconditional split in _rewrite_exps_pows() with the sign-carrying rewrite; record (s, w) pairs on the DifferentialExtension (a new slot, like backsubs).
  2. After integration, in risch_integrate(), substitute the signs back and apply the jump corrections. Keep this outside the tower machinery -- it is a post-processing pass on the final result.
  3. Breakpoints: real roots of numerator and denominator of each w.
  4. Where the sign of w is provable (w.is_nonnegative), skip the whole apparatus and split unconditionally, as now.

4.5 Open problems to solve during implementation

  • Complex jumps. When a breakpoint coincides with a singularity of the integrand, J comes out complex (our sqrt(x**2+2*x+1)/x case gives J = -1 + I*pi). Theorem 7 of the signum paper covers the integrable-singularity case; decide what to return when the jump is genuinely infinite.
  • Odd-order radicals. Cube roots need the sgn**(2/3) treatment of Jeffrey 1993 §5 -- this is the Timofeev case, and sympy's principal-branch cube root of a negative real is complex while Rubi's is real. Decide the convention explicitly.
  • Non-polynomial sign arguments. Breakpoint detection needs the real roots of arbitrary w; restrict to what roots/real_roots can do and fall back to leaving the case unsplit.
  • Interaction with the acceptance filter. With signs as symbolic constants, the filter's kernel test must not treat s as an ordinary constant that could vanish -- s**2 == 1 should be part of the relation set.
  • Output size. sgn factors multiply through; consider collecting them (sgn(w)*G rather than distributing) for readable answers.

5. Working notes for a new session

Infrastructure

Conventions learned the hard way

  • Never point long-running jobs at the live checkout. Aaron works in it and switches branches; a mid-run checkout silently invalidated a whole timing study (every call raised TypeError in 0.00 s, which looked like instant success). Use git worktree add <scratchpad>/wt-<name> <branch> and run against that. Remove worktrees when done so Aaron can check the branch out.
  • Timing claims come from serial runs only; parallel sweeps inflate timeouts through contention.
  • Sanity-check measurements: 0.00 s medians mean the code did not run.
  • Tests: python -m pytest -p no:pudb <files> -q (the pytest-pudb plugin is broken). Lint: ruff check sympy/integrals/.
  • Commit style: no comments about past code state or session context (that belongs in the commit message); known sympy idioms need no comment; keep fixes minimal; fix bugs at their source rather than working around them in a caller. No Claude-Session trailers; keep Co-Authored-By.
  • Commits are SSH-signed through Secretive and need Aaron's approval touch; a persistent monitor that retries commit-and-push is the smoothest way to handle it.

Bugs this corpus has already surfaced (all fixed, with tests)

Five in sympy master: the spde() infinite loop for SymPy Integer degree bounds; malformed Polys from is_log_deriv_k_t_radical_in_field(); ratint_logpart()'s PolynomialError for radical coefficients (issue #26502, partially fixed -- octic cases still crash); PolyMatrix scalar multiplication with a ground-domain fallback; plus the bound_degree Poly-vs-int family from Phase 0. Two in the algebraic branch itself, and one hole in the acceptance filter. Expect item 1 to find more.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment