diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index e00250237..f9f38feec 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -46,6 +46,9 @@ internal static class IndefiniteIntegralSolver if (asWritten.Count >= 2 && asWritten.All(term => term.Complexity <= LargestTermTakenAsWritten) && Integrated(asWritten) is { } termByTerm) return termByTerm; + // As written, where a power of a linear is not yet expanded into a polynomial; it asks nothing. + if (asWritten.Count >= 2 && AsTheDerivativeOfALowerPowerOfALinearBesideAFunction(asWritten, x) is { } writtenAsADerivative) + return writtenAsADerivative; } // A sum above the bar that is the derivative of a factor below it is not expanded: // `(e^x + e^-x)/(e^x - e^-x)^4` -- `cosh/sinh^4` as the library writes it -- is @@ -94,6 +97,9 @@ internal static class IndefiniteIntegralSolver }).ToList(); if (Integrated(cancelledTerms) is { } termByExpandedTerm) return termByExpandedTerm; + // Asks nothing, so at any depth. + if (AsTheDerivativeOfALowerPowerOfALinearBesideAFunction(splitted, x) is { } exactDerivative) + return exactDerivative; // The expanded terms gathered again by what is not polynomial in them, and each // group asked as one term: Hearn's // `(-2 sqrt(1 + x^3) + 5x^4 sqrt(1 + x^3) - 3x^2 sqrt(1 - 2x + x^5))/(2 sqrt(1 + x^3) sqrt(1 - 2x + x^5))` @@ -206,6 +212,164 @@ internal static class IndefiniteIntegralSolver } } + /// + /// The largest remainder expands and simplifies + /// where the point gave it no value to check first. + /// + private const int LargestRemainderExpandedUnseen = 120; + + /// + /// A sum whose terms are not elementary apart, integrated by parts the other way in some of the terms that are a + /// power L^(-k), k >= 2, of one linear L beside a function F with + /// only inside transcendental functions -- the power integrated, Q = L^(1 - k)/((1 - k) L'), and F + /// differentiated -- where what that leaves, -Q F', cancels against the other terms: then the sum is + /// (Q F)' summed over those terms. psi_1(a + b x)/x^2 - b psi_2(a + b x)/x is + /// (-psi_1(a + b x)/x)', Rubi's 8.6, and neither of its terms is elementary. + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// + /// + /// Some, not all: in psi_1(x)/(2 + 3 x)^3 - psi_2(x)/(6 (2 + 3 x)^2) the second term has a power of -2 as well, + /// and is what the first leaves. So each set of them is tried, the smallest first, and rejected at a point by adding + /// numbers before anything is simplified. Nothing is asked, so a remainder that does not cancel declines rather than + /// becoming a question: lowering the power toward L^(-1) raises the order of F, and a remainder in + /// F' is no nearer an answer. + /// + private static Entity? AsTheDerivativeOfALowerPowerOfALinearBesideAFunction(IReadOnlyList terms, Entity.Variable x) + { + // Each such term with what parts the other way makes of it: Q F, and the -Q F' it leaves. + var steps = new List<(int Term, Entity Parts, Entity Leaves)>(); + for (var i = 0; i < terms.Count; i++) + { + Entity coefficient = Number.Integer.One; + Entity function = Number.Integer.One; + Entity? linear = null; + var power = EInteger.Zero; + var read = true; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(terms[i])) + { + if (!factor.ContainsNode(x)) + { + coefficient = underneath ? coefficient / factor : coefficient * factor; + continue; + } + var (@base, exponent) = factor is Powf(var b, Number.Integer n) ? (b, n.EInteger) : (factor, EInteger.One); + if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slopeOfIt, out _) && !TreeAnalyzer.IsZero(slopeOfIt)) + { + if (linear is not null && linear != @base) + { + read = false; + break; + } + linear = @base; + power = underneath ? power.Subtract(exponent) : power.Add(exponent); + } + else if (OnlyInsideTranscendentalFunctions(factor, x)) + function = underneath ? function / factor : function * factor; + else + { + read = false; + break; + } + } + // And a function differentiable where it is defined: (Q F)' is Q' F + Q F' only there, and across the jump of + // a sign it is not, as the arcsecant's antiderivative, sgn(c x) ln(c x + sqrt(c^2 x^2 - 1)), has. + if (!read || linear is null || power.CompareTo(EInteger.FromInt32(-2)) > 0 || !function.ContainsNode(x) + || function.Nodes.Any(node => node is Signumf or Absf or Piecewise or Providedf) + || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var slope, out _)) + continue; + // Four at most: fifteen sets. + if (steps.Count == 4) + return null; + // Q F = -c F/((k - 1) L' L^(k - 1)), with no negative power to write, and -Q F' = c F'/((k - 1) L' L^(k - 1)). + var lowered = power.Negate().Subtract(EInteger.One); + var below = (Number.Integer.Create(lowered) * slope).InnerSimplified + * (lowered.Equals(EInteger.One) ? linear : MathS.Pow(linear, Number.Integer.Create(lowered))); + var derivative = function.Differentiate(x); + if (derivative is Providedf(var bare, _)) + derivative = bare; + steps.Add((i, -(coefficient * function) / below, coefficient * derivative / below)); + } + if (steps.Count == 0) + return null; + + // Every piece at one point, the same for all: a value each term has, and each step leaves. + var symbols = terms.SelectMany(term => term.Vars).Concat(steps.SelectMany(step => step.Leaves.Vars)).Distinct().ToList(); + var values = new[] { 1.37, 0.61, 2.23, 1.91, 0.83, 1.13 }; + System.Numerics.Complex? AtAPoint(Entity piece) + { + var pinned = piece; + for (var k = 0; k < symbols.Count; k++) + pinned = pinned.Substitute(symbols[k], values[k % values.Length]); + if (pinned.Evaled is not Number.Complex { IsNaN: false } value) + return null; + var number = new System.Numerics.Complex(value.RealPart.EDecimal.ToDouble(), value.ImaginaryPart.EDecimal.ToDouble()); + return double.IsInfinity(number.Real) || double.IsInfinity(number.Imaginary) || double.IsNaN(number.Real) || double.IsNaN(number.Imaginary) + ? null : number; + } + var termsAtAPoint = terms.Select(AtAPoint).ToList(); + var leavesAtAPoint = steps.Select(step => AtAPoint(step.Leaves)).ToList(); + foreach (var set in Enumerable.Range(1, (1 << steps.Count) - 1).OrderBy(set => Enumerable.Range(0, steps.Count).Count(k => (set >> k & 1) == 1))) + { + var taken = new HashSet(Enumerable.Range(0, steps.Count).Where(k => (set >> k & 1) == 1).Select(k => steps[k].Term)); + var atAPoint = System.Numerics.Complex.Zero; + var scale = 0.0; + var known = true; + for (var i = 0; i < terms.Count && known; i++) + if (!taken.Contains(i)) + known = Add(termsAtAPoint[i]); + for (var k = 0; k < steps.Count && known; k++) + if ((set >> k & 1) == 1) + known = Add(leavesAtAPoint[k]); + if (known && atAPoint.Magnitude > 1e-12 * (1 + scale)) + continue; + Entity left = Number.Integer.Zero; + Entity? parts = null; + for (var i = 0; i < terms.Count; i++) + if (!taken.Contains(i)) + left += terms[i]; + for (var k = 0; k < steps.Count; k++) + if ((set >> k & 1) == 1) + { + left += steps[k].Leaves; + parts = parts is null ? steps[k].Parts : parts + steps[k].Parts; + } + if (parts is not null && CancelsIdentically(left, known)) + return parts.InnerSimplified; + + bool Add(System.Numerics.Complex? value) + { + if (value is not { } number) + return false; + atAPoint += number; + scale = System.Math.Max(scale, number.Magnitude); + return true; + } + } + return null; + + // InnerSimplified leaves psi_2 b/x - b psi_2/x as two terms, and expanded they cancel. Where the point gave no + // value -- a function of a symbolic order, as polygamma(n, a + b x) is -- only a small remainder is expanded and + // simplified: doing so to what the arcsecant's integral leaves took eight seconds a time, and the integral + // timed out. + static bool CancelsIdentically(Entity remainder, bool numericallyZero) + { + var inner = remainder.InnerSimplified; + if (TreeAnalyzer.IsZero(inner)) + return true; + if (!numericallyZero && inner.Complexity > LargestRemainderExpandedUnseen) + return false; + var expanded = inner.Expand().InnerSimplified; + if (expanded is Providedf(var bareExpanded, _)) + expanded = bareExpanded; + if (TreeAnalyzer.IsZero(expanded)) + return true; + var simplified = inner.Simplify(); + if (simplified is Providedf(var bareSimplified, _)) + simplified = bareSimplified; + return TreeAnalyzer.IsZero(simplified); + } + } + /// /// Whether a term of has a factor that is a polynomial /// in of positive degree. diff --git a/Sources/Tests/UnitTests/Calculus/ASumThatIsTheDerivativeOfALowerPowerBesideAFunctionTest.cs b/Sources/Tests/UnitTests/Calculus/ASumThatIsTheDerivativeOfALowerPowerBesideAFunctionTest.cs new file mode 100644 index 000000000..6065eea04 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumThatIsTheDerivativeOfALowerPowerBesideAFunctionTest.cs @@ -0,0 +1,70 @@ +// +// 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 sum whose terms are not elementary apart, integrated by parts the other way in a term that is a power of a + /// linear, L^(-k) with k >= 2, beside a function: the power integrated and the function differentiated, + /// where what that leaves cancels against the other terms. psi_1(a + b x)/x^2 - b psi_2(a + b x)/x is + /// (-psi_1(a + b x)/x)'. Rubi's 8.6. + /// #2075 + /// + [Trait("Area", "Calculus")] + public sealed class ASumThatIsTheDerivativeOfALowerPowerBesideAFunctionTest + { + [Theory] + [InlineData("polygamma(1, a + b*x)/x^2 - b*polygamma(2, a + b*x)/x", "-polygamma(1, a + b*x)/x")] + [InlineData("digamma(a + b*x)/x^2 - b*polygamma(1, a + b*x)/x", "-digamma(a + b*x)/x")] + // The second term has a power of -2 as well, and is what the first leaves. + [InlineData("polygamma(1, x)/(a + b*x)^3 - polygamma(2, x)/(2*b*(a + b*x)^2)", "-polygamma(1, x)/(2*b*(a + b*x)^2)")] + // Two such derivatives at once. + [InlineData("polygamma(1, x)/x^2 - polygamma(2, x)/x + sin(x)/x^3 - cos(x)/(2*x^2)", "-polygamma(1, x)/x - sin(x)/(2*x^2)")] + public void IsTheDerivativeOfTheProduct(string integrand, string expected) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + var difference = integral.Substitute("C", 0) - expected.ToEntity(); + var derivative = integral.Substitute("C", 0).Differentiate("x"); + foreach (var (a, b, x) in new[] { (1.3, 0.7, 0.9), (2.1, 1.9, 1.6), (0.6, 2.4, 3.1) }) + { + Entity At(Entity e) => e.Substitute("a", a).Substitute("b", b).Substitute("x", x); + var off = At(difference).EvalNumerical(); + Assert.True(Math.Abs((double)off.RealPart) + Math.Abs((double)off.ImaginaryPart) < 1e-12, + $"the antiderivative of {integrand} is {integral}, not {expected}, at a = {a}, b = {b}, x = {x}"); + var want = At(integrand.ToEntity()).EvalNumerical(); + var got = At(derivative).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {x}, where the integrand is {want}"); + } + } + + /// + /// A symbolic order, which no point evaluates: the product is read the same way, and its derivative is the + /// integrand as written. + /// + [Fact] + public void ASymbolicOrder() + { + var integral = "polygamma(n, a + b*x)/x^2 - b*polygamma(1 + n, a + b*x)/x".ToEntity().Integrate("x"); + Assert.Equal("-polygamma(n, a + b * x) / x + C", integral.Stringize()); + } + + /// + /// psi_1(x)/x^2 + psi_2(x)/x is not such a derivative: what parts leaves is 2 psi_2(x)/x, and that is not + /// elementary. Declined, not answered. + /// + [Fact] + public void ARemainderThatDoesNotCancelIsDeclined() + => Assert.Contains("integral(", "polygamma(1, x)/x^2 + polygamma(2, x)/x".ToEntity().Integrate("x").Stringize()); + } +}