Skip to content

Commit 1fce7d8

Browse files
Reduce a rational function over a power of a symbolic quadratic a power at a time (#1891)
* Lower the power of a symbolic quadratic in a rational function one at a time A rational function over a whole power of a symbolic quadratic, beside at most one linear, is taken apart at the linear's root -- through Q - Q(p) = (x - p)(c (x + p) + b), so that every piece stays a short polynomial over a power of L and of Q -- and the power of the quadratic is lowered by the recurrence for the integral of 1/Q^k, down to the first, which is written in closed form. The discriminant, Q(p) and p are named while it runs. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Leave a quadratic that factors, and a question below the top, to the partial fractions b^2 - 4ac a square is a quadratic over two linears, which the partial fractions take apart cheaply: (d + e x)^10/((d + e x)(a e + c d x))^4 went from 2 seconds past the budget here. And below a substitution the routes this pre-empts answered shorter. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Leave a binomial in x^2 and a quadratic over the imaginary unit to the rules before Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the answers a power of a symbolic quadratic reduced gives Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --------- Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent 37b4798 commit 1fce7d8

4 files changed

Lines changed: 282 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -559,6 +559,23 @@ of the quadratic, parts leaves the quadratic one power lower, which is answered
559559
| `"(b + 2*c*x)/((a + b*x + c*x^2)^2*sqrt(d + k*x))".ToEntity().Integrate("x")` | `integral(...)` | 1,075 characters |
560560
| `"(b + 2*c*x)*(d + k*x)^(3/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,054 characters |
561561

562+
### A rational function over a power of a symbolic quadratic is reduced a power at a time
563+
564+
**Answers where there were none.** `1/((d + k x) (a + b x + c x^2)^4)` ran past a minute, the partial fractions taking the
565+
symbolic quadratic apart over its roots and its powers with them. The linear is taken off at its root now, and the power
566+
of the quadratic is lowered one at a time by the recurrence for the integral of `1/Q^k`, down to the first power, which
567+
is written in closed form. A quadratic that factors over linears is left to the partial fractions
568+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
569+
570+
| Input | Was (2.5.0) | Now |
571+
|---|---|---|
572+
| `"x^8/(a + b*x + c*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | 8,511 characters |
573+
| `"1/((d + k*x)*(a + b*x + c*x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 9,091 characters |
574+
| `"1/((d + k*x)^2*(a + b*x + c*x^2)^5)".ToEntity().Integrate("x")` | `integral(...)` | 25,767 characters |
575+
| `"1/((d + k*x)^2*(a + c*x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 9,599 characters |
576+
| `"(f + g*x)/((d + k*x)^2*(a + b*x + c*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 21,943 characters |
577+
| `"(A + B*x + C*x^2)/((d + k*x)^3*(a + c*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 25,201 characters |
578+
562579
### An odd power of a quadratic's derivative over a power of the quadratic is integrated in the quadratic
563580

564581
**Shorter answers, and answers where there were none.** `(b + 2c x)^9/(a + b x + c x^2)^3` was answered through the

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

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19342,6 +19342,208 @@ Entity TakeOut(Entity node)
1934219342
return inTermsOfU.Substitute(u, Number.Integer.One / x);
1934319343
}
1934419344

