A quotient of homogeneous polynomials in sine and cosine, and the Rubi sample has no timeouts left

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

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -1882,6 +1882,163 @@ node is Logf(var @base, var argument) && !@base.ContainsNode(x) && @base != Math
                 : null;
         }
 
+        /// <summary>
+        /// A quotient of two <b>homogeneous</b> polynomials in <c>sin(u)</c> and <c>cos(u)</c>,
+        /// integrated by <c>t = tan(u)</c> — which turns it into a rational function of <c>t</c>
+        /// whenever the two degrees differ by an even number.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <c>1/(cos(x) + sin(x))^6</c> took <b>231 seconds</b> to be declined and
+        /// <c>1/(b^2 cos(x)^2 + a^2 sin(x)^2)</c> was one of the Rubi sample's three timeouts.
+        /// Both are this shape, and so is <c>sin(x)/(cos(x)^3 + sin(x)^3)</c>.
+        /// </para>
+        /// <para>
+        /// <b>Why the degrees decide it.</b> Writing <c>t = tan(u)</c> and <c>c = cos(u)</c>, a
+        /// homogeneous polynomial of degree <c>n</c> in the pair is <c>c^n</c> times a polynomial
+        /// in <c>t</c> alone. So a quotient of degrees <c>n</c> over <c>d</c> is
+        /// <c>c^(n-d) N(t)/D(t)</c>, and with <c>c^2 = 1/(1 + t^2)</c> and
+        /// <c>dx = dt/(a(1 + t^2))</c>:
+        /// </para>
+        /// <code>
+        ///     f dx = [N(t)/D(t)] (1 + t^2)^((d - n)/2 - 1) dt / a
+        /// </code>
+        /// <para>
+        /// which is rational exactly when <c>d - n</c> is even. Odd is a genuine boundary rather
+        /// than a first cut: there the substitution leaves a square root of <c>1 + t^2</c> behind,
+        /// which is a different problem and not a rational one.
+        /// </para>
+        /// <para>
+        /// This is the third of Bioche's rules, and it is the one the other two do not cover:
+        /// <see cref="SolveByTangentSubstitution"/> answers an integrand that is a function of
+        /// <c>tan</c> <i>alone</i>, and <c>1/(b^2 cos^2 + a^2 sin^2)</c> is not — it is a function
+        /// of <c>tan</c> times <c>sec^2</c>, which is what the degree difference of two is saying.
+        /// </para>
+        /// <para>
+        /// <b>Where the answer holds.</b> <c>t = tan(u)</c> 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.
+        /// </para>
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        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;
+        }
+
+        /// <summary>The two placeholders a homogeneous polynomial is rewritten over.</summary>
+        [ConstantField] private static readonly Entity.Variable HomogeneousTangent =
+            Entity.Variable.CreateVariableOrConstant("__homogeneous_tangent");
+
+        /// <summary>See <see cref="HomogeneousTangent"/>.</summary>
+        [ConstantField] private static readonly Entity.Variable HomogeneousCosine =
+            Entity.Variable.CreateVariableOrConstant("__homogeneous_cosine");
+
+        /// <summary>
+        /// Reads <paramref name="expr"/> as a homogeneous polynomial in <c>sin(u)</c> and
+        /// <c>cos(u)</c> over one common argument, rewritten with <c>sin(u)</c> as
+        /// <c>tangent * cosine</c> and <c>cos(u)</c> as <c>cosine</c> so that the caller can put
+        /// <c>cosine</c> to one and read off the polynomial in the tangent.
+        /// </summary>
+        /// <remarks>
+        /// Homogeneous means every term has the same total degree, and that is checked rather
+        /// than assumed: <c>1 + cos(u)</c> 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.
+        /// </remarks>
+        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;
+        }
+
         /// <summary>
         /// An integrand that is a function of <c>tan(x)</c> and of nothing else, integrated by
         /// the substitution <c>u = tan(x)</c>, under which <c>dx</c> is <c>du/(1 + u^2)</c>.

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

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
+{
+    /// <summary>
+    /// A quotient of two <b>homogeneous</b> polynomials in <c>sin(u)</c> and <c>cos(u)</c>, under
+    /// <c>t = tan(u)</c> — the third of Bioche's rules.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    /// <remarks>
+    /// <para>
+    /// <c>1/(cos(x) + sin(x))^6</c> took <b>239 seconds</b> to be declined and is answered in
+    /// four; <c>1/(b^2 cos(x)^2 + a^2 sin(x)^2)</c> was one of the Rubi sample's three timeouts.
+    /// All three of those timeouts were this shape, and there are now none.
+    /// </para>
+    /// <para>
+    /// Writing <c>t = tan(u)</c> and <c>c = cos(u)</c>, a homogeneous polynomial of degree
+    /// <c>n</c> is <c>c^n</c> times a polynomial in <c>t</c>, so a quotient of degrees <c>n</c>
+    /// over <c>d</c> becomes <c>N(t)/D(t)</c> times <c>(1 + t^2)^((d-n)/2 - 1)</c> — rational
+    /// exactly when <c>d - n</c> is even.
+    /// </para>
+    /// <para>
+    /// 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.
+    /// </para>
+    /// </remarks>
+    [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");
+        }
+
+        /// <summary>
+        /// A power of a sum of a sine and a cosine, which is the case that took 239 seconds to be
+        /// declined.
+        /// </summary>
+        [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);
+
+        /// <summary>
+        /// 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.
+        /// </summary>
+        [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);
+
+        /// <summary>
+        /// A linear argument, which divides the whole answer by the rate.
+        /// </summary>
+        [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);
+
+        /// <summary>
+        /// 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 <c>1 + t^2</c>
+        /// behind, which is a different problem.
+        /// </summary>
+        [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);
+        }
+
+        /// <summary>
+        /// Two different arguments, which is not one substitution's — declined rather than
+        /// guessed at.
+        /// </summary>
+        [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());
+        }
+
+        /// <summary>
+        /// The neighbours, each answered by a rule in front of this one and each of which must
+        /// keep the shorter answer that rule gives.
+        /// </summary>
+        [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);
+    }
+}