Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))`
Expand Down Expand Up @@ -206,6 +212,164 @@ internal static class IndefiniteIntegralSolver
}
}

/// <summary>
/// The largest remainder <see cref="AsTheDerivativeOfALowerPowerOfALinearBesideAFunction"/> expands and simplifies
/// where the point gave it no value to check first.
/// </summary>
private const int LargestRemainderExpandedUnseen = 120;

/// <summary>
/// A sum whose terms are not elementary apart, integrated by parts the other way in some of the terms that are a
/// power <c>L^(-k)</c>, <c>k &gt;= 2</c>, of one linear <c>L</c> beside a function <c>F</c> with <paramref name="x"/>
/// only inside transcendental functions -- the power integrated, <c>Q = L^(1 - k)/((1 - k) L')</c>, and <c>F</c>
/// differentiated -- where what that leaves, <c>-Q F'</c>, cancels against the other terms: then the sum is
/// <c>(Q F)'</c> summed over those terms. <c>psi_1(a + b x)/x^2 - b psi_2(a + b x)/x</c> is
/// <c>(-psi_1(a + b x)/x)'</c>, Rubi's 8.6, and neither of its terms is elementary.
/// https://github.com/asc-community/AngouriMath/issues/2075
/// </summary>
/// <remarks>
/// Some, not all: in <c>psi_1(x)/(2 + 3 x)^3 - psi_2(x)/(6 (2 + 3 x)^2)</c> 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 <c>L^(-1)</c> raises the order of <c>F</c>, and a remainder in
/// <c>F'</c> is no nearer an answer.
/// </remarks>
private static Entity? AsTheDerivativeOfALowerPowerOfALinearBesideAFunction(IReadOnlyList<Entity> 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<int>(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);
}
}

/// <summary>
/// Whether a term of <paramref name="antiderivative"/> has a factor that is a polynomial
/// in <paramref name="x"/> of positive degree.
Expand Down
Original file line number Diff line number Diff line change
@@ -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
{
/// <summary>
/// A sum whose terms are not elementary apart, integrated by parts the other way in a term that is a power of a
/// linear, <c>L^(-k)</c> with <c>k &gt;= 2</c>, beside a function: the power integrated and the function differentiated,
/// where what that leaves cancels against the other terms. <c>psi_1(a + b x)/x^2 - b psi_2(a + b x)/x</c> is
/// <c>(-psi_1(a + b x)/x)'</c>. Rubi's 8.6.
/// <a href="https://github.com/asc-community/AngouriMath/issues/2075">#2075</a>
/// </summary>
[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}");
}
}

/// <summary>
/// A symbolic order, which no point evaluates: the product is read the same way, and its derivative is the
/// integrand as written.
/// </summary>
[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());
}

/// <summary>
/// <c>psi_1(x)/x^2 + psi_2(x)/x</c> is not such a derivative: what parts leaves is <c>2 psi_2(x)/x</c>, and that is not
/// elementary. Declined, not answered.
/// </summary>
[Fact]
public void ARemainderThatDoesNotCancelIsDeclined()
=> Assert.Contains("integral(", "polygamma(1, x)/x^2 + polygamma(2, x)/x".ToEntity().Integrate("x").Stringize());
}
}
Loading