From 0173545044518f88a2370a4608b27eef31c954f2 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 6 Oct 2026 09:54:21 +0000 Subject: [PATCH] A quotient in x^2 over a power of x and a biquadratic is split in x^2 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 18 ++++ .../Integration/IndefiniteIntegralSolver.cs | 92 ++++++++++++++----- ...verAPowerOfXAndABiquadraticIntegralTest.cs | 59 ++++++++++++ 3 files changed, 146 insertions(+), 23 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index b45aa2ed8..f6dd61c5c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -336,6 +336,24 @@ quotient of two such linears the sum is `(b - d t^2)^2 + (c t^2 - a)^2`. Rubi's | `"(c + d*tan(x))^(3/2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)` | the same | | `"1/((a + b*tan(x))^(3/2)*(c + d*tan(x))^(3/2))".ToEntity().Integrate("x")` | `integral(...)`; past a minute on the unreleased master | in the root of the quotient of the two, in a second | +### A quotient in `x^2` over a power of `x` and a biquadratic is split in `x^2` + +**Shorter answers, sooner.** `cot(x)^(13/2) (a + b tan(x))^(5/2) (A + B tan(x))` was declined in 2.5.0 +and answered in a million characters after eighteen seconds since; it is answered in two thousand, +in under a second. Under `u = tan(x)` and the root of the quotient `t = sqrt(u/(a + b u))`, Rubi's +half-odd powers of the cotangent beside one of `a + b tan(x)` leave a polynomial in `t^2` over +`t^(2k) Q(t^2)`, `Q` a biquadratic, which was split in `t`: the repeated `t^(2k)` beside the quartic +went to the conjugates of a sum of two squares, at length, or, with the power of a product +distributed, to the Hermite reduction, which declined it after seconds. With a power of `x` on both +sides cancelled, such a quotient is split in `w = x^2` now: the terms in `w^(-j)` are the expansion +of the remainder over `Q` at `w = 0`, and what is left is `(d + e w)/Q`, which the biquadratic's rule +answers by arctangents ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"cot(x)^(13/2)*(a + b*tan(x))^(5/2)*(A + B*tan(x))".ToEntity().Integrate("x")` | `integral(...)`; 1,018,377 characters after 18 seconds on the unreleased master | 2,305 characters, in under a second | +| `"cot(x)^(5/2)/(a + b*tan(x))^(3/2)".ToEntity().Integrate("x")` | `integral(...)`; 4,779 characters after 6 seconds on the unreleased master | 955 characters | + ### A function of `x^n` for a symbolic `n` beside a power of `x` is integrated in a power of `x` **Answers where there were none.** `x^(-1 + 4n)/(a + b x^n + c x^(2n))` and the rest of Rubi's diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 222c4f620..49cf0c67e 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -2422,7 +2422,15 @@ over is Entity.Powf(var @base, var power) ? /// (d + e r_1)/(q (x^2 - r_1)) - (d + e r_2)/(q (x^2 - r_2)), and /// 1/(x^2 - r) is atan(x/sqrt(-r))/sqrt(-r) for any complex r. A /// higher even degree is divided down first, and an odd numerator is left to the - /// substitution u = x^2. + /// substitution u = x^2. And over x^(2k) times the biquadratic, a power of + /// x on both sides cancelled first: in w = x^2 the remainder R over + /// w^k Q is the expansion of R/Q at w = 0 up to w^(k-1), over + /// w^k, plus (d + e w)/Q, so the terms in x^(-2j) go by the power + /// rule and the rest as above. Split in x, the repeated factor x^(2k) beside + /// the quartic went to the Hermite reduction or to the conjugates of a sum of two squares, + /// and 2 a (1 - b t^2)^4/(t^4 (1 - 2 b t^2 + (a^2 + b^2) t^4)), which the root of the + /// quotient makes of Rubi's cot(c + d x)^(5/2)/(a + b tan(c + d x))^(3/2), was + /// declined after seconds in one spelling and answered at length in the other. /// /// /// @@ -2441,14 +2449,30 @@ over is Entity.Powf(var @base, var power) ? { if (!TryReadAsQuotient(expr, out var numerator, out var denominator)) return null; - if (!TreeAnalyzer.TryGetPolynomial(denominator, x, out var below) || below.Count == 0 - || !below.ContainsKey(EInteger.FromInt32(4)) || below.Keys.Any(power => !(power.IsZero || power.Equals(EInteger.FromInt32(2)) || power.Equals(EInteger.FromInt32(4)))) - || below.Values.Any(coefficient => coefficient.ContainsNode(x)) + if (!TreeAnalyzer.TryGetPolynomial(denominator, x, out var belowAsWritten) || belowAsWritten.Count == 0 + || belowAsWritten.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32()) + || belowAsWritten.Values.Any(coefficient => coefficient.ContainsNode(x)) || !denominator.Vars.Any(v => v != x)) return null; - var a = below.TryGetValue(EInteger.Zero, out var a0) ? a0 : Number.Integer.Zero; - var b = below.TryGetValue(EInteger.FromInt32(2), out var b0) ? b0 : Number.Integer.Zero; - var c = below[EInteger.FromInt32(4)]; + if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var aboveAsWritten) || aboveAsWritten.Count == 0 + || aboveAsWritten.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32()) + || aboveAsWritten.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + // A power of x on both sides is cancelled first, and what is left below is x^(2k) times the + // biquadratic: `2 a (1 - b t^2)^4 t/(t^5 ((1 - b t^2)^2 + a^2 t^4))` is how the root of + // a quotient writes Rubi's `cot(c + d x)^(5/2)/(a + b tan(c + d x))^(3/2)`. + var common = EInteger.Min(aboveAsWritten.Keys.Min()!, belowAsWritten.Keys.Min()!).ToInt32Unchecked(); + var below = belowAsWritten.ToDictionary(pair => pair.Key.ToInt32Unchecked() - common, pair => pair.Value); + var above = aboveAsWritten.ToDictionary(pair => pair.Key.ToInt32Unchecked() - common, pair => pair.Value); + var lowest = below.Keys.Min(); + if (lowest % 2 != 0 || !below.ContainsKey(lowest + 4) + || below.Keys.Any(power => power != lowest && power != lowest + 2 && power != lowest + 4) + || above.Keys.Any(power => power % 2 != 0)) + return null; + var k = lowest / 2; + var a = below[lowest]; + var b = below.TryGetValue(lowest + 2, out var b0) ? b0 : Number.Integer.Zero; + var c = below[lowest + 4]; // Zero as a value, and not only as written: the coefficients are read off the // denominator as it is written, and `((c - u^2)/c - 1)(u^2 - c)` is `-u^2 (u^2 - c)/c`, // with `-(c/c - 1) c` for its constant term. Read as a biquadratic, one of its roots @@ -2456,31 +2480,45 @@ over is Entity.Powf(var @base, var power) ? // https://github.com/asc-community/AngouriMath/issues/1665 if (Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(c)) return null; - if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var above) || above.Count == 0 - || above.Keys.Any(power => power.Sign < 0 || !power.IsEven) || above.Values.Any(coefficient => coefficient.ContainsNode(x))) - return null; var discriminant = Functions.PartialFractions.Bare((b * b - 4 * a * c).Simplify()); if (Functions.PartialFractions.IsZeroAsAValue(discriminant)) return null; - // The numerator in w = x^2, divided down by c w^2 + b w + a to a remainder d + e w. - var degree = above.Keys.Max()!.ToInt32Checked() / 2; - var inW = new Entity[degree + 1]; - for (var k = 0; k <= degree; k++) - inW[k] = above.TryGetValue(EInteger.FromInt32(2 * k), out var at) ? at : Number.Integer.Zero; + // The numerator in w = x^2, divided down by w^k (c w^2 + b w + a) to a remainder of a + // degree below k + 2. + var degree = above.Keys.Max() / 2; + var inW = new Entity[System.Math.Max(degree, k + 1) + 1]; + for (var j = 0; j < inW.Length; j++) + inW[j] = above.TryGetValue(2 * j, out var at) ? at : Number.Integer.Zero; Entity polynomialPart = Number.Integer.Zero; - for (var k = degree; k >= 2; k--) + for (var j = degree; j >= k + 2; j--) { - var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k] / c); + var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j] / c); if (lead == Number.Integer.Zero) continue; - polynomialPart += lead * MathS.Pow(x, 2 * (k - 2)); - inW[k - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k - 1] - lead * b); - inW[k - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[k - 2] - lead * a); - } - var d = inW[0]; - var e = degree >= 1 ? inW[1] : Number.Integer.Zero; + polynomialPart += lead * MathS.Pow(x, 2 * (j - k - 2)); + inW[j - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 1] - lead * b); + inW[j - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 2] - lead * a); + } + // The remainder R over w^k Q is sum r_m w^(m - k) over m < k, the expansion of R/Q at + // w = 0, plus (d + e w)/Q: r_m = (R_m - b r_(m-1) - c r_(m-2))/a, and d + e w is what + // (R - Q sum r_m w^m)/w^k leaves, every lower coefficient cancelling by the recurrence. + var pole = new Entity[k]; + for (var m = 0; m < k; m++) + { + var term = inW[m]; + if (m >= 1) term -= b * pole[m - 1]; + if (m >= 2) term -= c * pole[m - 2]; + pole[m] = Functions.PartialFractions.InLowestTermsOverTheSymbols(term / a); + } + var dTerm = inW[k]; + if (k >= 1) dTerm -= b * pole[k - 1]; + if (k >= 2) dTerm -= c * pole[k - 2]; + var eTerm = inW[k + 1]; + if (k >= 1) eTerm -= c * pole[k - 1]; + var d = k == 0 ? inW[0] : Functions.PartialFractions.InLowestTermsOverTheSymbols(dTerm); + var e = k == 0 ? (degree >= 1 ? inW[1] : Number.Integer.Zero) : Functions.PartialFractions.InLowestTermsOverTheSymbols(eTerm); var q = MathS.Sqrt(discriminant); var firstRoot = (-b + q) / (2 * c); @@ -2492,6 +2530,14 @@ over is Entity.Powf(var @base, var power) ? return null; total += whole; } + // r_m x^(2m - 2k), whose exponent is odd and so never minus one. + for (var m = 0; m < k; m++) + { + if (pole[m] == Number.Integer.Zero) + continue; + var exponent = 2 * (m - k) + 1; + total += pole[m] * MathS.Pow(x, exponent) / exponent; + } // `1/(x^2 - r)` is `atan(x/s)/s` with `s = sqrt(-r)` for every complex `r` but zero: // `d/dx atan(x/s)/s` is `1/(s^2 + x^2)` whatever `s` is, and for a positive `r` the // arctangent of an imaginary argument is the hyperbolic one, `-atanh(x/sqrt(r))/sqrt(r)`, diff --git a/Sources/Tests/UnitTests/Calculus/EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest.cs new file mode 100644 index 000000000..7f30498ab --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest.cs @@ -0,0 +1,59 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A polynomial in x^2 over x^(2k) times a biquadratic with a symbol in it, a power + /// of x on both sides cancelled first, split in w = x^2: the terms in + /// w^(-j) are the expansion of the remainder over the biquadratic at w = 0, and + /// what is left is the biquadratic's own. A half-odd power of the cotangent beside one of + /// a + b tan(x) comes to one under u = tan(x) and the root of the quotient + /// t = sqrt(u/(a + b u)); split in t, the answer to the last row ran to a million + /// characters. Rubi's 4.3.2.1 and 4.3.3.1. + /// #718 + /// + /// + /// The cotangent's rows are real where tan(x) and a + b tan(x) are positive, + /// which the points are. The pinned answer is differentiated as it is, since simplifying + /// answers of this shape costs minutes. + /// + [Trait("Area", "Calculus")] + public sealed class EvenQuotientOverAPowerOfXAndABiquadraticIntegralTest + { + [Theory] + [InlineData("1/(x^2*(a + b*x^2 + c*x^4))")] + [InlineData("(1 + x^10)/(x^2*(a + b*x^2 + c*x^4))")] + [InlineData("2*a*(1 - b*x^2)^4*x/(x^5*(1 - 2*b*x^2 + (a^2 + b^2)*x^4))")] + [InlineData("cot(x)^(5/2)/(a + b*tan(x))^(3/2)")] + [InlineData("cot(x)^(9/2)*sqrt(a + b*tan(x))*(A + B*tan(x))")] + [InlineData("cot(x)^(13/2)*(a + b*tan(x))^(5/2)*(A + B*tan(x))")] + public void SplitInTheSquare(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 2.9) + .Substitute("A", 0.4).Substitute("B", 1.1); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { 0.2, 0.5, 0.9, 1.2 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + } +}