Skip to content

Commit 607fd57

Browse files
Integrate a sum of powers of an inverse or error function by parts down the powers (#1972)
A sum of terms, each a power of one function K whose derivative is free of it (an inverse trigonometric function or an error function) times a cofactor free of it, is now integrated by parts down the powers. The rule uses ∫ A K^p = F K^p - ∫ p F K^(p - 1) K', with F = ∫ A for the cofactor of the top power, and the last integral joins the terms a power below. Where every power cancels that way, the parts are the antiderivative: x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2) is (-x^3/(2 a atan(a x)^2))', in one step. Rubi's 5.1.4:534, 5.3.4:768 and 8.1:122/251/382 were declined or 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 56648c3 commit 607fd57

4 files changed

Lines changed: 219 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6434,6 +6434,29 @@ no integrand without one changes: Rubi's independent suites and a sample of its
64346434
| `"e^(b*x)*Ei(b*x)/x".Integrate("x")` | `Ei * e ^ (b * x) + C`, with `Ei` a variable | `Ei(b * x) ^ 2 / 2 + C` |
64356435
| `"Si(b*x)*sin(b*x)/x".Integrate("x")` | `Si * -cos(b * x) + C`, with `Si` a variable | `Si(b * x) ^ 2 / 2 + C` |
64366436

6437+
### A sum of powers of an inverse trigonometric or error function is integrated by parts down the powers
6438+
6439+
**Answers where there were none.** `x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2)`,
6440+
`x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x)))` and `erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x)` with its
6441+
`erfc` and `erfi` twins, Rubi's 5.1.4, 5.3.4 and 8.1, were declined: no term is elementary alone, and the split and
6442+
the parts each read one term at a time. A sum of terms each a power of one function `K`, whose derivative is free of
6443+
it, times a cofactor free of it is now integrated by parts from the top power,
6444+
`integral A K^p = F K^p - integral p F K^(p - 1) K'` with `F` the antiderivative of `A`, and the last integral joins
6445+
the terms a power below. Where every power cancels so, the
6446+
parts are the answer, `(-x^3/(2 a atan(a x)^2))'` for the first; where what is left is free of `K` it is asked. The
6447+
powers a whole number apart and at most five apart, a cancelling group decided at sampled points, and the answer
6448+
differentiated back. Three sums of the same kind the unreleased master answered come back shorter
6449+
([#233](https://github.com/asc-community/AngouriMath/issues/233)).
6450+
6451+
| Input | Was (2.5.0) | Now |
6452+
|---|---|---|
6453+
| `"x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)".ToEntity().Integrate("x")` | `integral(...)` | `-x^3/(2 a arctan(a x)^2)`, 51 characters |
6454+
| `"x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))".ToEntity().Integrate("x")` | `integral(...)` | 99 characters |
6455+
| `"erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | 219 characters |
6456+
| `"erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `integral(...)`, reading `erfi` as a symbol | 131 characters |
6457+
| `"erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`; 352 characters on the unreleased master | `x erf(x)^2`, 18 characters |
6458+
| `"2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)".ToEntity().Integrate("x")` | `integral(...)`; 121 characters on the unreleased master | `x^2/arccos(x)`, 33 characters |
6459+
64376460
### A square of a special function, and one beside a power of `x` and its elementary derivative, are integrated by parts
64386461

64396462
`x Si(b x) sin(b x)` was left unevaluated where `Si(b x) sin(b x)` was not: by parts it is

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

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15991,6 +15991,137 @@ static bool More(Entity difference)
1599115991
return changed ? rebuilt / bottom : null;
1599215992
}
1599315993

15994+
/// <summary>
15995+
/// A sum of terms each a power of one function <c>K</c> whose derivative is free of it -- an
15996+
/// inverse trigonometric function or an error function -- times a cofactor free of it,
15997+
/// integrated by parts down the powers: <c>∫ A K^p = F K^p - ∫ p F K^(p - 1) K'</c> for
15998+
/// <c>F = ∫ A</c>, the cofactor of the top power, and the last integral joins the terms a
15999+
/// power below. Where every power cancels so, the parts are the antiderivative:
16000+
/// <c>x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2)</c> is
16001+
/// <c>(-x^3/(2 a atan(a x)^2))'</c>, in one step, and
16002+
/// <c>x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x)))</c> in two. Where what is
16003+
/// left is free of <c>K</c>, it is asked: <c>erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x)</c> leaves
16004+
/// <c>b e^(-2 b^2 x^2)/(sqrt(pi) x^2)</c>.
16005+
/// </summary>
16006+
/// <remarks>
16007+
/// <para>
16008+
/// The terms are not elementary apart -- <c>∫ x/((1 - x^2) sqrt(asin(x)))</c> is
16009+
/// <c>∫ tan(t)/sqrt(t) dt</c> -- so the split declined each, and by parts reads one term
16010+
/// at a time. Rubi's 5.1.4, 5.3.4 and 8.1 write such sums to test that the parts meet.
16011+
/// </para>
16012+
/// <para>
16013+
/// One such function in the sum, every term a rational power of it times a cofactor free
16014+
/// of it, the powers a whole number apart, at most five apart, and each step's
16015+
/// <c>F</c> free of <c>K</c>. A group that cancels is decided at sampled points, and the
16016+
/// answer is checked by differentiating it back. A remainder free of <c>K</c> is asked
16017+
/// only for a top power of two terms or more with a cofactor among them that is not
16018+
/// algebraic, which is where the cofactors integrate together and not apart: a sum of
16019+
/// algebraic cofactors, <c>a + b asin(c x)</c>, is what the split answers term by term,
16020+
/// and this rule would only respell it. At the question asked only, and not inside its
16021+
/// own questions, which are free of <c>K</c> anyway.
16022+
/// </para>
16023+
/// https://github.com/asc-community/AngouriMath/issues/233
16024+
/// </remarks>
16025+
internal static Entity? SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(Entity expr, Entity.Variable x, bool integrateByParts)
16026+
{
16027+
if (!Integration.AnsweringTheQuestionAsked || expr is not (Sumf or Minusf) || descendingThePowers)
16028+
return null;
16029+
Entity? kernel = null;
16030+
foreach (var node in expr.Nodes)
16031+
if (node is Arcsinf or Arccosf or Arctanf or Arccotanf or Erff or Erfcf or Erfif && node.ContainsNode(x))
16032+
{
16033+
if (kernel is null)
16034+
kernel = node;
16035+
else if (node != kernel)
16036+
return null;
16037+
}
16038+
if (kernel is null)
16039+
return null;
16040+
var slope = kernel.Differentiate(x).InnerSimplified;
16041+
if (slope.ContainsNode(kernel))
16042+
return null;
16043+
Number.Rational? PowerOfTheKernel(Entity factor)
16044+
=> factor == kernel ? Number.Integer.One
16045+
: factor is Powf(var @base, Number.Rational power) && PowerOfTheKernel(@base) is { } inner ? inner * power
16046+
: null;
16047+
var terms = new List<(Number.Rational Power, Entity Cofactor)>();
16048+
foreach (var term in Sumf.LinearChildren(expr))
16049+
{
16050+
Number.Rational power = Number.Integer.Zero;
16051+
Entity cofactor = Number.Integer.One;
16052+
foreach (var factor in Mulf.LinearChildren(term))
16053+
if (!factor.ContainsNode(kernel))
16054+
cofactor *= factor;
16055+
else if (PowerOfTheKernel(factor) is { } read)
16056+
power += read;
16057+
else
16058+
return null;
16059+
terms.Add((power, cofactor));
16060+
}
16061+
var top = terms.Select(term => term.Power).Aggregate((l, r) => Above(r, l) ? r : l);
16062+
var lowest = terms.Select(term => term.Power).Aggregate((l, r) => Above(l, r) ? r : l);
16063+
if (terms.Any(term => !(top - term.Power).ERational.Denominator.Equals(EInteger.One))
16064+
|| Above(top - lowest, Number.Integer.Create(5)))
16065+
return null;
16066+
var groups = new List<(Number.Rational Power, Entity Cofactor)>();
16067+
foreach (var (power, cofactor) in terms)
16068+
Join(groups, power, cofactor);
16069+
var atTheTop = terms.Where(term => Same(term.Power, top)).ToList();
16070+
var mayAskWhatIsLeft = atTheTop.Count >= 2 && atTheTop.Any(term => !IsAlgebraicIn(term.Cofactor, x));
16071+
if (groups.Count < 2 && !mayAskWhatIsLeft)
16072+
return null;
16073+
descendingThePowers = true;
16074+
try
16075+
{
16076+
Entity parts = Number.Integer.Zero;
16077+
for (var step = 0; step < 8; step++)
16078+
{
16079+
groups.RemoveAll(group => Functions.PartialFractions.HoldsAtSampledPoints(group.Cofactor, Number.Integer.Zero, x));
16080+
if (groups.Count == 0)
16081+
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(parts, expr, x) ? parts : null;
16082+
var (power, cofactor) = groups.Aggregate((l, r) => Above(r.Power, l.Power) ? r : l);
16083+
if (power.ERational.IsZero)
16084+
{
16085+
// Free of K: asked where the top power's cofactors integrate together only.
16086+
if (groups.Count > 1 || !mayAskWhatIsLeft
16087+
|| Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } rest)
16088+
return null;
16089+
var answer = parts + rest;
16090+
return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
16091+
}
16092+
// Below zero the powers never reach a remainder free of K, so below the lowest
16093+
// written one nothing is left to cancel.
16094+
if (power.ERational.Sign < 0 && Above(lowest, power)
16095+
|| Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } integrated
16096+
|| integrated.ContainsNode(kernel))
16097+
return null;
16098+
parts += integrated * (Same(power, Number.Integer.One) ? kernel : MathS.Pow(kernel, power));
16099+
groups.RemoveAll(group => Same(group.Power, power));
16100+
Join(groups, power - Number.Integer.One, -power * integrated * slope);
16101+
for (var i = 0; i < groups.Count; i++)
16102+
groups[i] = (groups[i].Power, groups[i].Cofactor.InnerSimplified);
16103+
}
16104+
return null;
16105+
}
16106+
finally
16107+
{
16108+
descendingThePowers = false;
16109+
}
16110+
static bool Same(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) == 0;
16111+
static bool Above(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) > 0;
16112+
static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.Rational power, Entity cofactor)
16113+
{
16114+
var at = groups.FindIndex(group => Same(group.Power, power));
16115+
if (at < 0)
16116+
groups.Add((power, cofactor));
16117+
else
16118+
groups[at] = (power, groups[at].Cofactor + cofactor);
16119+
}
16120+
}
16121+
16122+
/// <summary>Set while <see cref="SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers"/> asks its questions.</summary>
16123+
[System.ThreadStatic] private static bool descendingThePowers;
16124+
1599416125
/// <summary>
1599516126
/// A product of powers times a sum that is the derivative of the product with some of
1599616127
/// the powers raised by one: <c>e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))</c> is

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

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -940,6 +940,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
940940
// recurrence that raises the whole power to -1: closed, where the substitution for the
941941
// root makes a rational function over a power of a symbolic quadratic.
942942
if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfALinearBesideAPowerOfAnotherByTheRecurrence(expr, x, integrateByParts)) is { }) return answer;
943+
// A sum of powers of one inverse trigonometric or error function beside cofactors free
944+
// of it, by parts down the powers from the top: `x^3/((1 + a^2 x^2) atan(a x)^3) -
945+
// 3 x^2/(2 a atan(a x)^2)` is the derivative of one product, and neither term is
946+
// elementary on its own. Before the split, which loses it.
947+
if ((answer = IndefiniteIntegralSolver.SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(expr, x, integrateByParts)) is { }) return answer;
943948
// A product of powers times a sum that is the derivative of the product with some
944949
// of the powers raised by one -- `e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))` -- read off
945950
// the sum, for symbolic exponents too. Before the split, which loses it.
Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
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 sum of powers of one inverse trigonometric or error function beside cofactors free of
16+
/// it, by parts down the powers from the top. Its terms are not elementary apart, so the
17+
/// split declined each, and the sum was declined or ran out of time. Rubi's 5.1.4, 5.3.4
18+
/// and 8.1.
19+
/// <a href="https://github.com/asc-community/AngouriMath/issues/233">#233</a>
20+
/// </summary>
21+
[Trait("Area", "Calculus")]
22+
public sealed class ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest
23+
{
24+
[Theory]
25+
// The parts cancel every power below the top: one step, and two.
26+
[InlineData("x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)")]
27+
[InlineData("x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))")]
28+
[InlineData("1/((1-x^2)*arcsin(x)^2)-x/((1-x^2)^(3/2)*arcsin(x))")]
29+
[InlineData("2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)")]
30+
[InlineData("erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)")]
31+
// What is left is free of the function, and the cofactors integrate together only.
32+
[InlineData("erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)")]
33+
[InlineData("erfc(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfc(b*x)/(e^(b^2*x^2)*x)")]
34+
[InlineData("erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)")]
35+
public void ByPartsDownThePowers(string integrand)
36+
{
37+
var integral = integrand.ToEntity().Integrate("x");
38+
var text = integral.Stringize();
39+
Assert.DoesNotContain("integral(", text);
40+
Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}");
41+
Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7);
42+
var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
43+
var original = Pinned(integrand.ToEntity());
44+
var compared = 0;
45+
// Inside (0, 1), where the arcsine and the arccosine are real and positive.
46+
foreach (var at in new[] { 0.15, 0.3, 0.45, 0.6, 0.75, 0.9 })
47+
{
48+
var want = original.Substitute("x", at).EvalNumerical();
49+
var got = derivative.Substitute("x", at).EvalNumerical();
50+
if (want.IsNaN)
51+
continue;
52+
compared++;
53+
Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
54+
< 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)),
55+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
56+
}
57+
Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}");
58+
}
59+
}
60+
}

0 commit comments

Comments
 (0)