Split an even quotient over a power of a biquadratic beside another factor in x^2

BREAKING-CHANGES.md

@@ -101,6 +101,25 @@ with the slope in another place, `F^(a + b/(c + d x)^2)/(c + d x)^2` in 15 ms wh
 | `"(k*g+f*g*x)^(p-1)/ln(d*(k+f*x)^p)".ToEntity().Integrate("x")` | `integral(...)` | 150 characters, `Ei(ln(d (k + f x)^p))` |
 | `"F^(a+b/(c+d*x)^2)/(c+d*x)^2".ToEntity().Integrate("x")` | `integral(...)` | `-F^a sqrt(pi) erf(sqrt(-b ln F)/(c + d x))/(2 d sqrt(-b ln F))` |
 
+### An even quotient over a power of a biquadratic beside another factor is split in `x^2`
+
+**Answers where there were none.** `1/((x^2 - 1)(c + d x^2 + f x^4)^2)` ran past the budget, and so did Rubi's 6.1.7
+`csch(c + d x)/(a - b sinh(c + d x)^4)^2`, which is that shape under `w = cosh(c + d x)`, and 1.2.2.3's
+`1/((d + e x^2)(a + b x^2 + c x^4)^2)` took eleven seconds and 10.8 GB. Split in x, the series at the
+roots of `x^2 - 1` came to coefficients of a hundred thousand nodes. An even rational function with a symbol in it, over a
+power of a biquadratic of three terms beside another factor in `x^2`, is now split in `w = x^2`, and each term is asked in
+x. While the split is solved, each factor is made monic and its compound coefficients are named. Every other factor
+stands once, and not where the numerator vanishes at a linear's root: under the tangent `sinh(u)^4/(a - b sinh(u)^4)^2`
+shares a power of `1 - w` with its numerator, and is left to the rules that answered it. Checked by its derivative
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"1/((x^2-1)*(c+d*x^2+f*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,164 characters |
+| `"x^2/((x^2+1)*(c+d*x^2+f*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,062 characters, where before this change master gave 13,714 in 3.8 s |
+| `"csch(c+d*x)/(a-b*sinh(c+d*x)^4)^2".ToEntity().Integrate("x")` | `integral(...)` | 2,511 characters |
+| `"1/((d+k*x^2)*(a+b*x^2+c*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 9,968 characters, in 0.9 s, where before this change master gave 18,154 in 11 s and 10.8 GB |
+
 ### A polynomial over a symbolic quartic that is two quadratics in a hidden variable is split over them
 
 **Answers where there were none.** `1/(a + b x + c x^2 + b x^3 + a x^4)` was declined, and so were Rubi's 1.3.1

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

@@ -4631,6 +4631,117 @@ over is Entity.Powf(var @base, var power) ?
             return total / (constant * MathS.Pow(l.First, k));
         }
 
