diff --git a/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs b/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs index 51e62e832..7744f6dc2 100644 --- a/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs +++ b/Sources/AngouriMath/Functions/Algebra/Polynomials/MultivariatePolynomial.cs @@ -329,6 +329,58 @@ internal MultivariatePolynomial DerivativeIn(int variable) return new(VariableCount, result); } + /// + /// The polynomial whose square this is, or where it is not the + /// square of one over Q. Of the two roots, the one whose leading coefficient in + /// each variable, read in turn, is positive. + /// + /// + /// 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. + /// + 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); diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 8811b1000..eb29c904c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/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; } + /// + /// with each written factor A x^4 + B x^2 + C with symbols + /// in it whose discriminant B^2 - 4AC is the square of a polynomial S in them + /// written as the two quadratics it is, (2A x^2 + B - S)(2A x^2 + B + S)/(4A); + /// where there is none. + /// + /// + /// The rules below read a written quadratic in x, and a symbolic quartic is nothing + /// they factor: x^2/((x^2 - a)(x^4 - 2a x^2 + a^2 - b^2)^2), which the root of a linear + /// makes of Rubi's cot(x)^3 sqrt(a + b sec(x)) in the secant, ran past a minute, and + /// written over (x^2 - a - b)(x^2 - a + b) it is answered in a fifth of a second. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + 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(); + 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; + } + /// /// with the content of every written sum among its /// factors taken out in front of it: (a u + a)(1 - u^2) is a (u + 1)(1 - u^2), diff --git a/Sources/Tests/UnitTests/Calculus/SymbolicBiquadraticSplitIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/SymbolicBiquadraticSplitIntegralTest.cs new file mode 100644 index 000000000..cb49af497 --- /dev/null +++ b/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 +{ + /// + /// A symbolic quartic in x^2 whose discriminant is a square, written as the two + /// quadratics it is before the partial fractions. + /// #718 + /// + [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}"); + } + } +}