Skip to content

Commit fc5d732

Browse files
Integrate a power of a function of one exponential beside x by an ansatz in that power (#2059)
Rubi's 6.7.1 sqrt(csch(x)) (x cosh(x) - 4 sech(x) tanh(x)) is the derivative of 2 x/sqrt(csch(x)) - 4 sech(x)/csch(x)^(3/2). By parts leaves -2/sqrt(csch(x)), which is elliptic, and so is the other term, so the split declined both. With t = e^(s x), a power that is not whole of a rational function of t, W = B(t)^p, beside a polynomial in x with coefficients rational in t, is W R with W'/W = p s t B_t/B rational in t. By Liouville's theorem an elementary antiderivative without a logarithmic part is W times a rational function of x and t, its denominator free of x and its degree in x R's, so F = W P/(t^m E) with P of unknown coefficients and E read off R's denominator makes F' = W R a linear system. Solved, the answer is checked by its derivative. The answer keeps the integrand's own power, and every other power of its base is written as W times a whole power of the base, which holds on every branch. The shelved respelling in powers of the hyperbolic sine wrote sqrt(1/S) as 1/sqrt(S), and was right only where sinh(x) > 0. Last in the dispatch, at the question asked. Part of #718. Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent 87be1ba commit fc5d732

4 files changed

Lines changed: 287 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -151,6 +151,21 @@ coefficients are the roots of that quadratic, and the polynomial over it is spli
151151
| `"1/(4*a*c+4*c^2*x^2+4*c*d*x^3+d^2*x^4)".ToEntity().Integrate("x")` | `integral(...)` | 3,920 characters |
152152
| `"x/(4*a*c+4*c^2*x^2+4*c*d*x^3+d^2*x^4)".ToEntity().Integrate("x")` | `integral(...)` | 2,315 characters |
153153

154+
### A power of a function of one exponential beside x is integrated by an ansatz in that power
155+
156+
**Answers where there were none.** Rubi's 6.7.1 `sqrt(csch(x)) (x cosh(x) - 4 sech(x) tanh(x))` is the derivative of
157+
`2 x/sqrt(csch(x)) - 4 sech(x)/csch(x)^(3/2)`. By parts leaves `-2/sqrt(csch(x))`, which is elliptic, and so is the
158+
other term, so the split declined both. A power that is not whole of a rational function of one exponential,
159+
`W = B(e^(s x))^p`, beside a polynomial in x whose coefficients are rational in `e^(s x)`, is now integrated by the
160+
ansatz `W` times a rational function of x and `e^(s x)`, whose coefficients are a linear system. The answer keeps the
161+
integrand's own power, holds where `sinh(x) < 0` too, and is checked by its derivative
162+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
163+
164+
| Input | Was (2.5.0) | Now |
165+
|---|---|---|
166+
| `"sqrt(csch(x))*(x*cosh(x)-4*sech(x)*tanh(x))".ToEntity().Integrate("x")` | `integral(...)` | 137 characters |
167+
| `"sqrt(csch(a+x))*((a+x)*cosh(a+x)-4*sech(a+x)*tanh(a+x))".ToEntity().Integrate("x")` | `integral(...)` | 191 characters |
168+
154169
### An even power of the sine beside the secant is written in the cosine
155170

156171
**Answers where there were none, shorter than master's.** Rubi's 4.5.1.3 `sin(c + d x)^2/(a + b sec(c + d x))` ran

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

Lines changed: 210 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30090,6 +30090,216 @@ Entity Named(Entity value)
3009030090
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
3009130091
}
3009230092

