A symbolic quartic in x^2 with a square discriminant, split before the partial fractions

Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs

@@ -329,6 +329,58 @@ internal MultivariatePolynomial DerivativeIn(int variable)
             return new(VariableCount, result);
         }
 
+        /// <summary>
+        /// The polynomial whose square this is, or <see langword="null"/> where it is not the
+        /// square of one over <c>Q</c>. Of the two roots, the one whose leading coefficient in
+        /// each variable, read in turn, is positive.
+        /// </summary>
+        /// <remarks>
+        /// In the first variable that occurs, as for one variable: the root's leading coefficient
+        /// is the root of this one's, found the same way in the variables left, and each lower
+        /// coefficient is what the square of the root so far leaves at the next power down,
+        /// divided exactly by twice that leading coefficient. Whatever is not a square fails
+        /// one of those divisions or the check at the end.
+        /// </remarks>
+        internal MultivariatePolynomial? TrySquareRoot()
+        {
+            if (IsZero)
+                return this;
+            if (IsConstant)
+            {
+                var value = ConstantValue.ToLowestTerms();
+                if (value.Sign < 0)
+                    return null;
+                var (above, below) = (value.Numerator.Sqrt(), value.Denominator.Sqrt());
+                return above.Multiply(above).Equals(value.Numerator) && below.Multiply(below).Equals(value.Denominator)
+                    ? Constant(VariableCount, ERational.Create(above, below))
+                    : null;
+            }
+            var variable = 0;
+            while (DegreeIn(variable) == 0)
+                variable++;
+            var degree = DegreeIn(variable);
+            if (degree % 2 != 0)
+                return null;
+            if (LeadingCoefficientIn(variable).TrySquareRoot() is not { } leading
+                || leading.ShiftedBy(variable, degree / 2) is not { } root)
+                return null;
+            var twice = leading.ScaleBy(ERational.FromInt32(2));
+            for (var step = 1; step <= degree / 2; step++)
+            {
+                if (root.Multiply(root) is not { } square)
+                    return null;
+                var left = Subtract(square);
+                if (left.IsZero)
+                    return root;
+                if (!left.CoefficientsIn(variable).TryGetValue(degree - step, out var next))
+                    continue;
+                if (next.DivideExact(twice) is not { } term || term.ShiftedBy(variable, degree / 2 - step) is not { } shifted)
+                    return null;
+                root = root.Add(shifted);
+            }
+            return root.Multiply(root) is { } last && Subtract(last).IsZero ? root : null;
+        }
+
         internal MultivariatePolynomial LeadingCoefficientIn(int variable)
         {
             var degree = DegreeIn(variable);

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

@@ -761,6 +761,12 @@ is var (multiple, leftover)
             // And the powers of x it takes out joined to the ones beside them: `x (a x + b x^3 + c x^5)^2`
             // is `x^3 (a + b x^2 + c x^4)^2`, and written `x x^2 (a + b x^2 + c x^4)^2` the splits
             // below read `x` and `x^2` as two factors and the search ran past a minute.
+            // A symbolic quartic in x^2 that is two quadratics in it, written as them first.
+            if (WithSymbolicBiquadraticsSplit(denominator, x) is { } split
+                && (SolveByPartialFractions(numerator / split, x, integrateByParts)
+                    ?? Integration.ComputeIndefiniteIntegral(numerator / split, x, integrateByParts)) is { } overTheQuadratics)
+                return overTheQuadratics;
+
             if (WithTheContentOutOfEachSumFactor(denominator, x) is { } primitive
                 && Patterns.GatherPowersOfOneBase(numerator / primitive) is var overThePrimitives
                 && (SolveByPartialFractions(overThePrimitives, x, integrateByParts)
@@ -14965,6 +14971,70 @@ private static bool IsASumOfMonomials(Entity expr, Entity.Variable x)
             return true;
         }
 
+        /// <summary>
+        /// <paramref name="denominator"/> with each written factor <c>A x^4 + B x^2 + C</c> with symbols
+        /// in it whose discriminant <c>B^2 - 4AC</c> is the square of a polynomial <c>S</c> in them
+        /// written as the two quadratics it is, <c>(2A x^2 + B - S)(2A x^2 + B + S)/(4A)</c>;
+        /// <see langword="null"/> where there is none.
+        /// </summary>
+        /// <remarks>
+        /// The rules below read a written quadratic in <c>x</c>, and a symbolic quartic is nothing
+        /// they factor: <c>x^2/((x^2 - a)(x^4 - 2a x^2 + a^2 - b^2)^2)</c>, which the root of a linear
+        /// makes of Rubi's <c>cot(x)^3 sqrt(a + b sec(x))</c> in the secant, ran past a minute, and
+        /// written over <c>(x^2 - a - b)(x^2 - a + b)</c> it is answered in a fifth of a second.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        private static Entity? WithSymbolicBiquadraticsSplit(Entity denominator, Entity.Variable x)
+        {
+            var changed = false;
+            Entity product = Number.Integer.One;
+            foreach (var factor in Mulf.LinearChildren(denominator))
+            {
+                var (@base, power) = factor is Powf(var b, Number.Integer { EInteger.Sign: > 0 } p) ? (b, p) : (factor, Number.Integer.One);
+                if (@base is not (Sumf or Minusf) || !@base.ContainsNode(x) || !@base.Vars.Any(v => v != x)
+                    || !TreeAnalyzer.TryGetPolynomial(@base, x, out var read)
+                    || !read.Keys.All(degree => degree.Equals(EInteger.Zero) || degree.Equals(EInteger.FromInt32(2)) || degree.Equals(EInteger.FromInt32(4)))
+                    || !read.ContainsKey(EInteger.FromInt32(4)) || !read.ContainsKey(EInteger.Zero))
+                {
+                    product *= factor;
+                    continue;
+                }
+                var variables = @base.Vars.OrderBy(v => v.Name, System.StringComparer.Ordinal).ToList();
+                if (variables.Count > MultivariatePolynomial.MaxVariables)
+                {
+                    product *= factor;
+                    continue;
+                }
+                var indices = new Dictionary<Variable, int>();
+                for (var i = 0; i < variables.Count; i++)
+                    indices[variables[i]] = i;
+                var at = indices[x];
+                if (MultivariatePolynomial.TryParse(@base, indices) is not { } polynomial
+                    || !polynomial.CoefficientsIn(at).TryGetValue(4, out var a)
+                    || !polynomial.CoefficientsIn(at).TryGetValue(0, out var c))
+                {
+                    product *= factor;
+                    continue;
+                }
+                var b2 = polynomial.CoefficientsIn(at).TryGetValue(2, out var middle) ? middle : MultivariatePolynomial.Zero(variables.Count);
+                if (b2.Multiply(b2) is not { } bSquared || a.Multiply(c) is not { } ac
+                    || bSquared.Subtract(ac.ScaleBy(ERational.FromInt32(4))).TrySquareRoot() is not { IsZero: false } root
+                    || a.ScaleBy(ERational.FromInt32(2)).ShiftedBy(at, 2) is not { } twiceA)
+                {
+                    product *= factor;
+                    continue;
+                }
+                var first = twiceA.Add(b2).Subtract(root).ToEntity(variables);
+                var second = twiceA.Add(b2).Add(root).ToEntity(variables);
+                var scale = a.ScaleBy(ERational.FromInt32(4)).ToEntity(variables);
+                product *= power == Number.Integer.One
+                    ? first * second / scale
+                    : MathS.Pow(first, power) * MathS.Pow(second, power) / MathS.Pow(scale, power);
+                changed = true;
+            }
+            return changed ? product : null;
+        }
+
         /// <summary>
         /// <paramref name="denominator"/> with the content of every written sum among its
         /// factors taken out in front of it: <c>(a u + a)(1 - u^2)</c> is <c>a (u + 1)(1 - u^2)</c>,

Sources/Tests/UnitTests/Calculus/SymbolicBiquadraticSplitIntegralTest.cs

@@ -0,0 +1,49 @@
+//
+// 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 symbolic quartic in <c>x^2</c> whose discriminant is a square, written as the two
+    /// quadratics it is before the partial fractions.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class SymbolicBiquadraticSplitIntegralTest
+    {
+        [Theory]
+        [InlineData("x^2/((x^2 - a)*(x^4 - 2*a*x^2 + a^2 - b^2)^2)")]
+        [InlineData("sqrt(a + b*x)/(x*(x^2 - 1)^2)")]
+        public void OverTheTwoQuadratics(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            var text = integral.Stringize();
+            Assert.DoesNotContain("integral(", text);
+            Assert.True(text.Length < 40000, $"{text.Length} characters of answer for {integrand}");
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.4);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            var compared = 0;
+            foreach (var at in new[] { -2.1, -0.4, 0.3, 0.6, 1.6, 2.5 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                if (want.IsNaN)
+                    continue;
+                compared++;
+                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}");
+            }
+            Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}");
+        }
+    }
+}