Skip to content

Commit e8f9ade

Browse files
Integrate x^m (d + e x^2)^p times an inverse trigonometric function by parts (#1973)
x^m (d + e x^2)^p (a + b F(c x)), with F an inverse trigonometric function and m + 2p + 3 = 0, is integrated by parts against the antiderivative of x^m (d + e x^2)^p that is one term, v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1)): the integral is v (a + b F(c x)) - b ∫ v F(c x)', and F(c x)' is algebraic. Rubi's 5.1.4:789, 5.3.4:1399, 5.5.1:144 and 5.6.1:149 were past the budget. Part of #233. Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent 607fd57 commit e8f9ade

4 files changed

Lines changed: 217 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6637,6 +6637,25 @@ leaves are read as quotients of the sine and the cosine
66376637
| `"1/(a+b*asin(c*x))".Integrate("x")` | `integral(1 / (a + b * arcsin(c * x)), x)` | `Si` and `Ci` of `(a + b arcsin(c x))/b` |
66386638
| `"x/((c+a^2*c*x^2)^2*atan(a*x))".Integrate("x")` | `integral(x / ((c + a ^ 2 * c * x ^ 2) ^ 2 * arctan(a * x)), x)` | `1/2 * Si(2 * arctan(a * x)) / (a ^ 2 * c ^ 2) + C` |
66396639