19345+
/// <summary>
19346+
/// A rational function over a whole power of a symbolic quadratic, beside at most one
19347+
/// linear to a whole power, with the power of the quadratic lowered one at a time rather
19348+
/// than the quadratic taken apart: the linear comes off first, as
19349+
/// <c>N/(L^j Q^n) = N(p)/(Q(p)^n L^j) + M/(L^(j - 1) Q^n)</c> at its root <c>p</c>, and
19350+
/// what is left over <c>Q^k</c> is <c>S/Q^(k - 1) + (alpha Q' + beta)/Q^k</c>, whose second
19351+
/// part is <c>-alpha/((k - 1) Q^(k - 1))</c> and the recurrence
19352+
/// <c>integral 1/Q^k = (2 c x + b)/((k - 1) D Q^(k - 1)) + 2 c (2 k - 3)/((k - 1) D) integral 1/Q^(k - 1)</c>
19353+
/// with <c>D = 4 a c - b^2</c>, down to the first power, which is the integrator's.
19354+
/// </summary>
19355+
/// <remarks>
19356+
/// Rubi's 1.2.1.2 <c>1/((d + e x) (a + b x + c x^2)^4)</c> and <c>x^8/(a + b x + c x^2)^4</c>
19357+
/// ran past a minute: the partial fractions take a symbolic quadratic apart over its
19358+
/// roots, and its powers with it. <c>D</c> and <c>Q(p)</c> are named while the reduction
19359+
/// runs and written back after, so that each stays one atom.
19360+
/// https://github.com/asc-community/AngouriMath/issues/718
19361+
/// </remarks>
19362+
internal static Entity? SolveARationalFunctionOverAPowerOfASymbolicQuadratic(Entity expr, Entity.Variable x, bool integrateByParts)
19363+
{
19364+
// At the question asked only: below a substitution the routes it pre-empts answered
19365+
// `(a + b x + c x^2)^3/((1 - d x)^(3/2) (1 + d x)^(3/2))` in 2,917 characters, and
19366+
// this in 67,746.
19367+
if (!Integration.AnsweringTheQuestionAsked)
19368+
return null;
19369+
var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
19370+
if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var top) || top.Keys.Any(k => k.Sign < 0 || !k.CanFitInInt32()))
19371+
return null;
19372+
Entity? quadratic = null;
19373+
var power = 0;
19374+
Entity? linear = null;
19375+
var linearPower = 0;
19376+
Entity constant = Number.Integer.One;
19377+
foreach (var factor in Mulf.LinearChildren(denominator))
19378+
{
19379+
if (!factor.ContainsNode(x))
19380+
{
19381+
constant = constant * factor;
19382+
continue;
19383+
}
19384+
var (@base, whole) = factor is Powf(var b0, Number.Integer { EInteger: var w }) && w.Sign > 0 && w.CanFitInInt32()
19385+
? (b0, w.ToInt32Unchecked()) : (factor, 1);
19386+
if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) && !VanishesIdentically(slope))
19387+
{
19388+
if (linear is not null && linear != @base)
19389+
return null;
19390+
linear = @base;
19391+
linearPower += whole;
19392+
continue;
19393+
}
19394+
if (quadratic is null && whole >= 2 && TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var leading, out _, out _) && !VanishesIdentically(leading))
19395+
{
19396+
quadratic = @base;
19397+
power = whole;
19398+
continue;
19399+
}
19400+
return null;
19401+
}
19402+
if (quadratic is null || !quadratic.Vars.Any(symbol => symbol != x)
19403+
|| !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var a1, out var a0))
19404+
return null;
19405+
var discriminant = (4 * a2 * a0 - a1 * a1).InnerSimplified;
19406+
if (VanishesIdentically(discriminant))
19407+
return null;
19408+
// Real coefficients only, whose cases the table's 1/Q reads; and with no linear, not a
19409+
// binomial in x^2, which the binomial rules answer: `x^22/(a + b x^2)^10` was theirs in
19410+
// a second and ten levels of this ran past the budget.
19411+
if (Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic) || linear is null && VanishesIdentically(a1))
19412+
return null;
19413+
// A quadratic that factors over linears is the partial fractions' to take apart, and
19414+
// they do it cheaply: `(d + e x)^10/((d + e x) (a e + c d x))^4`, written out, went from
19415+
// two seconds past the budget here.
19416+
// It does where b^2 - 4 a c is the square of a polynomial in the symbols, which is
19417+
// what is asked: the factoring does not find `(d + e x) (a e + c d x)` written out.
19418+
var opposite = Functions.PartialFractions.Bare((a1 * a1 - 4 * a2 * a0).Expand().InnerSimplified);
19419+
var symbols = opposite.Vars.Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray();
19420+
var indices = new Dictionary<Variable, int>();
19421+
for (var i = 0; i < symbols.Length; i++)
19422+
indices[symbols[i]] = i;
19423+
if (symbols.Length <= MultivariatePolynomial.MaxVariables
19424+
&& MultivariatePolynomial.TryParse(opposite, indices) is { } parsed && parsed.TrySquareRoot() is not null)
19425+
return null;
19426+
var named = Variable.CreateUnique(expr, "D_q");
19427+
var atTheRoot = Variable.CreateUnique(expr + named, "Q_p");
19428+
var theRoot = Variable.CreateUnique(expr + named + atTheRoot, "x_p");
19429+
Entity d = named;
19430+
19431+
static Entity[] Trim(Entity[] poly)
19432+
{
19433+
var length = poly.Length;
19434+
while (length > 0 && poly[length - 1] == Number.Integer.Zero)
19435+
length--;
19436+
return poly.Take(length).ToArray();
19437+
}
19438+
static Entity[] Plus(Entity[] left, Entity[] right)
19439+
{
19440+
var sum = new Entity[System.Math.Max(left.Length, right.Length)];
19441+
for (var i = 0; i < sum.Length; i++)
19442+
sum[i] = ((i < left.Length ? left[i] : Number.Integer.Zero) + (i < right.Length ? right[i] : Number.Integer.Zero)).InnerSimplified;
19443+
return sum;
19444+
}
19445+
static Entity[] Scaled(Entity[] poly, Entity by) => poly.Select(c => (c * by).InnerSimplified).ToArray();
19446+
// Each coefficient in lowest terms over the symbols, so that the next step works on it
19447+
// as written small.
19448+
static Entity[] Reduced(Entity[] poly) => poly.Select(Functions.PartialFractions.InLowestTermsOverTheSymbols).ToArray();
19449+
19450+
var degree = top.Keys.Count == 0 ? 0 : top.Keys.Max()!.ToInt32Unchecked();
19451+
var start = new Entity[degree + 1];
19452+
for (var i = 0; i <= degree; i++)
19453+
start[i] = top.TryGetValue(EInteger.FromInt32(i), out var coefficient) ? coefficient : Number.Integer.Zero;
19454+
19455+
// The integrand as a sum of polynomials over L^j Q^k, one per (j, k), merged as they
19456+
// arrive so that each stays of small degree.
19457+
var terms = new Dictionary<(int Linear, int Quadratic), Entity[]>();
19458+
void Add(int j, int k, Entity[] poly)
19459+
=> terms[(j, k)] = terms.TryGetValue((j, k), out var already) ? Plus(already, poly) : poly;
19460+
Add(linearPower, power, start);
19461+
Entity answer = Number.Integer.Zero;
19462+
Entity m = Number.Integer.One, p = theRoot;
19463+
if (linear is not null && !TreeAnalyzer.TryGetPolyLinear(linear, x, out m!, out _))
19464+
return null;
19465+
// Q - Q(p) = (x - p) (c (x + p) + b).
19466+
var cofactor = new Entity[] { (a2 * p + a1).InnerSimplified, a2 };
19467+
19468+
while (terms.Count > 0)
19469+
{
19470+
var (j, k) = terms.Keys.OrderByDescending(key => key.Linear).ThenByDescending(key => key.Quadratic).First();
19471+
var poly = Reduced(terms[(j, k)]);
19472+
terms.Remove((j, k));
19473+
if (Trim(poly).Length == 0)
19474+
continue;
19475+
if (j > 0)
19476+
{
19477+
// N/(L^j Q^k) = N(p)/(Q(p)^k m^j (x - p)^j) + N1/(m L^(j-1) Q^k)
19478+
// - N(p) (c (x + p) + b) sum_i Q(p)^(k-1-i)/(Q(p)^k m L^(j-1) Q^(k-i)),
19479+
// with N1 = (N - N(p))/(x - p).
19480+
Entity atP = Number.Integer.Zero;
19481+
for (var i = poly.Length - 1; i >= 0; i--)
19482+
atP = (atP * p + poly[i]).InnerSimplified;
19483+
var c = Functions.PartialFractions.InLowestTermsOverTheSymbols(atP / MathS.Pow(atTheRoot, Number.Integer.Create(k)) / MathS.Pow(m, Number.Integer.Create(j)));
19484+
answer = answer + (j == 1 ? c * MathS.Ln(linear!) : c * MathS.Pow(linear!, Number.Integer.Create(1 - j)) * MathS.Pow(m, Number.Integer.Create(j - 1)) / (1 - j));
19485+
var n1 = new Entity[System.Math.Max(poly.Length - 1, 0)];
19486+
Entity carry = Number.Integer.Zero;
19487+
for (var i = poly.Length - 1; i >= 1; i--)
19488+
{
19489+
carry = (poly[i] + carry * p).InnerSimplified;
19490+
n1[i - 1] = carry;
19491+
}
19492+
Add(j - 1, k, Scaled(n1, 1 / m));
19493+
for (var i = 0; i < k; i++)
19494+
Add(j - 1, k - i, Scaled(cofactor, (-atP * MathS.Pow(atTheRoot, Number.Integer.Create(-1 - i)) / m).InnerSimplified));
19495+
continue;
19496+
}
19497+
if (k == 0)
19498+
{
19499+
for (var i = 0; i < poly.Length; i++)
19500+
answer = answer + poly[i] * MathS.Pow(x, Number.Integer.Create(i + 1)) / (i + 1);
19501+
continue;
19502+
}
19503+
// P/Q^k: P = S Q + t1 x + t0, and t1 x + t0 = alpha Q' + beta.
19504+
var remainder = poly.ToArray();
19505+
var quotient = Enumerable.Repeat((Entity)Number.Integer.Zero, System.Math.Max(remainder.Length - 2, 0)).ToArray();
19506+
for (var i = remainder.Length - 1; i >= 2; i--)
19507+
{
19508+
var factor = (remainder[i] / a2).InnerSimplified;
19509+
quotient[i - 2] = factor;
19510+
remainder[i] = Number.Integer.Zero;
19511+
remainder[i - 1] = (remainder[i - 1] - factor * a1).InnerSimplified;
19512+
remainder[i - 2] = (remainder[i - 2] - factor * a0).InnerSimplified;
19513+
}
19514+
var t1 = remainder.Length > 1 ? remainder[1] : Number.Integer.Zero;
19515+
var t0 = remainder.Length > 0 ? remainder[0] : Number.Integer.Zero;
19516+
var alpha = Functions.PartialFractions.InLowestTermsOverTheSymbols(t1 / (2 * a2));
19517+
var beta = Functions.PartialFractions.InLowestTermsOverTheSymbols(t0 - alpha * a1);
19518+
Add(0, k - 1, quotient);
19519+
if (k == 1)
19520+
{
19521+
// Over the first power in closed form: alpha ln(Q), and beta times the table's
19522+
// integral of 1/Q. Asked of the integrator instead, a linear over Q with these
19523+
// coefficients took forty seconds.
19524+
if (IntegralPatterns.TryStandardIntegrals(1 / quadratic, x) is not { } overQ)
19525+
return null;
19526+
answer = answer + alpha * MathS.Ln(quadratic) + beta * overQ;
19527+
continue;
19528+
}
19529+
var lower = MathS.Pow(quadratic, Number.Integer.Create(k - 1));
19530+
answer = answer - alpha / ((k - 1) * lower) + beta * (2 * a2 * x + a1) / ((k - 1) * d * lower);
19531+
Add(0, k - 1, new[] { (beta * 2 * a2 * (2 * k - 3) / ((k - 1) * d)).InnerSimplified });
19532+
}
19533+
Entity withoutTheLinear = Number.Integer.Zero;
19534+
Entity scale = Number.Integer.One;
19535+
answer = ((withoutTheLinear + scale * answer) / constant).Substitute(named, discriminant);
19536+
if (linear is not null && TreeAnalyzer.TryGetPolyLinear(linear, x, out var mm, out var nn))
19537+
{
19538+
var root = (-nn / mm).InnerSimplified;
19539+
// Q(p) as one quotient in lowest terms, (c d^2 - b d e + a e^2)/e^2, rather than
19540+
// Q written at -d/e wherever it stood.
19541+
answer = answer.Substitute(atTheRoot, Functions.PartialFractions.InLowestTermsOverTheSymbols(a2 * root * root + a1 * root + a0)).Substitute(theRoot, root);
19542+
}
19543+
answer = answer.InnerSimplified;
19544+
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
19545+
}
19546+
1934519547
/// <summary>
1934619548
/// A power of a linear beside the derivative of a quadratic, over a whole power of that
1934719549
/// quadratic, by parts once: <c>Q' L^m/Q^n</c> with <c>u = L^m</c> and

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

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -669,6 +669,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
669669
// A power of a linear beside the derivative of a quadratic over a whole power of it, by
670670
// parts once, before anything takes the symbolic quadratic apart.
671671
if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearBesideTheDerivativeOfAQuadraticOverItsPower(expr, x, integrateByParts)) is { }) return answer;
672+
// A rational function over a whole power of a symbolic quadratic, the power lowered one
673+
// at a time, before the partial fractions take the quadratic apart over its roots.
674+
if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverAPowerOfASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer;
672675
// A root or a power of a product beside a power of one of its factors, the product
673676
// split: `(a + b x^2)^(5/2)/sqrt(a^2 - b^2 x^4)` is `K (a + b x^2)^2/sqrt(a - b x^2)`.
674677
// Early, since what reads the product whole is slower and longer: Rubi's 1.2.1.2:2218,
Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
//
2+
// Copyright (c) 2019-2026 Angouri.
3+
// AngouriMath is licensed under MIT.
4+
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
5+
// Website: https://am.angouri.org.
6+
//
7+
8+
using System;
9+
using AngouriMath.Extensions;
10+
using Xunit;
11+
12+
namespace AngouriMath.Tests.Calculus
13+
{
14+
/// <summary>
15+
/// A rational function over a whole power of a symbolic quadratic, beside one linear: the
16+
/// linear taken off at its root and the power of the quadratic lowered one at a time. Rubi's
17+
/// 1.2.1.2 <c>1/((d + e x) (a + b x + c x^2)^4)</c>, written here with <c>k</c> for <c>e</c>,
18+
/// which is Euler's number.
19+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
20+
/// </summary>
21+
[Trait("Area", "Calculus")]
22+
public sealed class ARationalFunctionOverAPowerOfASymbolicQuadraticIntegralTest
23+
{
24+
private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
25+
{
26+
var integral = integrand.ToEntity().Integrate("x");
27+
Assert.DoesNotContain("integral(", integral.Stringize());
28+
29+
// Pinned before it is differentiated: the answer is a piecewise on the sign of the
30+
// discriminant, read where it is a number.
31+
var pinned = integral.Substitute("C", 0);
32+
Entity original = integrand.ToEntity();
33+
foreach (var (name, value) in pins)
34+
{
35+
pinned = pinned.Substitute(name, value);
36+
original = original.Substitute(name, value);
37+
}
38+
var derivative = pinned.Differentiate("x");
39+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
40+
{
41+
var got = derivative.Substitute("x", at).EvalNumerical();
42+
var want = original.Substitute("x", at).EvalNumerical();
43+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
44+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
45+
Assert.True(difference / scale < 1e-8,
46+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
47+
}
48+
}
49+
50+
[Theory]
51+
[InlineData("x^8/(a+b*x+c*x^2)^4")]
52+
[InlineData("1/((d+k*x)*(a+b*x+c*x^2)^4)")]
53+
[InlineData("1/((d+k*x)^2*(a+c*x^2)^4)")]
54+
public void IsAnsweredAndDifferentiatesBack(string integrand)
55+
{
56+
DifferentiatesBack(integrand, ("a", 2), ("b", 1.5), ("c", 1.25), ("d", 1.25), ("k", 2));
57+
DifferentiatesBack(integrand, ("a", -1), ("b", 2), ("c", 2), ("d", -0.5), ("k", 1.5));
58+
}
59+
}
60+
}

0 commit comments

Comments
 (0)