30093+
/// <summary>
30094+
/// A power that is not whole of a rational function of one exponential, <c>W = B(t)^p</c> with
30095+
/// <c>t = e^(s x)</c>, beside a polynomial in x whose coefficients are rational in t, by the ansatz
30096+
/// <c>F = W P/(t^m E)</c>: <c>P</c> a polynomial in x and t of unknown coefficients, and <c>E</c> read
30097+
/// off the denominator of the rest, <c>R = f/W</c>. Rubi's 6.7.1
30098+
/// <c>sqrt(csch(x)) (x cosh(x) - 4 sech(x) tanh(x))</c> is
30099+
/// <c>(2 x/sqrt(csch(x)) - 4 sech(x)/csch(x)^(3/2))'</c>, <c>sqrt(csch(x))</c> times a rational
30100+
/// function of x and <c>e^x</c>. By parts leaves <c>-2/sqrt(csch(x))</c>, which is elliptic, and so is
30101+
/// the other term, and the split declined them both.
30102+
/// </summary>
30103+
/// <remarks>
30104+
/// <para>
30105+
/// <c>W'/W = p s t B_t/B</c> is rational in t, so <c>F' = W R</c> is
30106+
/// <c>(P/(t^m E))' + p s t (B_t/B) P/(t^m E) = R</c>, a polynomial identity in x and t over one
30107+
/// denominator, linear in the unknowns. By Liouville's theorem an elementary antiderivative with no
30108+
/// logarithmic part is W times a rational function of x and t, with a denominator free of x, as R's
30109+
/// is, and a degree in x no higher than R's: a higher coefficient would be a constant over W, which
30110+
/// is not rational. Its denominator is R's with each factor's power lowered by one, or R's itself,
30111+
/// with t to any power, as <see cref="SolveByRischNormanAnsatz"/> bounds its own.
30112+
/// </para>
30113+
/// <para>
30114+
/// W is the integrand's own power, and every other power of its base is written as W times a whole
30115+
/// power of the base, which holds on every branch: <c>B^r = B^p B^(r - p)</c> for a whole
30116+
/// <c>r - p</c>. Respelled in powers of the hyperbolic sine, <c>sqrt(1/S)</c> as <c>1/sqrt(S)</c>, the
30117+
/// same integrand was right only where <c>sinh(x) > 0</c>. Last, at the question asked, and checked
30118+
/// by its derivative.
30119+
/// </para>
30120+
/// https://github.com/asc-community/AngouriMath/issues/718
30121+
/// </remarks>
30122+
internal static Entity? SolveAPowerOfAFunctionOfAnExponentialBesideXByAnAnsatz(Entity expr, Entity.Variable x)
30123+
{
30124+
if (!Integration.AnsweringTheQuestionAsked)
30125+
return null;
30126+
// Every exponential of something of x is e^(c + k s x), for one s and a whole k.
30127+
Entity? slope = null;
30128+
var exponentials = new List<(Entity Node, Entity Intercept, int Times)>();
30129+
foreach (var node in expr.Nodes)
30130+
{
30131+
if (node is not Powf(var e, var exponent) || e != MathS.e || !exponent.ContainsNode(x))
30132+
continue;
30133+
if (!TreeAnalyzer.TryGetPolyLinear(exponent, x, out var s, out var c) || s.ContainsNode(x) || c.ContainsNode(x) || TreeAnalyzer.IsZero(s))
30134+
return null;
30135+
slope ??= s;
30136+
if (Functions.PartialFractions.InLowestTermsOverTheSymbols(s / slope) is not Number.Integer { EInteger: var k }
30137+
|| !k.CanFitInInt32() || k.Abs().CompareTo(EInteger.FromInt32(MaximumMultipleOfTheExponent)) > 0)
30138+
return null;
30139+
if (!exponentials.Any(known => known.Node == node))
30140+
exponentials.Add((node, c, k.ToInt32Unchecked()));
30141+
}
30142+
if (slope is null)
30143+
return null;
30144+
var t = Variable.CreateUnique(expr, "t_w");
30145+
Entity InT(Entity e) => e.Replace(node =>
30146+
{
30147+
var at = exponentials.FindIndex(known => known.Node == node);
30148+
if (at < 0)
30149+
return node;
30150+
var (_, intercept, times) = exponentials[at];
30151+
var magnitude = System.Math.Abs(times);
30152+
Entity raised = magnitude == 1 ? t : MathS.Pow(t, magnitude);
30153+
Entity power = times > 0 ? raised : 1 / raised;
30154+
return TreeAnalyzer.IsZero(intercept) ? power : MathS.Pow(MathS.e, intercept) * power;
30155+
});
30156+
var inT = InT(expr);
30157+
// With x itself beside the exponential: a function of the exponential alone is the substitution's.
30158+
if (!inT.ContainsNode(x))
30159+
return null;
30160+
// One base under powers that are not whole, free of x but for the exponential, the powers a whole
30161+
// number apart; the base by its value in lowest terms, since two spellings of `csch(x)` that print
30162+
// alike need not be one tree.
30163+
Entity? radicand = null;
30164+
Number.Rational? power = null;
30165+
var spellings = new List<Entity>();
30166+
foreach (var node in inT.Nodes)
30167+
{
30168+
if (node is not Powf(var @base, var exponent) || exponent is Number.Integer || !@base.ContainsNode(t) && !@base.ContainsNode(x))
30169+
continue;
30170+
if (exponent is not Number.Rational r || @base.ContainsNode(x))
30171+
return null;
30172+
if (radicand is null)
30173+
{
30174+
(radicand, power) = (@base, r);
30175+
spellings.Add(@base);
30176+
continue;
30177+
}
30178+
if (!r.ERational.Subtract(power!.ERational).Denominator.Equals(EInteger.One))
30179+
return null;
30180+
if (spellings.Contains(@base))
30181+
continue;
30182+
if (Functions.PartialFractions.InLowestTermsOverTheSymbols(Functions.SingleQuotient.Combine(@base - radicand)) != Number.Integer.Zero)
30183+
return null;
30184+
spellings.Add(@base);
30185+
}
30186+
if (radicand is null || power is null)
30187+
return null;
30188+
var written = expr.Nodes.FirstOrDefault(node => node is Powf(var b, var r) && r == power && b.ContainsNode(x) && spellings.Contains(InT(b)));
30189+
if (written is null)
30190+
return null;
30191+
// The rest, each power of the base written over W, a whole power of the base.
30192+
var rest = inT.Replace(node => node is Powf(var @base, Number.Rational r) && r is not Number.Integer && spellings.Contains(@base)
30193+
? r.ERational.Subtract(power.ERational) is var whole && whole.IsZero ? Number.Integer.One
30194+
: whole.Equals(ERational.One) ? radicand : MathS.Pow(radicand, Number.Rational.Create(whole))
30195+
: node);
30196+
var variables = new List<Variable> { x, t };
30197+
var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(rest).InnerSimplified);
30198+
if (Functions.PolynomialGcd.TryCancel(above.Expand().InnerSimplified, below.Expand().InnerSimplified, out var cancelled) && cancelled is not null)
30199+
(above, below) = Functions.SingleQuotient.Of(Functions.PartialFractions.Bare(cancelled));
30200+
if (!TryReadInAll(above, variables, out var aboveRead) || !TryReadInAll(below, variables, out var belowRead)
30201+
|| belowRead.Keys.Any(key => key[0] > 0))
30202+
return null;
30203+
var degreeInX = aboveRead.Keys.Max(key => key[0]);
30204+
var aboveInT = aboveRead.Keys.Max(key => key[1]);
30205+
var belowInT = belowRead.Keys.Max(key => key[1]);
30206+
if (degreeInX > MaximumDegreeInXBesideAnExponential || aboveInT > 2 * MaximumDegreeInTheExponential || belowInT > 2 * MaximumDegreeInTheExponential)
30207+
return null;
30208+
var (radicandAbove, radicandBelow) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(radicand).InnerSimplified);
30209+
if (!TryReadInAll(radicandAbove, variables, out var baseAbove) || !TryReadInAll(radicandBelow, variables, out var baseBelow)
30210+
|| baseAbove.Keys.Concat(baseBelow.Keys).Any(key => key[0] > 0))
30211+
return null;
30212+
// W'/W = p s t (B_n' B_d - B_n B_d')/(B_n B_d), in t.
30213+
var betaAbove = power * slope * t * (radicandAbove.Differentiate(t) * radicandBelow - radicandAbove * radicandBelow.Differentiate(t));
30214+
var betaBelow = radicandAbove * radicandBelow;
30215+
// The denominator's factors in t with their powers, t itself apart: it may stand to any power.
30216+
var factors = new List<(Entity Factor, int Power)>();
30217+
var powerOfT = 0;
30218+
foreach (var factor in Mulf.LinearChildren(below))
30219+
{
30220+
var (@base, raisedTo) = factor is Powf(var inner, Number.Integer whole) && whole.EInteger.CanFitInInt32() ? (inner, whole.EInteger.ToInt32Unchecked()) : (factor, 1);
30221+
if (!@base.ContainsNode(t))
30222+
continue;
30223+
if (@base == t)
30224+
{
30225+
powerOfT += raisedTo;
30226+
continue;
30227+
}
30228+
var index = factors.FindIndex(known => known.Factor == @base);
30229+
if (index < 0)
30230+
factors.Add((@base, raisedTo));
30231+
else
30232+
factors[index] = (@base, factors[index].Power + raisedTo);
30233+
}
30234+
var candidates = new[] { true, false }.Select(lowered => factors.Aggregate((Entity)Number.Integer.One, (product, pair) =>
30235+
(lowered ? pair.Power - 1 : pair.Power) is var p && p > 0 ? product * (p == 1 ? pair.Factor : MathS.Pow(pair.Factor, p)) : product)).Distinct().ToList();
30236+
foreach (var q in candidates)
30237+
for (var m = 0; m <= powerOfT + 1; m++)
30238+
{
30239+
var denominator = m == 0 ? q : q * (m == 1 ? t : MathS.Pow(t, m));
30240+
if (!TreeAnalyzer.TryGetPolynomial(denominator, t, out var denominatorRead) || denominatorRead.Count == 0)
30241+
continue;
30242+
var degreeInT = System.Math.Min(MaximumDegreeInTheExponential,
30243+
System.Math.Max(0, aboveInT - belowInT + denominatorRead.Keys.Max()!.ToInt32Unchecked() + 1));
30244+
if ((degreeInX + 1) * (degreeInT + 1) > MaximumUnknownsBesideAnExponential)
30245+
continue;
30246+
var denominatorPrime = denominator.Differentiate(t);
30247+
var generators = new List<Entity>();
30248+
var columns = new List<Dictionary<int[], Entity>>();
30249+
var failed = false;
30250+
for (var i = 0; i <= degreeInX && !failed; i++)
30251+
for (var j = 0; j <= degreeInT; j++)
30252+
{
30253+
Entity generator = i == 0 && j == 0 ? Number.Integer.One
30254+
: i == 0 ? MathS.Pow(t, j) : j == 0 ? MathS.Pow(x, i) : MathS.Pow(x, i) * MathS.Pow(t, j);
30255+
// (g/E)' + (W'/W) g/E over E^2 B_n B_d D_R, with g' = g_x + s t g_t.
30256+
var numerator = ((generator.Differentiate(x) + slope * t * generator.Differentiate(t)) * denominator
30257+
- slope * t * generator * denominatorPrime) * betaBelow * below
30258+
+ betaAbove * generator * denominator * below;
30259+
if (!TryReadInAll(Functions.PartialFractions.Bare(numerator.Expand().InnerSimplified), variables, out var column))
30260+
{
30261+
failed = true;
30262+
break;
30263+
}
30264+
generators.Add(generator);
30265+
columns.Add(column);
30266+
}
30267+
if (failed)
30268+
continue;
30269+
if (!TryReadInAll(Functions.PartialFractions.Bare((above * MathS.Sqr(denominator) * betaBelow).Expand().InnerSimplified), variables, out var target))
30270+
continue;
30271+
var keys = columns.SelectMany(column => column.Keys).Concat(target.Keys).Distinct(new IntArrayComparer()).ToList();
30272+
var matrix = keys.Select(key => columns.Select(column => column.TryGetValue(key, out var entry) ? entry : Number.Integer.Zero).ToArray()).ToArray();
30273+
var rhs = keys.Select(key => target.TryGetValue(key, out var wanted) ? wanted : Number.Integer.Zero).ToArray();
30274+
if ((!Functions.PartialFractions.TrySolveLinear(matrix, rhs, out var values) || values is null)
30275+
&& (!Functions.PartialFractions.TrySolveLinearWithSymbols(matrix, rhs, out values) || values is null))
30276+
continue;
30277+
Entity? sum = null;
30278+
for (var c = 0; c < values.Length; c++)
30279+
{
30280+
var value = values[c].InnerSimplified;
30281+
if (value.Evaled is Number.Complex { IsZero: true })
30282+
continue;
30283+
sum = sum is null ? value * generators[c] : sum + value * generators[c];
30284+
}
30285+
if (sum is null)
30286+
continue;
30287+
var answer = written * (sum / denominator).Substitute(t, MathS.Pow(MathS.e, slope * x));
30288+
if (Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x))
30289+
return answer;
30290+
}
30291+
return null;
30292+
}
30293+
30294+
/// <summary>
30295+
/// The largest multiple of the exponent <see cref="SolveAPowerOfAFunctionOfAnExponentialBesideXByAnAnsatz"/>
30296+
/// reads, the largest degrees of its ansatz in x and in the exponential, and the most unknowns.
30297+
/// </summary>
30298+
private const int MaximumMultipleOfTheExponent = 8;
30299+
private const int MaximumDegreeInXBesideAnExponential = 3;
30300+
private const int MaximumDegreeInTheExponential = 12;
30301+
private const int MaximumUnknownsBesideAnExponential = 60;
30302+
3009330303
/// <summary>
3009430304
/// A logarithm of constants and nested powers of one linear <c>L</c>, <c>ln(c (d L^p)^q)</c>, is
3009530305
/// <c>p q ln(L) + K</c> with <c>K</c> constant wherever it is continuous, since its derivative is

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1818,6 +1818,10 @@ internal static bool ASumOverRootsWouldAnswer
18181818
// And half-odd powers of a hyperbolic sine beside an exponential as the derivative of the same kind,
18191819
// last: each term alone is not elementary, and the split declined them.
18201820
if ((answer = IndefiniteIntegralSolver.SolveHalfOddPowersOfAHyperbolicSineBesideAnExponentialAsADerivative(expr, x)) is { }) return answer;
1821+
// And a power that is not whole of a function of one exponential beside x, by an ansatz in that
1822+
// power times a rational function of x and the exponential, last: by parts and the split leave
1823+
// terms that are not elementary apart.
1824+
if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAFunctionOfAnExponentialBesideXByAnAnsatz(expr, x)) is { }) return answer;
18211825
return null;
18221826
}
18231827
/// <summary>
Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
//
2+
// Copyright (c) 2019-2026 Angouri.
3+
// AngouriMath is licensed under MIT.
4+
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
5+
// Website: https://am.angouri.org.
6+
//
7+
8+
using System;
9+
using AngouriMath.Extensions;
10+
using Xunit;
11+
12+
namespace AngouriMath.Tests.Calculus
13+
{
14+
/// <summary>
15+
/// A power that is not whole of a function of one exponential, beside a polynomial in x, by the
16+
/// ansatz <c>F = W P/(t^m E)</c> with <c>t = e^(s x)</c>: Rubi's 6.7.1
17+
/// <c>sqrt(csch(x)) (x cosh(x) - 4 sech(x) tanh(x))</c> is
18+
/// <c>(2 x/sqrt(csch(x)) - 4 sech(x)/csch(x)^(3/2))'</c>, and by parts and the split left terms that are
19+
/// not elementary apart. The answer keeps the integrand's own power, and holds where <c>sinh(x) &lt; 0</c>
20+
/// too.
21+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
22+
/// </summary>
23+
[Trait("Area", "Calculus")]
24+
public sealed class APowerOfAFunctionOfAnExponentialBesideXIntegralTest
25+
{
26+
[Theory]
27+
// Rubi 6.7.1:1243.
28+
[InlineData("sqrt(csch(x))*(x*cosh(x) - 4*sech(x)*tanh(x))")]
29+
[InlineData("sqrt(csch(a + x))*((a + x)*cosh(a + x) - 4*sech(a + x)*tanh(a + x))")]
30+
public void IsAnsweredOnBothSidesOfZero(string integrand)
31+
=> Check(integrand.ToEntity());
32+
33+
[Fact]
34+
public void AnswersTheDerivativeOfOneOfTheSameShape()
35+
// Quadratic in x: (sqrt(csch(x)) (x^2 sinh(x) + x cosh(x)))', whose terms are not elementary apart.
36+
=> Check("sqrt(csch(x))*(x^2*sinh(x) + x*cosh(x))".ToEntity().Differentiate("x"));
37+
38+
private static void Check(Entity integrand)
39+
{
40+
var integral = integrand.Integrate("x");
41+
var text = integral.Stringize();
42+
Assert.DoesNotContain("integral(", text);
43+
Assert.True(text.Length < 2000, $"{integrand} came back in {text.Length} characters");
44+
Entity Pinned(Entity e) => e.Substitute("a", 0.3);
45+
var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
46+
var original = Pinned(integrand);
47+
foreach (var at in new[] { -2.2, -0.9, 0.4, 1.4, 2.1 })
48+
{
49+
var want = original.Substitute("x", at).EvalNumerical();
50+
var got = derivative.Substitute("x", at).EvalNumerical();
51+
Assert.False(got.IsNaN, $"d/dx of the antiderivative of {integrand} has no value at x = {at}, where the integrand is {want}");
52+
Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
53+
< 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)),
54+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
55+
}
56+
}
57+
}
58+
}

0 commit comments

Comments
 (0)