Benchmarking the Unified Substitution Method (USM) Against Mathematica Integrate, Part 3

This post documents the updated benchmark suite for the Unified Substitution Method (USM), a branch-consistent integration framework for algebraic and half-angle integrands involving quadratic radicals. The method is described in the USM article, A Unified Substitution Method for Integration, arXiv:2505.03754, together with the supplementary appendix on half-angle identities and continuation of USM antiderivatives.

The full Wolfram benchmark code is available here: Download the Wolfram benchmark code PDF.

The key point is that USM is not merely a real-variable substitution trick. It is built on the principal complex branches inherited from the principal logarithm and the principal square root. In the USM article, the transforms are derived from identities involving exponentials of principal inverse trigonometric functions such as

e^(± i ArcCos[y])
e^(± i ArcSec[y])

These identities lead to explicit substitutions that reduce large classes of radical and half-angle integrands to rational functions of one parameter. Depending on the transform, that parameter is usually denoted by t, r, or s. These parameters act as branch-aware conformal parametrizations on the relevant analytic components of the complex plane, allowing the same rational primitive to be continued across domains with the correct principal-branch behavior.

What Type of Integrands Are Being Benchmarked?

The benchmark suite targets integrands that combine polynomial terms, quadratic radicals, and inverse-trigonometric half-angle expressions. These are precisely the types of integrands where a general-purpose symbolic integrator can suffer from branch ambiguity, expression swell, or non-termination.

Typical benchmark integrands have the form

1 / (c0 + c1*x^n0 + c2*x^n1*RadicalTerm + c3*x^n2*HalfAngleTerm)

where the radical or half-angle component belongs to one of the five USM transform geometries. For example, Transform 2 tests integrands involving

Sqrt[(x + b)^2 - a^2]

and

Sqrt[(x + b - a)/(x + b + a)]

while Transform 5 tests integrands involving

Sqrt[a^2 - (x + b)^2]

and

Sqrt[(a + b + x)/(a - b - x)]

These expressions are not artificial stress tests in the sense of being random noise. They are structured algebraic integrands generated from the exact radical geometries that the USM transforms are designed to rationalize.

The Five Geometries

Let

y = (x + b)/a

with a > 0. The topological behavior of the integrand depends on the core radical geometry.

Transform	Core Geometry	Parameter	Complex Region on the Real Axis
Transforms 1 and 2	Sqrt[y^2 - 1]	t	Complex on the interior -1 < y < 1
Transform 3	Sqrt[y^2 + 1]	s	Never complex on the real line
Transforms 4 and 5	Sqrt[1 - y^2]	r	Complex on the exteriors y > 1 and y < -1

This distinction is essential. Transforms 1 and 2 are naturally real on the exterior intervals and complex on the interior. Transforms 4 and 5 are naturally real on the bounded circular interval and complex on the exterior intervals. Transform 3 is different: since

y^2 + 1 >= 1

for every real y, the principal square root stays positive and real on the whole real axis. The parameter

s = y + Sqrt[y^2 + 1]

is real, positive, and never crosses zero. This is why Transform 3 requires no separate componentwise continuation theorem. A single universal template works globally.

Why the Benchmark Does Not Use Simplify[D[F[x], x] - I[x]]

For small examples, one might verify an antiderivative by asking Mathematica to simplify

Simplify[D[anti, x] - integrand]

to zero. That approach is not reliable for this benchmark class. The antiderivatives can contain nested radicals, branch-dependent logarithms, and large RootSum objects. In those cases, symbolic simplification can hang, return unevaluated, or spend more time proving the result than computing the integral itself.

Instead, the benchmark engine validates each antiderivative numerically at branch-safe test points. If F(x) is the returned antiderivative and I(x) is the original integrand, the engine checks whether

F'(x0) ≈ I(x0)

at carefully selected points x0 in the appropriate domain. The values are evaluated to high precision, microscopic numerical noise is removed with Chop, and the relative error must be less than

10^-4

for the branch to be marked as correct.

Why Fractional Test Points Matter

The random generator often creates examples with integer branch boundaries. Testing exactly at an integer can accidentally hit a singularity, a branch point, or an undefined value. The benchmark avoids this by stepping away from each boundary using fractional offsets.

For example, consider

I(x) = 1 / Sqrt[x^2 - 4]

The branch boundaries are x = -2 and x = 2. Instead of testing at the boundaries, the engine uses points such as:

• Right exterior: x = 2 + 2.718 = 4.718
• Interior: x = 0 + 2*0.314 = 0.628
• Left exterior: x = -2 - 2.718 = -4.718

This makes the verification branch-aware. Each returned antiderivative is tested where it is supposed to be valid, rather than at a point where the principal branch is changing.

What the Updated Tables Show

The updated benchmark tables report more than raw timing. They explicitly display:

