From 826524a8489c0136f521742a90ba8a23b502121e Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 21 Oct 2026 03:44:40 +0000 Subject: [PATCH] Split an odd half power of a linear over the second power of x (b + c x) over its two linears (d + e x)^(7/2)/(b x + c x^2)^2 was answered by a rule late in the chain that writes the whole integrand as one rational function in t = sqrt(d + e x) with b e - c d written out, and whose reduction never put its coefficients in lowest terms: denominators like b^35 c^5 d^31 k^35, one coefficient of 155,000 characters copied into four arms, 937,000 characters in all. The rule that keeps x (b + c x) for a root above the bar now takes the second power as well as the third and up: in s = t^2 the integrand is a polynomial part and partial fractions over s - d and s + rho, rho = (b e - c d)/c, whose coefficients come from the Taylor series at each root, each linear's fractions asked as one question, with the atoms d, b e/c, e and P's coefficients named -- naming whole coefficients hid a cancellation. Rubi's 1.2.1.2:435 is 2,568 characters from 936,939, and 1.2.1.3:1386 3,933 from 939,897, in about 10 MB from 10-11 GB. Of 83 rows of 1.2.1.2 and 1.2.1.3 with a half-odd power of a linear over a power of b x + c x^2 or a + b x + c x^2, 20 are shorter and 63 unchanged, character for character, and no verdict is lost. The #1976 entry's two third-power figures follow. Part of #233. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 23 ++- .../Integration/IndefiniteIntegralSolver.cs | 183 ++++++++++++++---- .../Integration/Integration.Definition.cs | 4 +- ...OfAQuadraticWithARootAtZeroIntegralTest.cs | 49 ++++- 4 files changed, 211 insertions(+), 48 deletions(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 0ff0b1fc7..cbaef64ae 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -4211,6 +4211,25 @@ in both quadratics, the cube's coefficients grow past what the reduction takes, | `"(7 + 13*x)/((5 + x + 2*x^2)^3*sqrt(2 + x + 3*x^2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative | | `"sqrt(a + a*sec(x))/(c + d*sec(x))^2".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative in `tan(x/2)` | +### An odd half power of a linear over the second power of `x (b + c x)` is split over its two linears + +**Answers where there were none.** `(d + e x)^(7/2)/(b x + c x^2)^2` and the rest of Rubi's 1.2.1.2 and 1.2.1.3 +with a polynomial times an odd half power of a linear over the second power of `b x + c x^2` were declined. Under +`t = sqrt(d + e x)` the integrand is rational in `s = t^2`, and it is split there: its quotient by the denominator, +and its partial fractions over `s - d` and `s + (b e - c d)/c`, which are `x` and `b + c x`, with the coefficients +of the Taylor series at each root. On the way to this release these were answered as one rational function in +`t`, in 2,900 to 940,000 characters, and the third power and higher of the entry below in 10,000 to 66,000; split, +each is under 7,000, in megabytes where the second power took up to 11 GB +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(d + k*x)^(7/2)/(b*x + c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 2,568 characters | +| `"(A + B*x)*(d + k*x)^(7/2)/(b*x + c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 3,933 characters | +| `"(d + k*x)^(1/2)/(b*x + c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 1,995 characters | +| `"(d + k*x)^(1/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 3,519 characters | +| `"(A + B*x)*(d + k*x)^(9/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 6,723 characters | + ### Cube roots of two linears beside a power of a third, and an odd half power of a linear over a power of `x (b + c x)`, are integrated **Answers where there were none.** `(a + b x)^(1/3) (c + d x)^(2/3)/(e + f x)^2` and its kin, Rubi's 1.1.1.3, and @@ -4230,8 +4249,8 @@ higher. Each declines where its linears share a root, decided as a value | `"(a + b*x)^(1/3)*(c + d*x)^(2/3)/(k + f*x)^2".ToEntity().Integrate("x")` | `integral(...)` | 2,243 characters | | `"1/((a + b*x)^(1/3)*(c + d*x)^(2/3)*(k + f*x)^3)".ToEntity().Integrate("x")` | `integral(...)` | 2,242 characters | | `"(a + b*x)^(4/3)/((c + d*x)^(4/3)*(k + f*x)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,327 characters | -| `"(d + k*x)^(1/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 10,875 characters | -| `"(A + B*x)*(d + k*x)^(1/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 14,763 characters | +| `"(d + k*x)^(1/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 3,519 characters | +| `"(A + B*x)*(d + k*x)^(1/2)/(b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,383 characters | ### Powers of two linears whose exponents sum to a whole number below `-2` are integrated diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index fb9203ee3..66808da3d 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14320,26 +14320,35 @@ static Entity Reduced(Entity constant) /// /// A polynomial times an odd half power of a linear over a whole power of a quadratic with /// a root at zero, P(x) (d + e x)^(j/2)/(b x + c x^2)^m with j odd and - /// positive and m >= 3, under t = sqrt(d + e x) with the quadratic written - /// as its two linears, x (b + c x). - /// - /// - /// - /// With x = (t^2 - d)/e, b + c x = (c t^2 + r)/e for r = b e - c d, and - /// dx = 2 t dt/e, the integrand is - /// 2 e^(2m - 1) P((t^2 - d)/e) t^(j + 1)/((t^2 - d)^m (c t^2 + r)^m), and exactly so: - /// t^j is (d + e x)^(j/2) for the principal root and every complex x. - /// The two factors below the bar have the resultant b e, and a compound coefficient - /// is named by a symbol of its own and written back into the answer, as in - /// . - /// - /// - /// Rubi's (A + B x) (d + e x)^(j/2)/(b x + c x^2)^3 ran past the budget: the rational - /// function in t is answered in a second, and the substitution that reaches it is - /// late in the chain and multiplies the quadratic out. Only the root above the bar and a - /// third power or higher, with a symbol to name: the root below the bar, the second power, - /// and numbers for coefficients are answered as they are, and left to that. Declined where - /// d or r is zero as a value, the linear then sharing a root with the quadratic. + /// positive and m >= 2, under t = sqrt(d + e x) with the quadratic written + /// as its two linears, x (b + c x), and split over them. + /// + /// + /// + /// With s = t^2, x = (s - d)/e, b + c x = c (s + rho)/e for + /// rho = r/c and r = b e - c d, and dx = 2 t dt/e, the integrand is + /// 2 e^(2m - 1)/c^m N(s)/((s - d)^m (s + rho)^m) with + /// N(s) = s^((j + 1)/2) P((s - d)/e), and exactly so: t^j is + /// (d + e x)^(j/2) for the principal root and every complex x. That is its + /// quotient by the monic denominator, a polynomial in s, and its partial fractions + /// over s - d and s + rho, whose coefficients are those of the Taylor series + /// of N over the other power at each root. The roots are b e/c apart, which + /// is not zero. Each linear's fractions are asked as one question in t. Each of + /// d, b e/c, e and the coefficients of P that is compound is + /// named by a symbol of its own and written back into the answer, as in + /// , and rho is + /// written in them, so that a power of x in P cancels. + /// + /// + /// Rubi's (A + B x) (d + e x)^(j/2)/(b x + c x^2)^3 ran past the budget: the + /// substitution that reaches it is late in the chain and multiplies the quadratic out. + /// Asked as one rational function in t, the second power came back in up to 940,000 + /// characters and 11 GB, (d + e x)^(7/2)/(b x + c x^2)^2 among them, and the third + /// in 10,000 to 66,000; split, each is under 7,000 characters, in megabytes. The first + /// power and numbers for coefficients are answered compactly as they are, and left to + /// that, and so is the root below the bar, which puts a power of s below it too. + /// Declined where d or r is zero as a value, the linear then sharing a root + /// with the quadratic. /// https://github.com/asc-community/AngouriMath/issues/233 /// /// @@ -14370,7 +14379,7 @@ static Entity Reduced(Entity constant) } if (root is not var (linear, half) || quadratic is not var (q, n) || half.Sign <= 0 || half.Multiply(ERational.FromInt32(2)) is not { } twice || !twice.IsInteger() || !twice.ToEInteger().CanFitInInt32() - || n < 3 || n > MaximumQuadraticPowerUnderARoot) + || n < 2 || n > MaximumQuadraticPowerUnderARoot) return null; var j = twice.ToEInteger().ToInt32Checked(); if (j > MaximumHalfPowerOfTheLinear @@ -14385,40 +14394,136 @@ static Entity Reduced(Entity constant) || VanishesIdentically(e) || VanishesIdentically(d) || VanishesIdentically(r)) return null; + // P((t^2 - d)/e) e^degree as a polynomial in s = t^2, whose coefficient of s^l is the sum + // over i of p_i e^(degree - i) binomial(i, l) (-d)^(i - l), and with the quadratic's + // coefficients and r there is something to name, or the rules after this answer it. + var coefficients = new Entity[degree + 1]; + for (var l = 0; l <= degree; l++) + { + Entity coefficient = Number.Integer.Zero; + foreach (var (i, p) in read.Select(pair => (pair.Key.ToInt32Checked(), pair.Value)).Where(pair => pair.Item1 >= l)) + coefficient += p * MathS.Pow(e, degree - i) * Binomial(i, l) * MathS.Pow(-d, i - l); + coefficients[l] = VanishesIdentically(coefficient) ? Number.Integer.Zero : Reduced(coefficient); + } + if (coefficients.Concat(new[] { Reduced(-d), Reduced(c), Reduced(r) }).All(value => value is Number || IsAName(value))) + return null; + + // The integrand is 2 e^(2n - 1 - degree)/c^n times N(s)/((s - d)^n (s + rho)^n) in + // s = t^2, with N = s^h P((s - d)/e) e^degree and rho = r/c: its polynomial part, and + // its partial fractions over s - d and s + rho, whose coefficients are those of the + // Taylor series at each root of N over the other power. The atoms are named and rho is + // written in them, so that a power of x in P cancels there. + var h = (j + 1) / 2; var t = Variable.CreateUnique(expr, "t_root"); var named = new List<(Variable Name, Entity Value)>(); Entity Named(Entity value) { - if (value is Number or Variable || value is Mulf(Number, Variable) || value is Powf(Variable, Number)) + if (value is Number || IsAName(value)) return value; var name = Variable.CreateUnique(expr + t + named.Aggregate((Entity)Number.Integer.Zero, (sum, pair) => sum + pair.Name), "k_named"); named.Add((name, value)); return name; } - // P((t^2 - d)/e) e^degree as a polynomial in t^2: the coefficient of t^(2l) is the sum - // over i of p_i e^(degree - i) binomial(i, l) (-d)^(i - l). - Entity above = Number.Integer.Zero; + var atD = Named(Reduced(d)); + var delta = Named(Reduced(b * e / c)); + var rho = delta - atD; + var slope = Named(Reduced(e)); + var ofP = read.Select(pair => (Power: pair.Key.ToInt32Checked(), Coefficient: Named(Reduced(pair.Value)))).ToList(); + var numerator = new Entity[h + degree + 1]; + for (var l = 0; l < h; l++) + numerator[l] = Number.Integer.Zero; for (var l = 0; l <= degree; l++) { Entity coefficient = Number.Integer.Zero; - foreach (var (i, p) in read.Select(pair => (pair.Key.ToInt32Checked(), pair.Value)).Where(pair => pair.Item1 >= l)) - coefficient += p * MathS.Pow(e, degree - i) * Binomial(i, l) * MathS.Pow(-d, i - l); - // t^0 is written as nothing: InnerSimplified makes it 1 provided t is not 0. - if (!VanishesIdentically(coefficient)) - above += l == 0 ? Named(Reduced(coefficient)) : Named(Reduced(coefficient)) * MathS.Pow(t, 2 * l); + foreach (var (i, p) in ofP.Where(pair => pair.Power >= l)) + coefficient += p * PowerOf(slope, degree - i) * Binomial(i, l) * PowerOf(-atD, i - l); + numerator[h + l] = Reduced(coefficient); } + // The shift a question in t is asked over, named apart from what it is in the atoms. + var shift = Named(Reduced(r / c)); var square = MathS.Sqr(t); - var below = MathS.Pow(square + Named(Reduced(-d)), n) * MathS.Pow(Named(Reduced(c)) * square + Named(Reduced(r)), n); - // With nothing to name the coefficients are numbers, and the rules after this answer them. - if (named.Count == 0) - return null; - var inT = (MathS.Pow(t, j + 1) * above / below).InnerSimplified; - if (Integration.ComputeAsAQuestionOfItsOwn(inT, t, integrateByParts) is not { } answer - || answer.Nodes.Any(node => node == MathS.NaN)) - return null; + + // The quotient by the monic (s - d)^n (s + rho)^n, integrated term by term. + Entity inT = Number.Integer.Zero; + if (numerator.Length > 2 * n) + { + var monic = new Entity[2 * n + 1]; + for (var k = 0; k <= 2 * n; k++) + monic[k] = Number.Integer.Zero; + for (var first = 0; first <= n; first++) + for (var second = 0; second <= n; second++) + monic[first + second] += Binomial(n, first) * PowerOf(-atD, n - first) * Binomial(n, second) * PowerOf(rho, n - second); + var remainder = numerator.ToArray(); + for (var top = numerator.Length - 1; top >= 2 * n; top--) + { + var quotient = Reduced(remainder[top]); + if (quotient == Number.Integer.Zero) + continue; + for (var k = 0; k < 2 * n; k++) + remainder[top - 2 * n + k] -= quotient * monic[k]; + var power = 2 * (top - 2 * n) + 1; + inT += quotient * MathS.Pow(t, power) / power; + } + } + + // Each power of s - d and of s + rho over a coefficient named apart, and each sum + // asked as one question: the answer over one quadratic is one piecewise in the sign + // of its constant. + var parts = new List<(Variable Name, Entity Value)>(); + foreach (var (at, gap, linearInS) in new[] { (atD, delta, square - atD), (-rho, -delta, square + shift) }) + { + Entity? fractions = null; + foreach (var (power, value) in TaylorCoefficients(at, gap).Select((value, m) => (n - m, value))) + { + if (VanishesIdentically(value)) + continue; + var name = Variable.CreateUnique(expr + t + named.Concat(parts).Aggregate((Entity)Number.Integer.Zero, (sum, pair) => sum + pair.Name), "k_part"); + parts.Add((name, value)); + var term = name / (power == 1 ? linearInS : MathS.Pow(linearInS, power)); + fractions = fractions is null ? term : fractions + term; + } + if (fractions is null) + continue; + if (Integration.ComputeAsAQuestionOfItsOwn(fractions, t, integrateByParts) is not { } part + || part.Nodes.Any(node => node == MathS.NaN)) + return null; + inT += part; + } + var answer = inT; + foreach (var (name, value) in parts) + answer = answer.Substitute(name, value); foreach (var (name, value) in named) answer = answer.Substitute(name, value); - return Reduced(2 * constant * MathS.Pow(e, 2 * n - 1 - degree)) * answer.Substitute(t, MathS.Sqrt(linear)); + return Reduced(2 * constant * MathS.Pow(e, 2 * n - 1 - degree) / MathS.Pow(c, n)) * answer.Substitute(t, MathS.Sqrt(linear)); + + // The coefficients of u^0 to u^(n - 1) in N(at + u)/(u + gap)^n, in lowest terms. + Entity[] TaylorCoefficients(Entity at, Entity gap) + { + var shifted = new Entity[n]; + for (var q = 0; q < n; q++) + { + Entity sum = Number.Integer.Zero; + for (var l = q; l < numerator.Length; l++) + if (numerator[l] != Number.Integer.Zero) + sum += numerator[l] * Binomial(l, q) * PowerOf(at, l - q); + shifted[q] = sum; + } + var series = new Entity[n]; + for (var m = 0; m < n; m++) + { + Entity sum = Number.Integer.Zero; + for (var q = 0; q <= m; q++) + sum += shifted[q] * Binomial(n + m - q - 1, m - q) * Number.Integer.Create((m - q) % 2 == 0 ? 1 : -1) / MathS.Pow(gap, n + m - q); + series[m] = Reduced(sum); + } + return series; + } + + static bool IsAName(Entity value) + => value is Variable || value is Mulf(Number, Variable) || value is Powf(Variable, Number); + + static Entity PowerOf(Entity @base, int power) + => power == 0 ? Number.Integer.One : power == 1 ? @base : MathS.Pow(@base, power); static Entity Reduced(Entity constant) => constant.Vars.Any() ? Functions.PartialFractions.InLowestTermsOverTheSymbols(constant) : constant.InnerSimplified; diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 1f6e0233e..d5f911f82 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -743,8 +743,8 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // roots with the coefficients named, before the polynomial term takes it a level down to // the substitution there, whose rational function ran past the budget. if ((answer = IndefiniteIntegralSolver.SolveCubeRootsOfTwoLinearsBesideAWholePowerOfAThird(expr, x, integrateByParts)) is { }) return answer; - // And a root of a linear over a third power or higher of `b x + c x^2`, under the root, - // the quadratic kept as its two linears: the substitution there multiplies it out. + // And a root of a linear over a second power or higher of `b x + c x^2`, under the root, + // split over the quadratic's two linears: the substitution there multiplies it out. if ((answer = IndefiniteIntegralSolver.SolveAnOddHalfPowerOfALinearOverAPowerOfAQuadraticWithARootAtZero(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByWritingEachLinearOnce(expr, x, integrateByParts)) is { }) return answer; // The same with a symbol for an exponent, where the product is one power that is not diff --git a/Sources/Tests/UnitTests/Calculus/AnOddHalfPowerOfALinearOverAPowerOfAQuadraticWithARootAtZeroIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AnOddHalfPowerOfALinearOverAPowerOfAQuadraticWithARootAtZeroIntegralTest.cs index 13268e5cc..44233aba5 100644 --- a/Sources/Tests/UnitTests/Calculus/AnOddHalfPowerOfALinearOverAPowerOfAQuadraticWithARootAtZeroIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/AnOddHalfPowerOfALinearOverAPowerOfAQuadraticWithARootAtZeroIntegralTest.cs @@ -12,11 +12,12 @@ namespace AngouriMath.Tests.Calculus { /// - /// A polynomial times an odd half power of a linear over a third power or higher of - /// b x + c x^2, under t = sqrt(d + k x) with the quadratic kept as - /// x (b + c x): the substitution reached the other way multiplies it out, and these - /// ran past the budget. Rubi's 1.2.1.2 and 1.2.1.3. Points on both sides of the zeros, where - /// the root is imaginary. + /// A polynomial times an odd half power of a linear over a second power or higher of + /// b x + c x^2, under t = sqrt(d + k x) and split over the quadratic's two + /// linears x and b + c x: the substitution reached the other way multiplies it + /// out, and these ran past the budget or came back in up to 940,000 characters. Rubi's 1.2.1.2 + /// and 1.2.1.3. Points on both sides of the zeros where the root is imaginary, and in each + /// interval between them where it is real. /// #233 /// [Trait("Area", "Calculus")] @@ -52,5 +53,43 @@ Entity Pinned(Entity e) => e.Substitute("A", 1.1).Substitute("B", 0.6).Substitut } Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); } + + [Theory] + [InlineData("(d + k*x)^(7/2)/(b*x + c*x^2)^2", 3500)] + [InlineData("(d + k*x)^(9/2)/(b*x + c*x^2)^2", 3500)] + [InlineData("(d + k*x)^(1/2)/(b*x + c*x^2)^2", 3500)] + [InlineData("(A + B*x)*(d + k*x)^(7/2)/(b*x + c*x^2)^2", 4500)] + [InlineData("(A + B*x)*(d + k*x)^(5/2)/(b*x + c*x^2)^2", 4500)] + [InlineData("(A + B*x)*(d + k*x)^(9/2)/(b*x + c*x^2)^3", 8000)] + [InlineData("x^2*(d + k*x)^(5/2)/(b*x + c*x^2)^2", 2000)] + public void SplitOverTheTwoLinearsInEachIntervalWhereTheRootIsReal(string integrand, int longest) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < longest, $"{text.Length} characters of answer for {integrand}"); + // The linear's root left of both of the quadratic's, -4.44 against -1.86 and 0, and + // then -0.44 against 0 and 1.86: two points in each interval. + foreach (var (c, d, points) in new[] + { + (0.7, 4.0, new[] { -4.1, -3.3, -1.5, -0.6, 0.4, 2.9 }), + (-0.7, 0.4, new[] { -0.3, -0.1, 0.5, 1.3, 2.4, 3.9 }), + }) + { + Entity Pinned(Entity e) => e.Substitute("A", 1.1).Substitute("B", 0.6).Substitute("b", 1.3) + .Substitute("c", c).Substitute("d", d).Substitute("k", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in points) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.False(want.IsNaN, $"{integrand} has no value at x = {at}"); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, c = {c}, d = {d}, where the integrand is {want}"); + } + } + } } }