diff --git a/Sources/.editorconfig b/Sources/.editorconfig index beb9e9c40..7856516ed 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -366,3 +366,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Calculus/LogarithmSubstitutionTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Calculus/HomogeneousTrigonometricTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a08b9750b..e1d795d24 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -1882,6 +1882,163 @@ node is Logf(var @base, var argument) && !@base.ContainsNode(x) && @base != Math : null; } + /// + /// A quotient of two homogeneous polynomials in sin(u) and cos(u), + /// integrated by t = tan(u) — which turns it into a rational function of t + /// whenever the two degrees differ by an even number. + /// + /// + /// + /// 1/(cos(x) + sin(x))^6 took 231 seconds to be declined and + /// 1/(b^2 cos(x)^2 + a^2 sin(x)^2) was one of the Rubi sample's three timeouts. + /// Both are this shape, and so is sin(x)/(cos(x)^3 + sin(x)^3). + /// + /// + /// Why the degrees decide it. Writing t = tan(u) and c = cos(u), a + /// homogeneous polynomial of degree n in the pair is c^n times a polynomial + /// in t alone. So a quotient of degrees n over d is + /// c^(n-d) N(t)/D(t), and with c^2 = 1/(1 + t^2) and + /// dx = dt/(a(1 + t^2)): + /// + /// + /// f dx = [N(t)/D(t)] (1 + t^2)^((d - n)/2 - 1) dt / a + /// + /// + /// which is rational exactly when d - n is even. Odd is a genuine boundary rather + /// than a first cut: there the substitution leaves a square root of 1 + t^2 behind, + /// which is a different problem and not a rational one. + /// + /// + /// This is the third of Bioche's rules, and it is the one the other two do not cover: + /// answers an integrand that is a function of + /// tan alone, and 1/(b^2 cos^2 + a^2 sin^2) is not — it is a function + /// of tan times sec^2, which is what the degree difference of two is saying. + /// + /// + /// Where the answer holds. t = tan(u) is a bijection on each interval + /// between the poles of the tangent, so the answer is an antiderivative on each of them — + /// the standing caveat on this substitution, shared with the half-angle one, and not + /// something this writes as a condition. + /// + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByHomogeneousTrigonometricSubstitution(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!TryReadAsQuotient(expr, out var numerator, out var denominator)) + return null; + if (!TryReadAHomogeneousTrigonometricPolynomial(numerator, x, out var above, out var aboveDegree, out var argument) + || !TryReadAHomogeneousTrigonometricPolynomial(denominator, x, out var below, out var belowDegree, out var otherArgument)) + return null; + // A side free of the variable is degree zero and names no argument -- `1/(cos + sin)^2` + // is the common case and was declined for it. Only one side may be silent; if both + // are, there is no integrand here. + argument ??= otherArgument; + otherArgument ??= argument; + if (argument is null || argument != otherArgument) + return null; + var difference = belowDegree - aboveDegree; + if (difference % 2 != 0) + return null; // an odd difference leaves a root of 1 + t^2, which is not rational + if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) || rate.Evaled == 0) + return null; + + var t = Variable.CreateUnique(expr, "t_homog"); + var inT = Functions.SingleQuotient.Combine( + above.Substitute(HomogeneousTangent, t).Substitute(HomogeneousCosine, 1) + / below.Substitute(HomogeneousTangent, t).Substitute(HomogeneousCosine, 1) + * MathS.Pow(1 + MathS.Sqr(t), Number.Integer.Create(difference / 2 - 1)) + / rate).Simplify(); + if (inT is Providedf(var inner, _)) + inT = inner; + if (inT.ContainsNode(x)) + return null; + + return Integration.ComputeIndefiniteIntegral(inT, t, integrateByParts) is { } result + ? result.Substitute(t, MathS.Tan(argument)) + : null; + } + + /// The two placeholders a homogeneous polynomial is rewritten over. + [ConstantField] private static readonly Entity.Variable HomogeneousTangent = + Entity.Variable.CreateVariableOrConstant("__homogeneous_tangent"); + + /// See . + [ConstantField] private static readonly Entity.Variable HomogeneousCosine = + Entity.Variable.CreateVariableOrConstant("__homogeneous_cosine"); + + /// + /// Reads as a homogeneous polynomial in sin(u) and + /// cos(u) over one common argument, rewritten with sin(u) as + /// tangent * cosine and cos(u) as cosine so that the caller can put + /// cosine to one and read off the polynomial in the tangent. + /// + /// + /// Homogeneous means every term has the same total degree, and that is checked rather + /// than assumed: 1 + cos(u) has terms of degree zero and one and is declined, which + /// is right — it is not this substitution's, and the half-angle one answers it. + /// + private static bool TryReadAHomogeneousTrigonometricPolynomial( + Entity expr, Entity.Variable x, out Entity rewritten, out int degree, out Entity? argument) + { + rewritten = expr; + degree = 0; + argument = null; + + var found = (Entity?)null; + var theDegree = (int?)null; + foreach (var term in Entity.Sumf.LinearChildren(expr.Expand())) + { + var thisDegree = 0; + foreach (var factor in Entity.Mulf.LinearChildren(term)) + { + var (piece, power) = factor is Powf(var @base, Number.Integer whole) + && whole.EInteger.CanFitInInt32() + ? (@base, whole.EInteger.ToInt32Checked()) + : (factor, 1); + switch (piece) + { + case Sinf(var a): + thisDegree += power; + if (found is not null && found != a) return false; + found = a; + break; + case Cosf(var a): + thisDegree += power; + if (found is not null && found != a) return false; + found = a; + break; + default: + if (piece.ContainsNode(x)) + return false; // anything else in the term is not this shape + break; + } + } + if (theDegree is null) + theDegree = thisDegree; + else if (theDegree != thisDegree) + return false; + } + if (theDegree is null) + return false; + if (found is null) + { + // No sine or cosine at all: a constant, which is homogeneous of degree zero and + // needs no rewriting. It names no argument, and the caller takes the other side's. + degree = 0; + return !expr.ContainsNode(x); + } + + degree = theDegree.Value; + argument = found; + rewritten = expr.Expand().Replace(node => node switch + { + Sinf(var a) when a == found => HomogeneousTangent * HomogeneousCosine, + Cosf(var a) when a == found => HomogeneousCosine, + _ => node + }); + return true; + } + /// /// An integrand that is a function of tan(x) and of nothing else, integrated by /// the substitution u = tan(x), under which dx is du/(1 + u^2). diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 134f9c066..768875ea2 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -411,6 +411,12 @@ private static Entity Normalized(Entity expr) => // is why the general substitution does not find it and why it pays: the integrator // answers an exponential times almost anything. if ((answer = IndefiniteIntegralSolver.SolveByLogarithmSubstitution(expr, x, integrateByParts)) is { }) return answer; + // And the third of Bioche's rules: a quotient of homogeneous polynomials in sine and + // cosine, which the tangent turns into a rational function whenever the two degrees + // differ by an even number. After the tangent substitution above, which answers an + // integrand that is a function of the tangent *alone* and gives a shorter answer for + // it; this one is for the rest. + if ((answer = IndefiniteIntegralSolver.SolveByHomogeneousTrigonometricSubstitution(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByPartialFractions(expr, x, integrateByParts)) is { }) return answer; // The tangent substitution again, and this time for a *rational* function: a repeated // irreducible quadratic beside a negative power of the variable is `sin^p cos^q` with diff --git a/Sources/Tests/UnitTests/Calculus/HomogeneousTrigonometricTest.cs b/Sources/Tests/UnitTests/Calculus/HomogeneousTrigonometricTest.cs new file mode 100644 index 000000000..804c6e57f --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/HomogeneousTrigonometricTest.cs @@ -0,0 +1,150 @@ +// +// 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 quotient of two homogeneous polynomials in sin(u) and cos(u), under + /// t = tan(u) — the third of Bioche's rules. + /// #718 + /// + /// + /// + /// 1/(cos(x) + sin(x))^6 took 239 seconds to be declined and is answered in + /// four; 1/(b^2 cos(x)^2 + a^2 sin(x)^2) was one of the Rubi sample's three timeouts. + /// All three of those timeouts were this shape, and there are now none. + /// + /// + /// Writing t = tan(u) and c = cos(u), a homogeneous polynomial of degree + /// n is c^n times a polynomial in t, so a quotient of degrees n + /// over d becomes N(t)/D(t) times (1 + t^2)^((d-n)/2 - 1) — rational + /// exactly when d - n is even. + /// + /// + /// Checked by differentiating back and comparing at points. The sample points sit inside one + /// interval between the poles of the tangent, which is where the substitution is a bijection + /// and the answer is an antiderivative. + /// + /// + [Trait("Area", "Calculus")] + public sealed class HomogeneousTrigonometricTest + { + private static void DifferentiatesBack(string integrand, double[] points) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("NaN", integral.Stringize()); + Assert.DoesNotContain("integral(", integral.Stringize()); + + var derivative = integral.Substitute("C", 0).Differentiate("x"); + var original = integrand.ToEntity(); + 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)); + 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} of {points.Length} points were comparable for {integrand}, " + + "so this asserts almost nothing"); + } + + /// + /// A power of a sum of a sine and a cosine, which is the case that took 239 seconds to be + /// declined. + /// + [Theory] + [InlineData("1/(cos(x) + sin(x))^6", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/(cos(x) + sin(x))^2", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/(cos(x) + sin(x))^4", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("cos(x)^2/(cos(x) + sin(x))^4", new[] { 0.3, 0.8, 1.2, -0.5 })] + public void APowerOfASumOfTheTwo(string integrand, double[] points) + => DifferentiatesBack(integrand, points); + + /// + /// A quadratic form, which is the shape the Rubi sample carries twice and both times as a + /// timeout. The numerator is a constant there, so it names no argument of its own — the + /// read has to take the denominator's, and declining for that was the first thing this + /// rule got wrong. + /// + [Theory] + [InlineData("1/(4*cos(x)^2 + 9*sin(x)^2)", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/(3*cos(x)^2 - sin(x)^2)", new[] { 0.3, 0.8, -0.5, -0.2 })] + [InlineData("1/(sin(x)^2 + sin(x)*cos(x))", new[] { 0.3, 0.8, 1.2, 2.0 })] + [InlineData("sin(x)^2/cos(x)^4", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("sin(x)/(cos(x)^3 + sin(x)^3)", new[] { 0.3, 0.8, 1.2, -0.5 })] + public void AQuadraticOrCubicForm(string integrand, double[] points) + => DifferentiatesBack(integrand, points); + + /// + /// A linear argument, which divides the whole answer by the rate. + /// + [Theory] + [InlineData("1/(cos(2*x) + sin(2*x))^2", new[] { 0.2, 0.5, 0.9, -0.3 })] + [InlineData("1/(4*cos(3*x)^2 + sin(3*x)^2)", new[] { 0.1, 0.3, 0.4, -0.2 })] + public void ALinearArgument(string integrand, double[] points) + => DifferentiatesBack(integrand, points); + + /// + /// What the read must refuse: a polynomial that is not homogeneous, and a degree + /// difference that is odd — there the substitution leaves a square root of 1 + t^2 + /// behind, which is a different problem. + /// + [Theory] + [InlineData("1/(1 + cos(x))", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/(1 + sin(x))", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/(2 + cos(x))", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("1/cos(x)", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("sin(x)/(1 + cos(x)^2)", new[] { 0.3, 0.8, 1.2, -0.5 })] + public void NotHomogeneousOrAnOddDifference(string integrand, double[] points) + { + // Each of these is answered by some other rule, and the point is that this one does + // not take them and get them wrong. + DifferentiatesBack(integrand, points); + } + + /// + /// Two different arguments, which is not one substitution's — declined rather than + /// guessed at. + /// + [Theory] + [InlineData("1/(cos(x) + sin(2*x))")] + [InlineData("sin(x)/(cos(2*x)^2 + sin(x)^2)")] + public void TwoDifferentArgumentsAreDeclined(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("NaN", integral.Stringize()); + } + + /// + /// The neighbours, each answered by a rule in front of this one and each of which must + /// keep the shorter answer that rule gives. + /// + [Theory] + [InlineData("tan(x)^2", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("tan(x)^3*sec(x)^4", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("sin(x)^3*cos(x)^2", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("sec(x)^4", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("csc(x)^6", new[] { 0.3, 0.8, 1.2, 2.0 })] + [InlineData("sin(x)", new[] { 0.3, 0.8, 1.2, -0.5 })] + [InlineData("x*sin(x)", new[] { 0.3, 0.8, 1.2, -0.5 })] + public void TheNeighboursAreUntouched(string integrand, double[] points) + => DifferentiatesBack(integrand, points); + } +}