+        /// <summary>
+        /// An even rational function with a symbol in it, over a power past the first of a biquadratic of
+        /// three terms beside another factor in <c>x^2</c>, split over its written factors in <c>w = x^2</c> and each
+        /// term asked in x. <c>1/((x^2 - 1)(c + d x^2 + f x^4)^2)</c> is
+        /// <c>A/(w - 1) + (B + C w)/Q(w) + (D + E w)/Q(w)^2</c>, and every term is a shape the rules
+        /// answer at once; split in x, the series at the roots of <c>x^2 - 1</c> ran to coefficients
+        /// of a hundred thousand nodes, with powers like <c>c^40 f^58</c>, and past the budget. That
+        /// is Rubi's 6.1.7 <c>csch(u)/(a - b sinh(u)^4)^2</c> under <c>w = cosh(u)</c>. Checked by
+        /// its derivative.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </summary>
+        internal static Entity? SolveAnEvenQuotientOverAPowerOfABiquadraticBesideAnotherFactorInTheSquare(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || !(numerator + denominator).Vars.Any(v => v != x))
+                return null;
+            var w = Variable.CreateUnique(expr, "w");
+            // Each coefficient that is more than a symbol or a number a fresh symbol while the split is
+            // solved and the terms integrated, and put back after: in `a - b (x^2 - 1)^2` they are
+            // `a - b`, `2 b` and `-b`, and the split took 3 s and 1.4 GB where with three symbols it takes 0.1.
+            var named = new List<(Variable Name, Entity Value)>();
+            Entity everything = expr + w;
+            Entity Named(Entity coefficient)
+            {
+                if (coefficient is Variable or Number)
+                    return coefficient;
+                foreach (var (already, value) in named)
+                    if (value == coefficient)
+                        return already;
+                var name = Variable.CreateUnique(everything, "k_q");
+                everything += name;
+                named.Add((name, coefficient));
+                return name;
+            }
+            // A polynomial in x with only even powers, written in w, its degree in w, how many terms it
+            // has and its leading coefficient; monic, where asked, with that coefficient out in front and
+            // the others over it in lowest terms, so that the names come back as `(b - a)/b`.
+            (Entity InW, int Degree, int Terms, Entity Leading)? InTheSquare(Entity polynomial, bool monic)
+            {
+                if (!TreeAnalyzer.TryGetPolynomial(polynomial, x, out var read) || read.Count == 0)
+                    return null;
+                foreach (var pair in read)
+                    if (pair.Key.Sign < 0 || !pair.Key.IsEven || !pair.Key.CanFitInInt32() || pair.Value.ContainsNode(x))
+                        return null;
+                var top = read.Keys.Max()!;
+                var leading = monic ? read[top] : Number.Integer.One;
+                Entity inW = Number.Integer.Zero;
+                foreach (var (power, coefficient) in read.OrderBy(pair => pair.Key).Select(pair => (pair.Key, pair.Value)))
+                {
+                    var half = power.ToInt32Unchecked() / 2;
+                    // Folded first, so that a number the reading left as `3 (-2)` is not a name.
+                    var inNames = power == top && monic ? Number.Integer.One
+                        : Named(monic ? Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficient / leading) : coefficient.InnerSimplified);
+                    inW += half == 0 ? inNames : inNames * (half == 1 ? w : MathS.Pow(w, half));
+                }
+                return (inW, top.ToInt32Unchecked() / 2, read.Count, leading);
+            }
+            Entity denominatorInW = Number.Integer.One;
+            Entity leadingOfTheDenominator = Number.Integer.One;
+            var linearsInW = new List<Entity>();
+            var factorsInX = 0;
+            var aRepeatedBiquadratic = false;
+            foreach (var factor in Mulf.LinearChildren(denominator))
+            {
+                if (!factor.ContainsNode(x))
+                {
+                    denominatorInW *= factor;
+                    continue;
+                }
+                var (@base, power) = factor is Powf(var raised, Number.Integer { EInteger: var n }) && n.Sign > 0 && n.CanFitInInt32()
+                    ? (raised, n.ToInt32Unchecked()) : (factor, 1);
+                if (InTheSquare(@base, monic: true) is not { } read || read.Degree is not (1 or 2))
+                    return null;
+                factorsInX++;
+                // With the middle term: over a power of `a + b x^4` the split in x answers at once.
+                var repeated = read.Degree == 2 && read.Terms == 3 && power >= 2;
+                aRepeatedBiquadratic |= repeated;
+                // Every other factor once: under the tangent Rubi's `sinh(u)^4/(a - b sinh(u)^4)^2` is
+                // `w^2 (w - 1)^4/((1 - w)^3 Q^2)`, which shares the power of `1 - w` with its numerator,
+                // and split as written it was 1.46 million nodes after twenty-two seconds.
+                if (!repeated && power != 1)
+                    return null;
+                if (read.Degree == 1)
+                    linearsInW.Add(read.InW);
+                denominatorInW *= power == 1 ? read.InW : MathS.Pow(read.InW, power);
+                leadingOfTheDenominator *= power == 1 ? read.Leading : MathS.Pow(read.Leading, power);
+            }
+            if (factorsInX < 2 || !aRepeatedBiquadratic || InTheSquare(numerator, monic: false) is not { } above)
+                return null;
+            // Nor where the numerator vanishes at the root of one of the linears, which it then shares.
+            foreach (var linear in linearsInW)
+                if (TreeAnalyzer.TryGetPolyLinear(linear, w, out var slope, out var intercept)
+                    && Functions.PartialFractions.IsZeroAsAValue(above.InW.Substitute(w, -intercept / slope)))
+                    return null;
+            if (!Functions.PartialFractions.TrySplitOverWrittenFactors(above.InW, denominatorInW, w, out var decomposition, aRepeatedQuadraticIsRead: true)
+                || decomposition.Complexity > 20_000)
+                return null;
+            // The split is a sum over a constant: each term of the sum asked in x, as a question of its own.
+            var (terms, overAll) = decomposition is Divf(var sum, var by) && !by.ContainsNode(w) ? (sum, by) : (decomposition, Number.Integer.One);
+            Entity total = Number.Integer.Zero;
+            foreach (var term in Sumf.LinearChildren(terms))
+            {
+                if (Integration.ComputeAsAQuestionOfItsOwn(term.Substitute(w, MathS.Pow(x, 2)), x, integrateByParts) is not { } part)
+                    return null;
+                total += part;
+            }
+            total /= overAll * leadingOfTheDenominator;
+            foreach (var (name, value) in named)
+                total = total.Substitute(name, value);
+            return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(total, expr, x) ? total : null;
+        }
+
         /// <summary>
         /// <c>N(s)/(s^k (a + b s + c s^2))</c>, the coefficients of <c>N</c> in <paramref name="above"/>,
         /// as a polynomial in <c>s</c>, the terms in <c>s^(m - k)</c> for <c>m</c> below <c>k</c>, and

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

