Skip to content

A repeated quadratic denominator is integrated (#180) #1155

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Sep 3, 2026
Merged

A repeated quadratic denominator is integrated (#180) #1155

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 55 additions & 0 deletions BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,61 @@ read first.
| | `"sum(k, k, 1, n)".ToEntity().Simplify()`, and every summation whose body is a polynomial in the index | `sum(k, k, 1, n)` — carried | `piecewise((n + n ^ 2) / 2 provided n >= 0, 0)` |
| | `"sum(k, k, 1, 100000)".ToEntity().Simplify()`, and every concrete range past a hundred terms | `sum(k, k, 1, 100000)` — carried | `5000050000` |
| | `"product(k, k, 1, n)".ToEntity().Simplify()`, and every product whose body is a monomial in the index | `product(k, k, 1, n)` — carried | `piecewise(n! provided n >= 1, 1)` |
| | `"1/(x^2 + 1)^2".Integrate("x")`, and every proper rational function over a repeated irreducible quadratic | `integral(1 / (1 + x ^ 2) ^ 2, x)` — left unevaluated, after 2.9 s | `arctan(x) / 2 + C + 1/2 * x / (x ^ 2 + 1)`, in 0.11 s |
| | `"x^2/(x^2 + 2)^2".Integrate("x")`, and every polynomial numerator over one | `integral(x ^ 2 / (x ^ 2 + 2) ^ 2, x)` — left unevaluated, after 93 s | the antiderivative, in 0.57 s |
| | `"1/(x^2 - 1)^2".Integrate("x")`, and every repeated quadratic whose roots are real | `C - ln(x + -1) / 4 + ln(1 + 2 * x + x ^ 2) / 8 + -1/2 * x / (x ^ 2 + -1)` | `C - ln(1 + (-2) / (x + 1)) / 4 + -1/2 * x / (x ^ 2 + -1)` — the same value, one logarithm rather than two |

### A repeated quadratic denominator is integrated

`1/(x^2 + 1)^2` had no antiderivative, while `1/(x^2 - 1)^2` and `1/(x^2 + 2x + 1)^2` both did — a
denominator with real roots comes apart into linear factors and never reached the gap. What was
missing was the irreducible case
([#180](https://github.com/asc-community/AngouriMath/issues/180)).

```
"1/(x^2 + 1)^2".Integrate("x") was integral(1 / (1 + x ^ 2) ^ 2, x)
is arctan(x) / 2 + C + 1/2 * x / (x ^ 2 + 1)
"x^2/(x^2 + 2)^2".Integrate("x") was integral(x ^ 2 / (x ^ 2 + 2) ^ 2, x)
is (sqrt(2) * arctan(sqrt(2) * x / 2) * 2 + (-4) * x / (x ^ 2 + 2)) / 8
+ C provided not x ^ 2 + 2 = 0
```

The condition on that second one is the library's own and it is true: the antiderivative has no
value where the denominator vanishes. It is attached to some of these answers and not others, and I
did not pin down what decides which — recorded as observed rather than explained.

Writing `Q` for the quadratic, `u` for its derivative `2ax + b` and `D` for `4ac - b^2`, the
identity `u^2 = 4aQ - D` turns the derivative of `u/Q^(m-1)` into an equation for the integral of
`1/Q^m`, which is unrolled down to the single power the table already answered. It is an identity in
`a`, `b` and `c` rather than a fact about signs, so one line serves both signs of the discriminant;
`D = 0` is the sole exclusion, and the division by `D` is why. A numerator of higher degree is
divided by `Q` first.

**A denominator with real roots now takes this route rather than partial fractions.** The value is
unchanged and the form is shorter — one logarithm where there were two — which is why the
node-count metric that ranks candidates selects it. No test pinned the old spelling.

**The cost, measured rather than estimated.** Where the integrand is rational this is between 26 and
165 times faster, because `TryStandardIntegrals` runs before every search and a shape answered
outright never enters the candidate exploration. Where the integrand is a transcendental function
over a repeated quadratic — `sin(x)/(x^2 + 1)^2`, `e^x/(x^2 + 1)^2`, `ln(x)/(x^2 + 1)^2`, none of
which has an elementary antiderivative — the same search now takes **3.4 to 10.8 times longer to
conclude nothing**, because the sub-integral it used to fail on immediately now succeeds and
integration by parts carries on from a larger expression. Both arms measured on one machine, master
and branch, the same way:

| | master | now |
|---|---|---|
| `1/(x^2 + 1)^2` | 2879 ms, unevaluated | **111 ms, answered** |
| `x^2/(x^2 + 2)^2` | 92 990 ms, unevaluated | **565 ms, answered** |
| `1/(x^4 + 1)^2` | 3811 ms | 3870 ms — the base is not quadratic, and nothing changes |
| `ln(x)/(x^2 + 1)^2` | 15 034 ms, unevaluated | 51 072 ms, unevaluated |
| `sin(x)/(x^2 + 1)^2` | 26 492 ms, unevaluated | 114 851 ms, unevaluated |
| `e^x/(x^2 + 1)^2` | 7944 ms, unevaluated | 85 683 ms, unevaluated |

That last group is a growth in integration by parts rather than in this rule, which those integrands
never match; it is filed separately rather than fixed here, since bounding by-parts is a change to
machinery every integral goes through.

### A nested radical comes apart

Expand Down
3 changes: 3 additions & 0 deletions Sources/.editorconfig
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,3 +189,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed

[Tests/UnitTests/PatternsTest/DenestRadicalTest.cs]
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n

[Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs]
file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
177 changes: 177 additions & 0 deletions Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,18 @@ _ when TryReadOverSquareTimesRoot(expr, x, out var overFactor, out var overRadic
&& a.Evaled is Entity.Number.Complex { IsZero: false }
=> IntegrateLinearOverQuadratic(p, q, a, b, c, denominator, x),

// ∫ N(x)/(ax^2 + bx + c)^n dx for a whole n of two or more and a proper fraction.
// After the single-power arms above, which already read a denominator like
// (x + 1)^2 -- expanded, it is a quadratic, so the perfect-square arm there answers
// it and nothing about that changes.
Entity.Divf(var numerator, Entity.Powf(var repeated, Entity.Number.Integer repetitions)) when
repetitions.EInteger.CanFitInInt32()
&& repetitions.EInteger.ToInt32Checked() is var n and >= 2
&& TreeAnalyzer.TryGetPolyQuadratic(repeated, x, out var qa, out var qb, out var qc)
&& (!numerator.ContainsNode(x) || qa.Evaled is Entity.Number.Complex { IsZero: false })
&& IntegrateRationalOverPowerOfQuadratic(numerator, qa, qb, qc, repeated, n, x) is { } repeatedAnswer
=> repeatedAnswer,

// ∫ (px + q)/(bx + c) dx, the same rewrite one degree down: the quotient is the
// constant p/b plus a remainder over the divisor. x/(x + 1) had no antiderivative.
Entity.Divf(var numerator, var denominator) when
Expand Down Expand Up @@ -529,5 +541,170 @@ private static Entity IntegrateRationalQuadratic(Entity numerator, Entity a, Ent
new Entity.Providedf(lnCase, discriminant < 0)
]).InnerSimplified;
}

/// <summary>
/// ∫ k / (a x^2 + b x + c)^n dx for a whole n of two or more, by the reduction that
/// takes one power off the denominator at a time.
/// </summary>
/// <remarks>
/// <para>
/// Write <c>Q</c> for the quadratic, <c>u</c> for its derivative <c>2ax + b</c> and
/// <c>D</c> for <c>4ac - b^2</c>. Then <c>u^2</c> is <c>4aQ - D</c>, so differentiating
/// <c>u / Q^(m-1)</c> gives <c>2a(3 - 2m)/Q^(m-1) + (m-1)D/Q^m</c>, and reading that as
/// an equation for the integral of <c>1/Q^m</c> leaves
/// <code>
/// J_m = (u / Q^(m-1) + 2a(2m - 3) J_(m-1)) / ((m - 1) D)
/// </code>
/// which is unrolled down to <c>J_1</c>, the case the table above already answers. It is
/// an identity in a, b and c rather than a fact about signs, so one line serves a
/// positive and a negative discriminant alike; <c>D = 0</c> is the only exclusion, and
/// the division by <c>D</c> is why.
/// </para>
/// <para>
/// The two shapes the reduction cannot speak for get their own arms, and each is
/// elementary: with <c>a = 0</c> the denominator is a linear power, and with
/// <c>D = 0</c> the quadratic is <c>u^2/(4a)</c> so the whole integrand is a power of
/// <c>u</c>. Both are already answered for a numeric quadratic — by the power rule and
/// by partial fractions — and the arms are here for a symbolic one, where the
/// discriminant's sign is not known and dropping an arm is not available.
/// </para>
/// <para>
/// What was missing was the irreducible case. <c>1/(x^2 - 1)^2</c> had an antiderivative
/// and <c>1/(x^2 + 2x + 1)^2</c> had one, because a denominator with real roots comes
/// apart into linear factors and never reaches here; <c>1/(x^2 + 1)^2</c> had none.
/// https://github.com/asc-community/AngouriMath/issues/180
/// </para>
/// </remarks>
private static Entity IntegrateOverPowerOfQuadratic(
Entity numerator, Entity a, Entity b, Entity c, int power, Entity.Variable x)
{
var quadratic = a * x * x + b * x + c;
var derivative = 2 * a * x + b;
var discriminant = 4 * a * c - b * b;

// a = 0: k/(bx + c)^n, an ordinary power -- and the constant k/c^n when b is zero
// too, where writing the power would divide by it. Same reason as the rule above.
var linearCase = TreeAnalyzer.IsZero(b)
? numerator * x / MathS.Pow(c, power)
: numerator * MathS.Pow(b * x + c, 1 - power) / (b * (1 - power));

// D = 0: the quadratic is u^2/(4a), so the integrand is k(4a)^n / u^(2n).
var perfectSquareCase =
numerator * MathS.Pow(4 * a, power) * MathS.Pow(derivative, 1 - 2 * power)
/ (2 * a * (1 - 2 * power));

var sqrtDiscriminant = MathS.Sqrt(discriminant);
var sqrtNegDiscriminant = MathS.Sqrt(-discriminant);

return MathS.Piecewise([
new Entity.Providedf(linearCase, a.EqualTo(0)),
new Entity.Providedf(perfectSquareCase, discriminant.EqualTo(0)),
new Entity.Providedf(
Reduce(2 * MathS.Arctan(derivative / sqrtDiscriminant) / sqrtDiscriminant),
discriminant > 0),
new Entity.Providedf(
Reduce(AntiderivativeLog((derivative - sqrtNegDiscriminant) / (derivative + sqrtNegDiscriminant))
/ sqrtNegDiscriminant),
discriminant < 0)
]).InnerSimplified;

Entity Reduce(Entity firstPower)
{
var accumulated = firstPower;
for (var m = 2; m <= power; m++)
accumulated =
(derivative / MathS.Pow(quadratic, m - 1) + 2 * a * (2 * m - 3) * accumulated)
/ ((m - 1) * discriminant);
return numerator * accumulated;
}
}

/// <summary>
/// ∫ (px + q)/(a x^2 + b x + c)^n dx, by the same rewrite of the numerator the single
/// power uses: <c>px + q = (p/2a)(2ax + b) + (q - pb/2a)</c>. The first part is the
/// denominator's own derivative over a power of it, which integrates as a power; the
/// second is the constant-numerator case above.
/// </summary>
/// <remarks>
/// A numerator that is already a multiple of that derivative leaves nothing of the
/// second part, and it is dropped rather than multiplied by zero. Multiplying is not
/// harmless: the reduction is a quotient by the quadratic, so <c>0 * (2x/(x^2 + 1))</c>
/// is <c>0 provided not 1 + x^2 = 0</c> and not <c>0</c> — correctly, since the factor
/// has no value there. Carrying that condition would attach it to the answer for
/// <c>x/(x^2 + 1)^2</c>, which had none before and is owed none.
/// </remarks>
private static Entity IntegrateLinearOverPowerOfQuadratic(
Entity p, Entity q, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x)
{
var alongTheDerivative = p / (2 * a) * MathS.Pow(quadratic, 1 - power) / (1 - power);
var constantPart = q - p * b / (2 * a);
return TreeAnalyzer.IsZero(constantPart)
? alongTheDerivative
: alongTheDerivative + IntegrateOverPowerOfQuadratic(constantPart, a, b, c, power, x);
}

/// <summary>
/// ∫ N(x)/(a x^2 + b x + c)^n dx for a numerator of any degree, or <see langword="null"/>
/// where the fraction is improper and this is not the rule for it.
/// </summary>
/// <remarks>
/// <para>
/// Dividing the numerator by the quadratic writes <c>N = qQ + r</c> with <c>r</c> linear
/// at worst, so <c>N/Q^n</c> is <c>q/Q^(n-1) + r/Q^n</c> — one power off the denominator
/// and two degrees off the numerator at each step, which ends. The remainder is answered
/// by the two helpers above and the quotient goes round again.
/// </para>
/// <para>
/// <b>Improper fractions are declined rather than divided out here.</b> The recursion
/// bottoms out at a single power, and a numerator still of degree two or more there is
/// exactly the improper case, which
/// <see cref="IndefiniteIntegralSolver.SolveByPartialFractions"/> already opens with
/// <see cref="TreeAnalyzer.PolynomialLongDivision"/>. Declining hands it there instead of
/// answering it twice, in two shapes, from two places.
/// </para>
/// <para>
/// This is not only reach. <c>TryStandardIntegrals</c> runs before every search, so a
/// shape answered here never reaches the candidate exploration that re-integrates
/// rewritten forms — which is where this family's cost has always been. Both
/// <c>x^2/(x^2 + 2)^2</c> and <c>(x^2 + 1)/(x^2 + 2)^2</c> were tens of seconds spent to
/// return the integral unevaluated.
/// </para>
/// </remarks>
private static Entity? IntegrateRationalOverPowerOfQuadratic(
Entity numerator, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x)
{
if (!numerator.ContainsNode(x))
return power == 1
? IntegrateRationalQuadratic(numerator, a, b, c, x)
: IntegrateOverPowerOfQuadratic(numerator, a, b, c, power, x);

if (TreeAnalyzer.TryGetPolyLinear(numerator, x, out var p, out var q))
return power == 1
? IntegrateLinearOverQuadratic(p, q, a, b, c, quadratic, x)
: IntegrateLinearOverPowerOfQuadratic(p, q, a, b, c, quadratic, power, x);

if (power < 2)
return null;

var division = TreeAnalyzer.PolynomialLongDivision(numerator, quadratic);
if (division is null)
return null;

// The quotient is what is wanted; the remainder that comes back with it is already
// divided by the divisor -- x^2 over x^2 + 2 answers (1, -2/(x^2 + 2)) -- so taking
// it as a polynomial and dividing it again is a numerator that is not one, and the
// next division declines. Subtracting gives the polynomial remainder whatever shape
// the division chose to report.
// Reduced before either is used: the division reports its quotient built up term by
// term, so x^2 over x^2 + 2 comes back as `0 + 1` rather than `1`, and the arms below
// ask whether a numerator is constant.
var quotient = division.Value.Divided.InnerSimplified;
var remainder = (numerator - quotient * quadratic).InnerSimplified;

return IntegrateRationalOverPowerOfQuadratic(remainder, a, b, c, quadratic, power, x) is { } head
&& IntegrateRationalOverPowerOfQuadratic(quotient, a, b, c, quadratic, power - 1, x) is { } rest
? head + rest
: null;
}
}
}
Loading
Loading