diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index cdc3c4a0a..ed52626dd 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -665,6 +665,23 @@ taken apart into pieces each closed the same way | `"(2 + x)/((2 + 4*x - 3*x^2)*(1 + 3*x + 2*x^2)^(3/2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative | | `"1/((x^2 + 1)*sqrt(x^2 + x + 1))".ToEntity().Integrate("x")` | `integral(...)` | a logarithm and an arctangent | +### A rational function of `x^n` beside a power of a symbolic binomial is integrated without a search + +**Answers where there were none.** A rational function of `x^n` beside `(c + d x^n)^(k - 1/n)` is a +rational function of `u = x/(c + d x^n)^(1/n)`, and the integrand in `u` was simplified before it +was integrated. With symbols in it the search for a simpler form grew with the power of the +binomial: `1/((a + b x^2)^(5/2) (1 + x^2))` took 5 s on the unreleased master, the `7/2` power ran +past 30 s, and both were declined on 2.5.0. What the search found was the factor the substitution +leaves on both sides of the bar, which is cancelled as it is written now, the other factors left as +they are written ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/((a + b*x^2)^(5/2)*(1 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | powers of `x/sqrt(a + b x^2)`, and an arctangent or a logarithm by the sign of `a - b`, in 0.6 s | +| `"1/((a + b*x^2)^(7/2)*(1 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | the same, in 0.7 s | +| `"1/((a + b*x^2)^(9/2)*(1 + x^2))".ToEntity().Integrate("x")` | `integral(...)` | the same, in 1 s | +| `"x^2/((a + b*x^2)^(7/2)*(c + d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | the same by the sign of `c (a d - b c)`, in 1 s | + ### A power of a multiple of a quadratic's derivative beside a power of the quadratic is a binomial **Improvement, not silent.** `(b d + 2 c d x)^m (a + b x + c x^2)^p`, with a power that is not whole, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 60538086d..a268e555e 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -4207,10 +4207,17 @@ node is Powf(var xAgain, Number.Integer power) && xAgain == x && power.EInteger. var uToTheN = MathS.Pow(u, n); var oneMinusDUToTheN = 1 - d * uToTheN; var binomialInU = c / oneMinusDUToTheN; - var integrand = Functions.SingleQuotient.Combine( + // The factor the substitution leaves on both sides cancelled as it is written, not simplified: + // `(1 - u^4)/((1 + u^4)(1 - u^4))` for Timofeev's above, which the rational integrator does + // not cancel, and the search for a simpler form that did grows with the power of the binomial + // while the symbols are in it: `1/((a + b x^2)^(7/2) (1 + x^2))` ran past 90 s on it. Each + // base is its polynomial in u first, which collects what the substitution writes twice -- + // `1 - (-u^2) - u^2` is 1 under a secant's half angle -- and the factors stay as they are + // written, which the rules after this read better than their product expanded. + var integrand = WithTheFactorsWrittenOnBothSidesCancelled(u, Functions.SingleQuotient.Combine( inXToTheN.Substitute(xToTheN, c * uToTheN / oneMinusDUToTheN) * (k.IsZero ? Number.Integer.One : MathS.Pow(binomialInU, Number.Integer.Create(k))) - / oneMinusDUToTheN).Simplify(); + / oneMinusDUToTheN).InnerSimplified); if (integrand is Providedf(var inner, _)) integrand = inner; if (integrand.ContainsNode(x) || integrand.Nodes.Any(node => node == MathS.NaN)) @@ -4221,6 +4228,66 @@ node is Powf(var xAgain, Number.Integer power) && xAgain == x && power.EInteger. return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// with each base that is a sum written as its polynomial in + /// , and each factor then written alike above and below the bar + /// cancelled, power for power; the factors are kept as factors. + /// + private static Entity WithTheFactorsWrittenOnBothSidesCancelled(Entity.Variable u, Entity quotient) + { + var (top, bottom) = Functions.SingleQuotient.Of(quotient); + Entity AsItsPolynomial(Entity @base) + { + if (@base is not (Sumf or Minusf) || !@base.ContainsNode(u) + || !TreeAnalyzer.TryGetPolynomial(@base, u, out var terms) + || terms.Any(term => term.Key.Sign < 0 || term.Value.ContainsNode(u))) + return @base; + Entity written = Number.Integer.Zero; + foreach (var term in terms.OrderBy(term => term.Key)) + { + var coefficient = term.Value.InnerSimplified; + if (Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + var monomial = term.Key.IsZero ? coefficient : coefficient * MathS.Pow(u, Number.Integer.Create(term.Key)); + written = written == Number.Integer.Zero ? monomial : written + monomial; + } + return written.InnerSimplified; + } + List<(Entity Base, EInteger Power)> Factors(Entity side) + { + var factors = new List<(Entity Base, EInteger Power)>(); + foreach (var factor in Mulf.LinearChildren(side)) + { + var (@base, power) = factor is Powf(var raised, Number.Integer exponent) && exponent.EInteger.Sign > 0 + ? (raised, exponent.EInteger) : (factor, EInteger.One); + @base = AsItsPolynomial(@base); + var at = factors.FindIndex(written => written.Base == @base); + if (at < 0) + factors.Add((@base, power)); + else + factors[at] = (@base, factors[at].Power.Add(power)); + } + return factors; + } + var above = Factors(top); + var below = Factors(bottom); + for (var i = 0; i < above.Count; i++) + { + if (above[i].Base is Number) + continue; + var j = below.FindIndex(written => written.Base == above[i].Base); + if (j < 0) + continue; + var common = EInteger.Min(above[i].Power, below[j].Power); + above[i] = (above[i].Base, above[i].Power.Subtract(common)); + below[j] = (below[j].Base, below[j].Power.Subtract(common)); + } + static Entity Product(List<(Entity Base, EInteger Power)> factors) + => factors.Where(written => !written.Power.IsZero).Aggregate((Entity)Number.Integer.One, + (product, written) => product * (written.Power.Equals(EInteger.One) ? written.Base : MathS.Pow(written.Base, Number.Integer.Create(written.Power)))); + return (Product(above) / Product(below)).InnerSimplified; + } + /// /// Reads as c + d x^n with c and d free of /// . diff --git a/Sources/Tests/UnitTests/Calculus/DividedByTheRootIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/DividedByTheRootIntegralTest.cs new file mode 100644 index 000000000..6c95653a8 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/DividedByTheRootIntegralTest.cs @@ -0,0 +1,69 @@ +// +// 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 rational function of x^n times (c + d x^n)^(k - 1/n), which + /// u = x/(c + d x^n)^(1/n) makes a rational function of u. + /// #718 + /// + /// + /// Checked by differentiating back at points where the integrand is real, with the symbols + /// pinned after the integral is taken. + /// + [Trait("Area", "Calculus")] + public sealed class DividedByTheRootIntegralTest + { + private static void DifferentiatesBack(string integrand, double[] points, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + var derivative = integral.Substitute("C", 0).Differentiate("x"); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + derivative = derivative.Substitute(name, value); + original = original.Substitute(name, value); + } + var compared = 0; + foreach (var at in points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + if (got.IsNaN || want.IsNaN) + continue; + compared++; + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)); + Assert.True(difference / scale < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 3, $"only {compared} points could be compared for {integrand}"); + } + + /// Timofeev's, which is 1/(1 + u^4) under the substitution. + [Fact] + public void TimofeevsQuarticBesideAFourthRoot() + => DifferentiatesBack("1/((1 + x^4)*(2 + x^4)^(1/4))", new[] { -1.3, -0.4, 0.3, 0.9, 1.7 }); + + /// + /// A higher power of the root with symbols in it, beside another quadratic: the rational + /// function of u grows with the power, and the larger ones ran past 90 s. + /// + [Theory] + [InlineData("1/((a + b*x^2)^(5/2)*(1 + x^2))")] + [InlineData("1/((a + b*x^2)^(7/2)*(1 + x^2))")] + [InlineData("1/((a + b*x^2)^(9/2)*(1 + x^2))")] + public void AHalfOddPowerOfASymbolicQuadratic(string integrand) + => DifferentiatesBack(integrand, new[] { -1.3, -0.4, 0.3, 0.9, 1.7 }, ("a", 1.3), ("b", 0.7)); + } +}