@@ -1423,6 +1423,8 @@ internal static bool ASumOverRootsWouldAnswer
             if ((answer = IndefiniteIntegralSolver.SolveAnEvenPolynomialOverASymbolicBiquadratic(expr, x, integrateByParts)) is { }) return answer;
             // And over a power of a linear in x^2 beside the biquadratic, split in x^2 minus its root.
             if ((answer = IndefiniteIntegralSolver.SolveAnEvenQuotientOverAPowerOfALinearInTheSquareAndABiquadratic(expr, x, integrateByParts)) is { }) return answer;
+            // And over a power of a biquadratic beside another factor in x^2, split in x^2 itself.
+            if ((answer = IndefiniteIntegralSolver.SolveAnEvenQuotientOverAPowerOfABiquadraticBesideAnotherFactorInTheSquare(expr, x, integrateByParts)) is { }) return answer;
             if ((answer = IndefiniteIntegralSolver.SolveByPartialFractions(expr, x, integrateByParts)) is { }) return answer;
             // A whole negative power of a polynomial of several terms among the factors,
             // written below the bar and asked again: the gathering on the way in writes

Sources/Tests/UnitTests/Calculus/AnEvenQuotientOverAPowerOfABiquadraticIntegralTest.cs

@@ -0,0 +1,55 @@
+//
+// 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 System.Threading;
+using System.Threading.Tasks;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// An even rational function over a power of a biquadratic of three terms beside another factor
+    /// in x^2, split in <c>w = x^2</c>: <c>1/((x^2 - 1)(c + d x^2 + f x^4)^2)</c> ran past the budget,
+    /// the series at the roots of <c>x^2 - 1</c> coming to coefficients of a hundred thousand nodes,
+    /// and so did Rubi's 6.1.7 <c>csch(u)/(a - b sinh(u)^4)^2</c>, which is it under <c>w = cosh(u)</c>.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class AnEvenQuotientOverAPowerOfABiquadraticIntegralTest
+    {
+        [Theory]
+        [InlineData("1/((x^2 - 1)*(c + d*x^2 + f*x^4)^2)")]
+        [InlineData("x^2/((x^2 + 1)*(c + d*x^2 + f*x^4)^2)")]
+        // Rubi 6.1.7:285.
+        [InlineData("csch(c + d*x)/(a - b*sinh(c + d*x)^4)^2")]
+        public void IsAnsweredWithinTheGuard(string integrand)
+        {
+            Entity? integral = null;
+            // A thread of its own, so that the guard is not spent waiting behind other tests' work.
+            var task = Task.Factory.StartNew(() => integral = integrand.ToEntity().Integrate("x"),
+                CancellationToken.None, TaskCreationOptions.LongRunning, TaskScheduler.Default);
+            Assert.True(task.Wait(LimitTermination.Guard), $"{integrand} did not finish within {LimitTermination.Guard}");
+            var text = integral!.Stringize();
+            Assert.DoesNotContain("integral(", text);
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4)
+                .Substitute("d", 0.9).Substitute("f", 0.8);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            foreach (var at in new[] { -2.2, -0.9, 0.3, 1.4, 2.1 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                Assert.False(got.IsNaN, $"d/dx of the antiderivative of {integrand} has no value at x = {at}, where the integrand is {want}");
+                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}, where the integrand is {want}");
+            }
+        }
+    }
+}