Skip to content

Commit 4e302c2

Browse files
A repeated quadratic denominator is integrated (#1155)
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. 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: J_m = (u / Q^(m-1) + 2a(2m - 3) J_(m-1)) / ((m - 1) D) unrolled down to J_1, which 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. The two shapes the reduction cannot speak for get their own arms -- a zero leading coefficient leaves a linear power, a zero discriminant leaves a power of u -- which is what lets a symbolic quadratic be answered as a piecewise instead of carrying a division by zero. A numerator of higher degree is divided by the quadratic first: N = qQ + r takes one power off the denominator and two degrees off the numerator, and ends. An improper fraction is declined rather than divided out here, since the recursion bottoms out at a single power with a numerator of degree two or more -- exactly the case SolveByPartialFractions already opens with a long division. The remainder is computed as N - qQ rather than read from the division, which reports it already divided by the divisor: x^2 over x^2 + 2 answers (1, -2/(x^2 + 2)), and taking that as a polynomial gives a numerator that is not one, so every polynomial-numerator case declined. A numerator already along the denominator's derivative leaves no second integral and the term is dropped rather than multiplied by zero. That is not tidiness: 0 * (2x/(x^2 + 1)) is `0 provided not 1 + x^2 = 0` and not 0, so multiplying would attach a condition to x/(x^2 + 1)^2, which had none and is owed none. 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 metric that ranks candidates selects it. Measured, both arms on one machine, the same way. Where the integrand is rational this is 26 to 165 times faster, because TryStandardIntegrals runs before every search and a shape answered outright never enters the candidate exploration: 1/(x^2 + 1)^2 goes from 2879 ms unevaluated to 111 ms answered, and x^2/(x^2 + 2)^2 from 92,990 ms unevaluated to 565 ms answered. Where the integrand is a transcendental function over a repeated quadratic, none of which has an elementary antiderivative, the same search takes 3.4 to 10.8 times longer to conclude nothing -- e^x/(x^2 + 1)^2 from 7944 ms to 85,683 ms. That is growth in integration by parts, which those integrands reach and this rule never matches; it is filed separately rather than bounded here, since by-parts is machinery every integral goes through. Full suite 9515 passed, 0 failed. 26 new tests, each answer differentiated back and compared with the integrand. Part of #180. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura
1 parent 0e643f3 commit 4e302c2

4 files changed

Lines changed: 390 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -111,6 +111,61 @@ read first.
111111
| | `"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)` |
112112
| | `"sum(k, k, 1, 100000)".ToEntity().Simplify()`, and every concrete range past a hundred terms | `sum(k, k, 1, 100000)` — carried | `5000050000` |
113113
| | `"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)` |
114+
| | `"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 |
115+
| | `"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 |
116+
| | `"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 |
117+
118+
### A repeated quadratic denominator is integrated
119+
120+
`1/(x^2 + 1)^2` had no antiderivative, while `1/(x^2 - 1)^2` and `1/(x^2 + 2x + 1)^2` both did — a
121+
denominator with real roots comes apart into linear factors and never reached the gap. What was
122+
missing was the irreducible case
123+
([#180](https://github.com/asc-community/AngouriMath/issues/180)).
124+
125+
```
126+
"1/(x^2 + 1)^2".Integrate("x") was integral(1 / (1 + x ^ 2) ^ 2, x)
127+
is arctan(x) / 2 + C + 1/2 * x / (x ^ 2 + 1)
128+
"x^2/(x^2 + 2)^2".Integrate("x") was integral(x ^ 2 / (x ^ 2 + 2) ^ 2, x)
129+
is (sqrt(2) * arctan(sqrt(2) * x / 2) * 2 + (-4) * x / (x ^ 2 + 2)) / 8
130+
+ C provided not x ^ 2 + 2 = 0
131+
```
132+
133+
The condition on that second one is the library's own and it is true: the antiderivative has no
134+
value where the denominator vanishes. It is attached to some of these answers and not others, and I
135+
did not pin down what decides which — recorded as observed rather than explained.
136+
137+
Writing `Q` for the quadratic, `u` for its derivative `2ax + b` and `D` for `4ac - b^2`, the
138+
identity `u^2 = 4aQ - D` turns the derivative of `u/Q^(m-1)` into an equation for the integral of
139+
`1/Q^m`, which is unrolled down to the single power the table already answered. It is an identity in
140+
`a`, `b` and `c` rather than a fact about signs, so one line serves both signs of the discriminant;
141+
`D = 0` is the sole exclusion, and the division by `D` is why. A numerator of higher degree is
142+
divided by `Q` first.
143+
144+
**A denominator with real roots now takes this route rather than partial fractions.** The value is
145+
unchanged and the form is shorter — one logarithm where there were two — which is why the
146+
node-count metric that ranks candidates selects it. No test pinned the old spelling.
147+
148+
**The cost, measured rather than estimated.** Where the integrand is rational this is between 26 and
149+
165 times faster, because `TryStandardIntegrals` runs before every search and a shape answered
150+
outright never enters the candidate exploration. Where the integrand is a transcendental function
151+
over a repeated quadratic — `sin(x)/(x^2 + 1)^2`, `e^x/(x^2 + 1)^2`, `ln(x)/(x^2 + 1)^2`, none of
152+
which has an elementary antiderivative — the same search now takes **3.4 to 10.8 times longer to
153+
conclude nothing**, because the sub-integral it used to fail on immediately now succeeds and
154+
integration by parts carries on from a larger expression. Both arms measured on one machine, master
155+
and branch, the same way:
156+
157+
| | master | now |
158+
|---|---|---|
159+
| `1/(x^2 + 1)^2` | 2879 ms, unevaluated | **111 ms, answered** |
160+
| `x^2/(x^2 + 2)^2` | 92 990 ms, unevaluated | **565 ms, answered** |
161+
| `1/(x^4 + 1)^2` | 3811 ms | 3870 ms — the base is not quadratic, and nothing changes |
162+
| `ln(x)/(x^2 + 1)^2` | 15 034 ms, unevaluated | 51 072 ms, unevaluated |
163+
| `sin(x)/(x^2 + 1)^2` | 26 492 ms, unevaluated | 114 851 ms, unevaluated |
164+
| `e^x/(x^2 + 1)^2` | 7944 ms, unevaluated | 85 683 ms, unevaluated |
165+
166+
That last group is a growth in integration by parts rather than in this rule, which those integrands
167+
never match; it is filed separately rather than fixed here, since bounding by-parts is a change to
168+
machinery every integral goes through.
114169

115170
### A nested radical comes apart
116171

‎Sources/.editorconfig‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -189,3 +189,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
189189

190190
[Tests/UnitTests/PatternsTest/DenestRadicalTest.cs]
191191
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
192+
193+
[Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs]
194+
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

‎Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs‎

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,18 @@ _ when TryReadOverSquareTimesRoot(expr, x, out var overFactor, out var overRadic
223223
&& a.Evaled is Entity.Number.Complex { IsZero: false }
224224
=> IntegrateLinearOverQuadratic(p, q, a, b, c, denominator, x),
225225

226+
// ∫ N(x)/(ax^2 + bx + c)^n dx for a whole n of two or more and a proper fraction.
227+
// After the single-power arms above, which already read a denominator like
228+
// (x + 1)^2 -- expanded, it is a quadratic, so the perfect-square arm there answers
229+
// it and nothing about that changes.
230+
Entity.Divf(var numerator, Entity.Powf(var repeated, Entity.Number.Integer repetitions)) when
231+
repetitions.EInteger.CanFitInInt32()
232+
&& repetitions.EInteger.ToInt32Checked() is var n and >= 2
233+
&& TreeAnalyzer.TryGetPolyQuadratic(repeated, x, out var qa, out var qb, out var qc)
234+
&& (!numerator.ContainsNode(x) || qa.Evaled is Entity.Number.Complex { IsZero: false })
235+
&& IntegrateRationalOverPowerOfQuadratic(numerator, qa, qb, qc, repeated, n, x) is { } repeatedAnswer
236+
=> repeatedAnswer,
237+
226238
// ∫ (px + q)/(bx + c) dx, the same rewrite one degree down: the quotient is the
227239
// constant p/b plus a remainder over the divisor. x/(x + 1) had no antiderivative.
228240
Entity.Divf(var numerator, var denominator) when
@@ -529,5 +541,170 @@ private static Entity IntegrateRationalQuadratic(Entity numerator, Entity a, Ent
529541
new Entity.Providedf(lnCase, discriminant < 0)
530542
]).InnerSimplified;
531543
}
544+
545+
/// <summary>
546+
/// ∫ k / (a x^2 + b x + c)^n dx for a whole n of two or more, by the reduction that
547+
/// takes one power off the denominator at a time.
548+
/// </summary>
549+
/// <remarks>
550+
/// <para>
551+
/// Write <c>Q</c> for the quadratic, <c>u</c> for its derivative <c>2ax + b</c> and
552+
/// <c>D</c> for <c>4ac - b^2</c>. Then <c>u^2</c> is <c>4aQ - D</c>, so differentiating
553+
/// <c>u / Q^(m-1)</c> gives <c>2a(3 - 2m)/Q^(m-1) + (m-1)D/Q^m</c>, and reading that as
554+
/// an equation for the integral of <c>1/Q^m</c> leaves
555+
/// <code>
556+
/// J_m = (u / Q^(m-1) + 2a(2m - 3) J_(m-1)) / ((m - 1) D)
557+
/// </code>
558+
/// which is unrolled down to <c>J_1</c>, the case the table above already answers. It is
559+
/// an identity in a, b and c rather than a fact about signs, so one line serves a
560+
/// positive and a negative discriminant alike; <c>D = 0</c> is the only exclusion, and
561+
/// the division by <c>D</c> is why.
562+
/// </para>
563+
/// <para>
564+
/// The two shapes the reduction cannot speak for get their own arms, and each is
565+
/// elementary: with <c>a = 0</c> the denominator is a linear power, and with
566+
/// <c>D = 0</c> the quadratic is <c>u^2/(4a)</c> so the whole integrand is a power of
567+
/// <c>u</c>. Both are already answered for a numeric quadratic — by the power rule and
568+
/// by partial fractions — and the arms are here for a symbolic one, where the
569+
/// discriminant's sign is not known and dropping an arm is not available.
570+
/// </para>
571+
/// <para>
572+
/// What was missing was the irreducible case. <c>1/(x^2 - 1)^2</c> had an antiderivative
573+
/// and <c>1/(x^2 + 2x + 1)^2</c> had one, because a denominator with real roots comes
574+
/// apart into linear factors and never reaches here; <c>1/(x^2 + 1)^2</c> had none.
575+
/// https://github.com/asc-community/AngouriMath/issues/180
576+
/// </para>
577+
/// </remarks>
578+
private static Entity IntegrateOverPowerOfQuadratic(
579+
Entity numerator, Entity a, Entity b, Entity c, int power, Entity.Variable x)
580+
{
581+
var quadratic = a * x * x + b * x + c;
582+
var derivative = 2 * a * x + b;
583+
var discriminant = 4 * a * c - b * b;
584+
585+
// a = 0: k/(bx + c)^n, an ordinary power -- and the constant k/c^n when b is zero
586+
// too, where writing the power would divide by it. Same reason as the rule above.
587+
var linearCase = TreeAnalyzer.IsZero(b)
588+
? numerator * x / MathS.Pow(c, power)
589+
: numerator * MathS.Pow(b * x + c, 1 - power) / (b * (1 - power));
590+
591+
// D = 0: the quadratic is u^2/(4a), so the integrand is k(4a)^n / u^(2n).
592+
var perfectSquareCase =
593+
numerator * MathS.Pow(4 * a, power) * MathS.Pow(derivative, 1 - 2 * power)
594+
/ (2 * a * (1 - 2 * power));
595+
596+
var sqrtDiscriminant = MathS.Sqrt(discriminant);
597+
var sqrtNegDiscriminant = MathS.Sqrt(-discriminant);
598+
599+
return MathS.Piecewise([
600+
new Entity.Providedf(linearCase, a.EqualTo(0)),
601+
new Entity.Providedf(perfectSquareCase, discriminant.EqualTo(0)),
602+
new Entity.Providedf(
603+
Reduce(2 * MathS.Arctan(derivative / sqrtDiscriminant) / sqrtDiscriminant),
604+
discriminant > 0),
605+
new Entity.Providedf(
606+
Reduce(AntiderivativeLog((derivative - sqrtNegDiscriminant) / (derivative + sqrtNegDiscriminant))
607+
/ sqrtNegDiscriminant),
608+
discriminant < 0)
609+
]).InnerSimplified;
610+
611+
Entity Reduce(Entity firstPower)
612+
{
613+
var accumulated = firstPower;
614+
for (var m = 2; m <= power; m++)
615+
accumulated =
616+
(derivative / MathS.Pow(quadratic, m - 1) + 2 * a * (2 * m - 3) * accumulated)
617+
/ ((m - 1) * discriminant);
618+
return numerator * accumulated;
619+
}
620+
}
621+
622+
/// <summary>
623+
/// ∫ (px + q)/(a x^2 + b x + c)^n dx, by the same rewrite of the numerator the single
624+
/// power uses: <c>px + q = (p/2a)(2ax + b) + (q - pb/2a)</c>. The first part is the
625+
/// denominator's own derivative over a power of it, which integrates as a power; the
626+
/// second is the constant-numerator case above.
627+
/// </summary>
628+
/// <remarks>
629+
/// A numerator that is already a multiple of that derivative leaves nothing of the
630+
/// second part, and it is dropped rather than multiplied by zero. Multiplying is not
631+
/// harmless: the reduction is a quotient by the quadratic, so <c>0 * (2x/(x^2 + 1))</c>
632+
/// is <c>0 provided not 1 + x^2 = 0</c> and not <c>0</c> — correctly, since the factor
633+
/// has no value there. Carrying that condition would attach it to the answer for
634+
/// <c>x/(x^2 + 1)^2</c>, which had none before and is owed none.
635+
/// </remarks>
636+
private static Entity IntegrateLinearOverPowerOfQuadratic(
637+
Entity p, Entity q, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x)
638+
{
639+
var alongTheDerivative = p / (2 * a) * MathS.Pow(quadratic, 1 - power) / (1 - power);
640+
var constantPart = q - p * b / (2 * a);
641+
return TreeAnalyzer.IsZero(constantPart)
642+
? alongTheDerivative
643+
: alongTheDerivative + IntegrateOverPowerOfQuadratic(constantPart, a, b, c, power, x);
644+
}
645+
646+
/// <summary>
647+
/// ∫ N(x)/(a x^2 + b x + c)^n dx for a numerator of any degree, or <see langword="null"/>
648+
/// where the fraction is improper and this is not the rule for it.
649+
/// </summary>
650+
/// <remarks>
651+
/// <para>
652+
/// Dividing the numerator by the quadratic writes <c>N = qQ + r</c> with <c>r</c> linear
653+
/// at worst, so <c>N/Q^n</c> is <c>q/Q^(n-1) + r/Q^n</c> — one power off the denominator
654+
/// and two degrees off the numerator at each step, which ends. The remainder is answered
655+
/// by the two helpers above and the quotient goes round again.
656+
/// </para>
657+
/// <para>
658+
/// <b>Improper fractions are declined rather than divided out here.</b> The recursion
659+
/// bottoms out at a single power, and a numerator still of degree two or more there is
660+
/// exactly the improper case, which
661+
/// <see cref="IndefiniteIntegralSolver.SolveByPartialFractions"/> already opens with
662+
/// <see cref="TreeAnalyzer.PolynomialLongDivision"/>. Declining hands it there instead of
663+
/// answering it twice, in two shapes, from two places.
664+
/// </para>
665+
/// <para>
666+
/// This is not only reach. <c>TryStandardIntegrals</c> runs before every search, so a
667+
/// shape answered here never reaches the candidate exploration that re-integrates
668+
/// rewritten forms — which is where this family's cost has always been. Both
669+
/// <c>x^2/(x^2 + 2)^2</c> and <c>(x^2 + 1)/(x^2 + 2)^2</c> were tens of seconds spent to
670+
/// return the integral unevaluated.
671+
/// </para>
672+
/// </remarks>
673+
private static Entity? IntegrateRationalOverPowerOfQuadratic(
674+
Entity numerator, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x)
675+
{
676+
if (!numerator.ContainsNode(x))
677+
return power == 1
678+
? IntegrateRationalQuadratic(numerator, a, b, c, x)
679+
: IntegrateOverPowerOfQuadratic(numerator, a, b, c, power, x);
680+
681+
if (TreeAnalyzer.TryGetPolyLinear(numerator, x, out var p, out var q))
682+
return power == 1
683+
? IntegrateLinearOverQuadratic(p, q, a, b, c, quadratic, x)
684+
: IntegrateLinearOverPowerOfQuadratic(p, q, a, b, c, quadratic, power, x);
685+
686+
if (power < 2)
687+
return null;
688+
689+
var division = TreeAnalyzer.PolynomialLongDivision(numerator, quadratic);
690+
if (division is null)
691+
return null;
692+
693+
// The quotient is what is wanted; the remainder that comes back with it is already
694+
// divided by the divisor -- x^2 over x^2 + 2 answers (1, -2/(x^2 + 2)) -- so taking
695+
// it as a polynomial and dividing it again is a numerator that is not one, and the
696+
// next division declines. Subtracting gives the polynomial remainder whatever shape
697+
// the division chose to report.
698+
// Reduced before either is used: the division reports its quotient built up term by
699+
// term, so x^2 over x^2 + 2 comes back as `0 + 1` rather than `1`, and the arms below
700+
// ask whether a numerator is constant.
701+
var quotient = division.Value.Divided.InnerSimplified;
702+
var remainder = (numerator - quotient * quadratic).InnerSimplified;
703+
704+
return IntegrateRationalOverPowerOfQuadratic(remainder, a, b, c, quadratic, power, x) is { } head
705+
&& IntegrateRationalOverPowerOfQuadratic(quotient, a, b, c, quadratic, power - 1, x) is { } rest
706+
? head + rest
707+
: null;
708+
}
532709
}
533710
}

0 commit comments

Comments
 (0)