diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index cb3188b36..eabde4cfc 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -294,6 +294,20 @@ angle, each of which the integrator answers | `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))".ToEntity().Integrate("x")` | `integral(...)` | 508 characters | | `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))^2".ToEntity().Integrate("x")` | `integral(...)` | 1,470 characters | +### A constant common to a half-angle denominator is divided out + +**Answers where there were none.** `x/(k + k sin(x))^2` was declined, and so were `x/(d c + c d sin(x))^2` and Rubi's +5.1.4 `(a + b asin(c x))/((d + c d x)^(5/2) sqrt(f - c f x))`. After the half-angle substitution, such a denominator +is a constant times a polynomial with number coefficients, `k^2 (t + 1)^4`. The constant is now divided out before the +partial fractions +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x/(k+k*sin(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 285 characters | +| `"x/(d*c+c*d*sin(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 329 characters | +| `"(a+b*asin(c*x))/((d+c*d*x)^(5/2)*sqrt(f-c*f*x))".ToEntity().Integrate("x")` | `integral(...)` | 927 characters | + ### A piecewise is defined where the case it takes is **Domain conditions that were false.** The domain condition of a piecewise asked every case's predicate diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a7c60b5d9..903dd96a1 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -32678,6 +32678,33 @@ bool IsAPolynomialInBoth(Entity polynomial) ?.Pipe(answer => underneath ? answer / constant : answer * constant); } + /// + /// A rational function of as a polynomial over a polynomial, both divided by + /// the denominator's leading coefficient with each coefficient in lowest terms over the symbols, + /// so that a constant common to the denominator's coefficients cancels; null where either side is + /// not a polynomial in with coefficients free of it. + /// + private static Entity? WithTheLeadingCoefficientDividedOut(Entity rational, Entity.Variable t) + { + var (above, below) = Functions.SingleQuotient.Of(rational); + if (!TreeAnalyzer.TryGetPolynomial(above.Expand(), t, out var top) || !TreeAnalyzer.TryGetPolynomial(below.Expand(), t, out var bottom) + || bottom.Count == 0 || top.Keys.Concat(bottom.Keys).Any(key => key.Sign < 0 || !key.CanFitInInt32()) + || top.Values.Concat(bottom.Values).Any(coefficient => coefficient.ContainsNode(t))) + return null; + var lead = bottom[bottom.Keys.Max()!]; + // Only where the symbols in the denominator are one factor common to all of it: there the + // rules below formed `k^2/k^2`, and anywhere else this respelling sent fifteen of Rubi's + // 4.x problems the half-angle answered in a second to a timeout. + if (!lead.Vars.Any() || bottom.Values.Any(coefficient => Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficient / lead) is not Number)) + return null; + Entity Rebuilt(Dictionary monomials) + => monomials.Aggregate((Entity)Number.Integer.Zero, (sum, pair) + => sum + Functions.PartialFractions.InLowestTermsOverTheSymbols(pair.Value / lead) + // t^0 written as 1: as a power it brings `provided not t = 0` with it. + * (pair.Key.IsZero ? Number.Integer.One : MathS.Pow(t, Number.Integer.Create(pair.Key)))); + return (Rebuilt(top) / Rebuilt(bottom)).InnerSimplified; + } + /// /// A rational function of sin(x) and cos(x), turned into a rational function /// of one variable by the half-angle substitution t = tan(x/2) and handed to the @@ -32869,6 +32896,14 @@ static bool NoRemainder(Entity rest) if (divided) integrand = (numerator / denominator).InnerSimplified; + // Over the denominator's leading coefficient, every coefficient in lowest terms over the + // symbols: `1/(k + k sin(x))^2` is `2 (1 + t^2)/(k^2 (t + 1)^4)`, and with `k^2` left in the + // coefficients the rules below formed `k^2/k^2`, which the simplification does not cancel, + // so the quadratic's discriminant was not decidably zero and `x/(k + k sin(x))^2` was + // answered `NaN` through a branch of the piecewise that can never hold. + if (integrand.Vars.Any(symbol => symbol != t) && WithTheLeadingCoefficientDividedOut(integrand, t) is { } normalized) + integrand = normalized; + return Integration.ComputeIndefiniteIntegral(integrand, t, integrateByParts) is { } result ? result.Substitute(t, MathS.Tan(argument / 2)) : null; diff --git a/Sources/Tests/UnitTests/Calculus/AConstantCommonToAHalfAngleDenominatorIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AConstantCommonToAHalfAngleDenominatorIntegralTest.cs new file mode 100644 index 000000000..1a64ef34b --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AConstantCommonToAHalfAngleDenominatorIntegralTest.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 constant common to the whole denominator after the half-angle substitution, divided out: + /// x/(k + k sin(x))^2 was NaN, and so was Rubi's 5.1.4 + /// (a + b asin(c x))/((d + c d x)^(5/2) sqrt(f - c f x)), which comes to it. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class AConstantCommonToAHalfAngleDenominatorIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + + var pinned = integral.Substitute("C", 0); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + pinned = pinned.Substitute(name, value); + original = original.Substitute(name, value); + } + var derivative = pinned.Differentiate("x"); + foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 }) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + 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-8, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + [Theory] + [InlineData("x/(k+k*sin(x))^2")] + [InlineData("1/(d*q+q*d*sin(x))^2")] + [InlineData("(a+b*asin(q*x))/((d+q*d*x)^(5/2)*sqrt(f-q*f*x))")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + // |q x| below 1 at every sampled x, where the inverse sine and the roots are real. + DifferentiatesBack(integrand, ("k", 1.3), ("a", 0.7), ("b", -0.4), ("d", 1.1), ("f", 1.5), ("q", 0.3)); + DifferentiatesBack(integrand, ("k", -0.7), ("a", -2), ("b", 1.5), ("d", 0.6), ("f", 2), ("q", -0.2)); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/AHalfAngleSubstitutionWithACommonConstantIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AHalfAngleSubstitutionWithACommonConstantIntegralTest.cs new file mode 100644 index 000000000..2065a0790 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AHalfAngleSubstitutionWithACommonConstantIntegralTest.cs @@ -0,0 +1,63 @@ +// +// 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 +{ + /// + /// The half-angle substitution with a constant common to the denominator's coefficients: it is + /// divided out in lowest terms over the symbols, which the simplification does not do for + /// k^2/k^2. Kept, it left a quadratic whose discriminant was not decidably zero, and + /// x/(k + k sin(x))^2 -- and Rubi's 5.1.4 (a + b asin(c x))/((d + c d x)^(5/2) sqrt(f - c f x)), + /// which comes to it -- were answered NaN. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class AHalfAngleSubstitutionWithACommonConstantIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + + var pinned = integral.Substitute("C", 0); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + pinned = pinned.Substitute(name, value); + original = original.Substitute(name, value); + } + var derivative = pinned.Differentiate("x"); + foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 }) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + 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-8, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + // d with the sign that keeps the argument real at every sampled x, and |d| below 2/2.3^2 + // so that the inverse sine and cosine stay inside [-1, 1]; off it the exact checks hold too. + [Theory] + [InlineData("x/(k+k*sin(x))^2")] + [InlineData("1/(d*q+q*d*sin(x))^2")] + [InlineData("x/(k-k*cos(x))^2")] + [InlineData("(a+b*asin(q*x))/((d+q*d*x)^(5/2)*sqrt(f-q*f*x))")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + Assert.DoesNotContain("NaN", integrand.ToEntity().Integrate("x").Stringize()); + DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25), ("q", 0.3), ("d", 0.75), ("f", 1.4), ("k", 1.5)); + DifferentiatesBack(integrand, ("a", -0.75), ("b", 0.4), ("q", -0.4), ("d", -2), ("f", 0.5), ("k", -0.3)); + } + } +}