6640+
### A power of x over a power of a quadratic times an inverse trigonometric function is integrated by parts against one term
6641+
6642+
**Answers where there were none.** `x^3 (a + b F(c x))/(d + e x^2)^3` for `F` the inverse sine, tangent, secant and
6643+
cosecant, Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1, was declined. By parts against the integrator's own antiderivative of
6644+
`x^3/(d + e x^2)^3`, a piecewise on `e = 0`, spent the budget on every arm. Where `m + 2p + 3 = 0`,
6645+
`x^m (d + e x^2)^p` has the one-term antiderivative `v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1))`, Rubi's, which
6646+
vanishes at 0, and the integral is now `v (a + b F(c x)) - b integral v F(c x)'` with `F(c x)'` algebraic.
6647+
`v' = x^m (d + e x^2)^p` is an identity of rational functions, so the parts hold wherever the integrand is defined,
6648+
on any branch of `F`. Only
6649+
for a whole `p` of `-3` or less, and not where `d + e x^2` is a multiple of the radicand `F` holds, which the
6650+
substitution under `F` answers ([#233](https://github.com/asc-community/AngouriMath/issues/233)).
6651+
6652+
| Input | Was (2.5.0) | Now |
6653+
|---|---|---|
6654+
| `"x^3*(a + b*arcsin(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 6,013 characters |
6655+
| `"x^3*(a + b*arctan(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,066 characters |
6656+
| `"x^3*(a + b*arcsec(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 1,772 characters |
6657+
| `"x^5*(a + b*arctan(c*x))/(d + k*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | 6,176 characters |
6658+
66406659
### The arcsecant beside a root of `x^2 - 1` is integrated right below `-1` too
66416660

66426661
**Answers where there were none.** Under `x = sec(u)`, the substitution that undoes the arcsecant,

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

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32357,6 +32357,134 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
3235732357
return (alpha / (p0 * slope) * product).InnerSimplified;
3235832358
}
3235932359

32360+
/// <summary>
32361+
/// <c>x^m (d + e x^2)^p (a + b F(c x))</c>, <c>F</c> an inverse trigonometric function and
32362+
/// <c>m + 2p + 3 = 0</c>, by parts against the antiderivative of <c>x^m (d + e x^2)^p</c> that is
32363+
/// one term, <c>v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1))</c>:
32364+
/// <c>integral = v (a + b F(c x)) - b integral v F(c x)'</c>, and <c>F(c x)'</c> is algebraic.
32365+
/// </summary>
32366+
/// <remarks>
32367+
/// <para>
32368+
/// <c>v' = x^m (d + e x^2)^p ((m + 1)(d + e x^2) + 2 (p + 1) e x^2)/(d (m + 1))</c>, and the bracket
32369+
/// is <c>(m + 1) d + (m + 2p + 3) e x^2</c>, which is <c>(m + 1) d</c> for every <c>e</c>: an
32370+
/// identity of rational functions, so the parts are exact wherever the integrand is defined, on
32371+
/// any branch of <c>F</c>. It is Rubi's choice of <c>v</c>, and the one that vanishes at
32372+
/// <c>x = 0</c>.
32373+
/// </para>
32374+
/// <para>
32375+
/// The integrator's own antiderivative of <c>x^3/(d + e x^2)^3</c> is a piecewise on <c>e = 0</c>,
32376+
/// and by parts against it spent the budget on every arm: Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1
32377+
/// <c>x^3 (a + b F(c x))/(d + e x^2)^3</c>, for <c>F</c> the inverse sine, tangent, secant and
32378+
/// cosecant, ran out of time, where the remainder against this <c>v</c> is answered in under a
32379+
/// second.
32380+
/// </para>
32381+
/// <para>
32382+
/// Narrow on purpose. <c>p</c> whole and at most <c>-3</c>, so that <c>v</c> is rational and
32383+
/// <c>m = -2p - 3</c> at least 3: <c>m = 1</c>, <c>x/(d + e x^2)^2</c>, is answered by parts as it
32384+
/// was, and a half-odd <c>p</c> is answered as it was beside the inverse sine and tangent and is
32385+
/// elliptic beside the secant. And not for <c>e = +-c^2 d</c>, the quadratic a multiple of
32386+
/// <c>F</c>'s own radicand, which the substitution under <c>F</c> answers.
32387+
/// https://github.com/asc-community/AngouriMath/issues/233
32388+
/// </para>
32389+
/// </remarks>
32390+
internal static Entity? SolveAPowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByParts(Entity expr, Entity.Variable x, bool integrateByParts)
32391+
{
32392+
if (!Integration.AnsweringTheQuestionAsked || integratingAPowerOverAQuadraticBesideAnInverse)
32393+
return null;
32394+
// The one inverse trigonometric function, of c x.
32395+
Entity? function = null;
32396+
foreach (var node in expr.Nodes)
32397+
if ((node is Arcsinf or Arccosf or Arctanf or Arccotanf or Arcsecantf or Arccosecantf) && node.ContainsNode(x))
32398+
{
32399+
if (function is not null && function != node)
32400+
return null;
32401+
function = node;
32402+
}
32403+
var argument = function switch
32404+
{
32405+
Arcsinf(var u) => u, Arccosf(var u) => u, Arctanf(var u) => u,
32406+
Arccotanf(var u) => u, Arcsecantf(var u) => u, Arccosecantf(var u) => u,
32407+
_ => null
32408+
};
32409+
if (function is null || argument is null || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out var offset)
32410+
|| rate.ContainsNode(x) || !Functions.PartialFractions.IsZeroAsAValue(offset) || Functions.PartialFractions.IsZeroAsAValue(rate))
32411+
return null;
32412+
// Above the bar one factor linear in it, a + b F, and the rest; below, nothing with it.
32413+
var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
32414+
if (denominator.ContainsNode(function))
32415+
return null;
32416+
Entity? linear = null;
32417+
Entity rest = Number.Integer.One;
32418+
foreach (var factor in Mulf.LinearChildren(numerator))
32419+
{
32420+
if (!factor.ContainsNode(function))
32421+
rest *= factor;
32422+
else if (linear is null)
32423+
linear = factor;
32424+
else
32425+
return null;
32426+
}
32427+
var t = Variable.CreateUnique(expr, "t_inv");
32428+
if (linear is null || !TreeAnalyzer.TryGetPolyLinear(linear.Substitute(function, t), t, out var slope, out var intercept)
32429+
|| slope.ContainsNode(x) || intercept.ContainsNode(x) || Functions.PartialFractions.IsZeroAsAValue(slope))
32430+
return null;
32431+
// The rest over the denominator as K x^m Q^p, Q = d + e x^2.
32432+
Entity constant = Number.Integer.One;
32433+
Entity? quadratic = null;
32434+
var m = 0;
32435+
var p = 0;
32436+
foreach (var (side, sign) in new[] { (rest, 1), (denominator, -1) })
32437+
foreach (var factor in Mulf.LinearChildren(side))
32438+
{
32439+
if (!factor.ContainsNode(x))
32440+
{
32441+
constant = sign > 0 ? constant * factor : constant / factor;
32442+
continue;
32443+
}
32444+
var (@base, power) = factor is Powf(var b0, Number.Integer { EInteger: var e0 }) && e0.CanFitInInt32()
32445+
? (b0, e0.ToInt32Unchecked()) : (factor, 1);
32446+
if (@base == x)
32447+
m += sign * power;
32448+
else if (quadratic is null || quadratic == @base)
32449+
{
32450+
quadratic = @base;
32451+
p += sign * power;
32452+
}
32453+
else
32454+
return null;
32455+
}
32456+
if (quadratic is null || p > -3 || m != -2 * p - 3
32457+
|| !TreeAnalyzer.TryGetPolynomial(quadratic, x, out var read)
32458+
|| read.Keys.Any(key => !key.IsZero && !key.Equals(EInteger.FromInt32(2)))
32459+
|| !read.TryGetValue(EInteger.Zero, out var d) || !read.TryGetValue(EInteger.FromInt32(2), out var e)
32460+
|| d.ContainsNode(x) || e.ContainsNode(x)
32461+
|| Functions.PartialFractions.IsZeroAsAValue(d) || Functions.PartialFractions.IsZeroAsAValue(e)
32462+
|| Functions.PartialFractions.IsZeroAsAValue((e - rate * rate * d).Expand())
32463+
|| Functions.PartialFractions.IsZeroAsAValue((e + rate * rate * d).Expand()))
32464+
return null;
32465+
var v = constant * MathS.Pow(x, Number.Integer.Create(m + 1)) * MathS.Pow(quadratic, Number.Integer.Create(p + 1))
32466+
/ (d * (m + 1));
32467+
// The remainder holds no inverse function, so it is not this integrand again; the flag is for
32468+
// whatever its own search writes.
32469+
var remainder = (slope * v * function.Differentiate(x)).InnerSimplified;
32470+
Entity? integral;
32471+
integratingAPowerOverAQuadraticBesideAnInverse = true;
32472+
try
32473+
{
32474+
integral = Integration.ComputeAsAQuestionOfItsOwn(remainder, x, integrateByParts);
32475+
}
32476+
finally
32477+
{
32478+
integratingAPowerOverAQuadraticBesideAnInverse = false;
32479+
}
32480+
if (integral is null)
32481+
return null;
32482+
var answer = v * linear - integral;
32483+
return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer;
32484+
}
32485+
32486+
[System.ThreadStatic] private static bool integratingAPowerOverAQuadraticBesideAnInverse;
32487+
3236032488
/// <summary>
3236132489
/// <c>x^m (a + b atan(c x))^p</c> over a power of <c>d (1 +- i c x)</c>, under <c>x = tan(u)/c</c>, where
3236232490
/// <c>1 + i tan(u)</c> is <c>e^(i u)/cos(u)</c>: the integrand in <c>u</c> is

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

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -707,6 +707,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
707707
// A rational function times the logarithm of a nested product of powers of linears, by parts,
708708
// its logarithm's derivative read off as written: `ln(e (f (a + b x)^p (c + d x)^q)^r)`.
709709
if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionTimesALogarithmOfPowersOfLinears(expr, x, integrateByParts)) is { }) return answer;
710+
// x^m (d + e x^2)^p (a + b F(c x)) with m + 2p + 3 = 0 and F an inverse trigonometric function, by
711+
// parts against v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1)), before by parts against the
712+
// integrator's piecewise v timed out on `x^3 (a + b asin(c x))/(d + e x^2)^3`.
713+
if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByParts(expr, x, integrateByParts)) is { }) return answer;
710714
// A function of an inverse of something not linear times its derivative up to a factor
711715
// constant though written with x: `atan(c x/sqrt(a - c^2 x^2))^2/sqrt(d - c^2 d x^2/a)`. Early,
712716
// since it fires only where that factor holds x, and the rules before it spent 28 seconds on
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
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 System.Globalization;
10+
using System.Linq;
11+
using AngouriMath.Extensions;
12+
using Xunit;
13+
14+
namespace AngouriMath.Tests.Calculus
15+
{
16+
/// <summary>
17+
/// <c>x^m (d + k x^2)^p (a + b F(c x))</c> with <c>m + 2p + 3 = 0</c> and <c>F</c> an inverse
18+
/// trigonometric function, by parts against <c>v = x^(m + 1) (d + k x^2)^(p + 1)/(d (m + 1))</c>.
19+
/// By parts against the integrator's own antiderivative, a piecewise on <c>k = 0</c>, ran out of
20+
/// time. Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1.
21+
/// <a href="https://github.com/asc-community/AngouriMath/issues/233">#233</a>
22+
/// </summary>
23+
[Trait("Area", "Calculus")]
24+
public sealed class APowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByPartsIntegralTest
25+
{
26+
private const string Inside = "-1.3 -0.9 -0.4 0.3 0.8 1.2";
27+
private const string Outside = "-4.1 -2.9 -2.1 2.1 2.9 4.1";
28+
29+
[Theory]
30+
[InlineData("x^3*(a + b*arcsin(c*x))/(d + k*x^2)^3", Inside)]
31+
[InlineData("x^3*(a + b*arctan(c*x))/(d + k*x^2)^3", Inside)]
32+
[InlineData("x^3*(a + b*arcsec(c*x))/(d + k*x^2)^3", Outside)]
33+
[InlineData("x^3*(a + b*arccsc(c*x))/(d + k*x^2)^3", Outside)]
34+
[InlineData("x^3*(a + b*arccos(c*x))/(d + k*x^2)^3", Inside)]
35+
[InlineData("x^3*(a + b*arccot(c*x))/(d + k*x^2)^3", Inside)]
36+
[InlineData("x^3*arcsin(c*x)/(d + k*x^2)^3", Inside)]
37+
[InlineData("x^5*(a + b*arctan(c*x))/(d + k*x^2)^4", Inside)]
38+
public void ByPartsAgainstTheOneTermAntiderivative(string integrand, string points)
39+
{
40+
var integral = integrand.ToEntity().Integrate("x");
41+
var text = integral.Stringize();
42+
Assert.DoesNotContain("integral(", text);
43+
Assert.True(text.Length < 20000, $"{text.Length} characters of answer for {integrand}");
44+
foreach (var (d, k) in new[] { (2.1, 0.9), (-1.7, 1.3), (1.9, -0.8) })
45+
{
46+
Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.6)
47+
.Substitute("d", d).Substitute("k", k);
48+
var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
49+
var original = Pinned(integrand.ToEntity());
50+
var compared = 0;
51+
foreach (var at in points.Split(' ').Select(point => double.Parse(point, CultureInfo.InvariantCulture)))
52+
{
53+
var want = original.Substitute("x", at).EvalNumerical();
54+
var got = derivative.Substitute("x", at).EvalNumerical();
55+
if (want.IsNaN)
56+
continue;
57+
compared++;
58+
Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
59+
< 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)),
60+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, d = {d}, k = {k}, where the integrand is {want}");
61+
}
62+
Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand} at d = {d}, k = {k}");
63+
}
64+
}
65+
}
66+
}

0 commit comments

Comments
 (0)