• the generated integrand,
• the USM antiderivative,
• the interval or branch domain where the formula is used,
• the fractional test points,
• whether the derivative check passed,
• cold execution time,
• warm execution time,
• and byte count of the output expression.

The byte count is important because it measures expression bloat. A system may technically find an antiderivative but return a huge expression that is difficult to inspect, differentiate, simplify, or reuse. USM often avoids this by rationalizing the integrand first and then back-substituting through a controlled branch formula.

Where Mathematica Integrate Struggled

The screenshot below shows a native Mathematica Integrate benchmark table for a difficult family of algebraic radical integrands.

In this run, Mathematica returned several antiderivatives as unevaluated integral expressions. One case produced a large antiderivative with a byte count of 55296, and another case hit the configured 15-second timeout. The table also reports the aggregate timing:

MATH AVERAGES | Cold: 2.2397 s | Warm: 0.1952 s

This is exactly the type of behavior the benchmark is designed to expose. The issue is not that Mathematica cannot integrate simple radicals. The issue is that when polynomial factors, shifted quadratic radicals, quotient radicals, and branch-sensitive algebraic structure are mixed together, the general-purpose integrator may return an unevaluated result, a very large expression, or a timeout.

Why USM Is Different

USM attacks the structure directly. Instead of asking a general integrator to discover the right substitution, the benchmark applies the appropriate USM transform first. The radical and half-angle terms are converted into rational functions of a single parameter.

For Transforms 1 and 2, the parameter t satisfies reconstruction formulas of the form

y = (1/2) (t + 1/t)

x = (a/2) (t + 1/t) - b

dx = a*(t^2 - 1)/(2*t^2) dt

For Transforms 4 and 5, the parameter r satisfies

y = 2*r/(1 + r^2)

x = 2*a*r/(1 + r^2) - b

dx = a*2*(1 - r^2)/(1 + r^2)^2 dr

For Transform 3, the hyperbolic parameter s satisfies

y = (s - 1/s)/2

x = a*(s^2 - 1)/(2*s) - b

Sqrt[(x + b)^2 + a^2] = a*(s^2 + 1)/(2*s)

dx = a*(s^2 + 1)/(2*s^2) ds

After this substitution, the transformed integrand is rational in the parameter. Mathematica is then used for the rational integration step, which is a much more controlled problem than the original branch-heavy expression.

Dynamic Formula Merging

The updated benchmark also prevents redundant output. Some transforms require separate formulas on different real components. However, when the analytical continuation of the USM parameter aligns with Mathematica’s principal root across a boundary, the code checks whether the two returned antiderivatives are identical. If they are, the table collapses them into a single formula.

This is why the output sometimes displays one global antiderivative and sometimes displays separate formulas for the left, interior, or right branches. The goal is not to force piecewise output unnecessarily. The goal is to preserve branch correctness while avoiding redundant formulas.

Changing the Sample Size

The scripts are standalone. To run larger or smaller batches, change the final iterator count in the integrand generation block and then change the matching iterator in the results block.

For example:

SeedRandom[102];

integrands2 = Table[
  Module[
    {
      a = RandomInteger[{1, 15}],
      b = RandomInteger[{-10, 10}],
      c0, c1, c2, c3, n0, n1, n2, expr
    },

    c0 = RandomChoice[DeleteCases[Range[-5, 5], 0]];
    c1 = RandomChoice[DeleteCases[Range[-5, 5], 0]];
    c2 = RandomChoice[DeleteCases[Range[-5, 5], 0]];
    c3 = RandomChoice[DeleteCases[Range[-5, 5], 0]];

    n0 = RandomInteger[{1, 3}];
    n1 = RandomInteger[{0, 2}];
    n2 = RandomInteger[{0, 2}];

    expr =
      c0 +
      c1*x^n0 +
      c2*x^n1*Sqrt[Expand[(x + b)^2 - a^2]] +
      c3*x^n2*Sqrt[(x + b - a)/(x + b + a)];

    1/expr
  ],
  {10}
];  (* Change this number *)

Then update the matching results iterator:

results2 = Table[
  Module[
    {
      eq = integrands2[[i]],
      antiUSM, coldUSM, warmUSM, bytesUSM
    },

    (* benchmark computation here *)

  ],
  {i, 10}
];  (* Change this number to match *)

Conclusion

These benchmarks are not only measuring speed. They are measuring whether a branch-aware substitution framework can turn difficult symbolic integration problems into predictable rational integrations.

For the integrand families tested here, the central advantage of USM is structural. The method knows in advance which radical geometry is present, which parameter should be used, where the branch transitions occur, and how to back-substitute using the principal complex branches. That is why the benchmark records not only cold and warm timing, but also domains, numerical derivative verification, and expression byte count.

The attached Mathematica screenshot shows why this matters: native Integrate can return unevaluated integrals, very large expressions, or timeouts on exactly the kind of mixed algebraic radical integrands that USM is designed to rationalize.