Skip to content

Commit 2d17cd1

Browse files
Integrate the reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one over its roots (#1905)
* Integrate the reciprocal of a binomial in a power of a sine over its roots With q = 2 for an even n and 1 for an odd one, m = n/q and r_k the m-th roots of -a/b, 1/(a + b s^n) = (1/(m a)) sum r_k/(r_k - s^q), and each term is an arctangent of tan(y) or of tan(y/2), written F(K u)/K so that only K^2 is fixed. The roots of unity are a cosine and a sine: (-1)^(2/3) is the real cube root of 1 here. Rubi's 4.1.7 and 4.2.7 1/(a + b sin(x)^5), 1/(a - b cos(x)^6) and the like, with a symbol in a or b, were declined. Part of #718 Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Integrate the reciprocal of a binomial in a power of a hyperbolic sine or cosine over its roots The same split as for the sine and cosine, each term an inverse hyperbolic tangent of tanh(y) or of tanh(y/2), the functions read back from (e^y -+ e^(-y))/2 as the parser writes them. Rubi's 6.1.7 and 6.2.7 1/(a + b sinh(x)^5), 1/(a - b cosh(x)^3) and the like were declined. Part of #718 Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Leave a hyperbolic cube to the rules before the split over roots They write 1/(a + b cosh(x)^3) through e^x in 655 characters, where the split writes 1491. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the split over a binomial's roots in the breaking changes Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --------- Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent 5b5252b commit 2d17cd1

4 files changed

Lines changed: 179 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -396,6 +396,26 @@ in Rubi's 1.1.3.4 come out shorter as well, and those that took two seconds take
396396
| `"1/(x^8*(a + b*x^3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 1,217 characters |
397397
| `"1/(x^6*(a + b*x^4)*(c + d*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 795 characters |
398398

399+
### The reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one is split over its roots
400+
401+
**Answers where there were none.** `1/(a + b sin(x)^5)` was declined, and so were twenty more of Rubi's 4.1.7, 4.2.7,
402+
6.1.7 and 6.2.7. With `s` the function, `n` from 5 to 12, `q` 2 for an even `n` and 1 for an odd one, and `r_k` the
403+
`n/q`-th roots of `-a/b`, `1/(a + b s^n)` is `(q/(n a))` times the sum of `r_k/(r_k - s^q)`, and each term is an
404+
arctangent, or an inverse hyperbolic tangent, of the tangent of the argument or of its half
405+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
406+
407+
| Input | Was (2.5.0) | Now |
408+
|---|---|---|
409+
| `"1/(a+b*sin(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,110 characters |
410+
| `"1/(a-b*sin(x)^6)".ToEntity().Integrate("x")` | `integral(...)` | 648 characters |
411+
| `"1/(a+b*sin(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 569 characters |
412+
| `"1/(a+b*cos(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,362 characters |
413+
| `"1/(a-b*cos(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 841 characters |
414+
| `"1/(a+b*sinh(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,678 characters |
415+
| `"1/(a+b*sinh(x)^6)".ToEntity().Integrate("x")` | `integral(...)` | 1,139 characters |
416+
| `"1/(a+b*cosh(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 3,607 characters |
417+
| `"1/(a-b*cosh(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 1,441 characters |
418+
399419
### A whole power of an inverse function of `c + d x^2` is integrated by its recurrence
400420

401421
**Answers where there were none.** `(a + b acos(1 + d x^2))^4` was declined, and so were twenty-five more of Rubi's

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

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25940,6 +25940,98 @@ static bool IsASumOfPowersOfX(Entity sum, Entity.Variable x)
2594025940
return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
2594125941
}
2594225942

25943+
/// <summary>
25944+
/// The reciprocal of <c>a + b s^n</c>, <c>s</c> a sine or a cosine of a linear and a symbol in
25945+
/// <c>a</c> or <c>b</c>, over the roots of the binomial: with <c>q = 2</c> for an even <c>n</c>
25946+
/// and <c>1</c> for an odd one, <c>m = n/q</c> and <c>r_k</c> the <c>m</c>-th roots of
25947+
/// <c>-a/b</c>, <c>1/(a + b s^n) = (1/(m a)) sum r_k/(r_k - s^q)</c>, and each term is an
25948+
/// arctangent of <c>tan(y)</c> or of <c>tan(y/2)</c>.
25949+
/// </summary>
25950+
/// <remarks>
25951+
/// Each term is written <c>F(K u)/K</c> with only <c>K^2</c> fixed, so its derivative is the
25952+
/// same whichever root <c>K</c> is: <c>(atan(K tan(y))/K)' = 1/(cos(y)^2 + K^2 sin(y)^2)</c>.
25953+
/// A root that makes <c>K^2</c> zero or infinite -- <c>1/(1 + sin(y)^5)</c> has <c>-1</c> -- is
25954+
/// another shape, and declines. <c>n</c> up to <c>4</c>, and numbers for <c>a</c> and
25955+
/// <c>b</c>, are answered by the rules before it. Rubi's 4.1.7 and 4.2.7
25956+
/// <c>1/(a + b sin(x)^5)</c>, <c>1/(a - b cos(x)^6)</c> and the like were declined.
25957+
/// https://github.com/asc-community/AngouriMath/issues/718
25958+
/// </remarks>
25959+
internal static Entity? SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(Entity expr, Entity.Variable x)
25960+
{
25961+
if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || numerator.ContainsNode(x))
25962+
return null;
25963+
// The hyperbolic sine and cosine are parsed as `(e^y -+ e^(-y))/2`, and read back from it.
25964+
static (Entity Argument, bool Sine, bool Hyperbolic)? Read(Entity node) => node switch
25965+
{
25966+
Sinf(var y) => (y, true, false),
25967+
Cosf(var y) => (y, false, false),
25968+
Divf(Sumf(Powf(var e1, var y), Powf(var e2, var minusY)), Number.Integer { EInteger: var two })
25969+
when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e
25970+
&& Functions.PartialFractions.IsZeroAsAValue((y + minusY).InnerSimplified) => (y, false, true),
25971+
Divf(Minusf(Powf(var e1, var y), Powf(var e2, var minusY)), Number.Integer { EInteger: var two })
25972+
when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e
25973+
&& Functions.PartialFractions.IsZeroAsAValue((y + minusY).InnerSimplified) => (y, true, true),
25974+
_ => null
25975+
};
25976+
var functions = denominator.Nodes.Where(node => node.ContainsNode(x) && Read(node) is not null).Distinct().ToList();
25977+
if (functions.Count != 1)
25978+
return null;
25979+
var function = functions[0];
25980+
var (argument, sine, hyperbolic) = Read(function)!.Value;
25981+
if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out _) || Functions.PartialFractions.IsZeroAsAValue(slope))
25982+
return null;
25983+
var s = Variable.CreateUnique(expr, "s_bin");
25984+
var inS = denominator.Substitute(function, s);
25985+
if (inS.ContainsNode(x) || !TreeAnalyzer.TryGetPolynomial(inS, s, out var read) || read.Count != 2
25986+
|| !read.ContainsKey(EInteger.Zero) || read.Keys.Max() is not { } top || !top.CanFitInInt32())
25987+
return null;
25988+
var n = top.ToInt32Unchecked();
25989+
// Past what the rules before it answer: they write `1/(a + b cosh(x)^3)` through `e^x` in
25990+
// 655 characters, where this writes 1491.
25991+
if (n < 5 || n > 12)
25992+
return null;
25993+
var a = read[EInteger.Zero];
25994+
var b = read[top];
25995+
if (!(a + b).Vars.Any() || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b))
25996+
return null;
25997+
var q = n % 2 == 0 ? 2 : 1;
25998+
var m = n / q;
25999+
var half = q == 2 ? argument : argument / 2;
26000+
var u = hyperbolic ? MathS.Hyperbolic.Tanh(half) : MathS.Tan(half);
26001+
var radius = MathS.Pow(-a / b, Number.Rational.Create(1, m));
26002+
Entity sum = Number.Integer.Zero;
26003+
for (var k = 0; k < m; k++)
26004+
{
26005+
// The root of unity as a cosine and a sine: `(-1)^(2/3)` is 1 here, the real cube
26006+
// root of 1, where it means the principal `e^(2 pi i/3)` to Rubi.
26007+
var turn = 2 * MathS.pi * k / m;
26008+
var r = (k == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn))).InnerSimplified;
26009+
// K^2, the factor in front, and the argument of the arctangent over K.
26010+
var (square, front, shift) = (q, sine, hyperbolic) switch
26011+
{
26012+
(2, true, false) => ((r - 1) / r, (Entity)Number.Integer.One, (Entity)Number.Integer.Zero),
26013+
(2, false, _) => (r / (r - 1), r / (r - 1), Number.Integer.Zero),
26014+
(1, true, false) => ((MathS.Sqr(r) - 1) / MathS.Sqr(r), 2, -1 / r),
26015+
(2, true, true) => ((r + 1) / r, Number.Integer.One, Number.Integer.Zero),
26016+
(1, true, true) => (1 + 1 / MathS.Sqr(r), 2, 1 / r),
26017+
_ => ((r - 1) / (r + 1), 2 * r / (r + 1), Number.Integer.Zero),
26018+
};
26019+
square = square.InnerSimplified;
26020+
if (Functions.PartialFractions.IsZeroAsAValue(square) || Functions.PartialFractions.IsZeroAsAValue((1 / square).InnerSimplified))
26021+
return null;
26022+
var root = MathS.Sqrt(square);
26023+
// An even power under tan(y), where K multiplies u; an odd one under tan(y/2), where it divides.
26024+
// The hyperbolic ones are the same with the inverse hyperbolic tangent.
26025+
Entity Inverse(Entity z) => hyperbolic ? MathS.Hyperbolic.Artanh(z) : MathS.Arctan(z);
26026+
sum += q == 2 ? front * Inverse(root * u) / root : front * Inverse((u + shift) / root) / root;
26027+
}
26028+
var answer = (numerator / (m * a * slope) * sum).InnerSimplified;
26029+
// Refused where it differs, not where it cannot be told: the intervals do not settle
26030+
// the arctangents of complex arguments at every sampled point, and `1/(a + b sin(x)^5)`,
26031+
// right at every point checked exactly, was declined for too few comparisons.
26032+
return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer;
26033+
}
26034+
2594326035
/// <summary>
2594426036
/// A rational function of a sine, or of a cosine, over a quadratic in it with a symbol
2594526037
/// among the coefficients: <c>sin(x)/(a + b sin(x) + c sin(x)^2)</c>, and the same with

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

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -968,6 +968,8 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
968968
// among the coefficients, split over the quadratic's two roots: under the
969969
// half-angle it is a quartic in t with a symbol in every coefficient.
970970
if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOfASineOverASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer;
971+
// And over a binomial in a higher power of it, over the binomial's roots.
972+
if ((answer = IndefiniteIntegralSolver.SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer;
971973
// Fractional powers of `a ± a sin(y)` beside anything rational in the sine and
972974
// cosine, by the half angle at which they are squares: `1 + sin(y)` is `2 sin(u)^2`. Before
973975
// the substitution search, which spent twenty seconds on the radicals of the sine.
Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
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+
/// The reciprocal of <c>a + b s^n</c>, <c>s</c> a sine or a cosine and <c>n</c> past 4, or a
16+
/// hyperbolic sine or cosine, over the roots of the binomial: Rubi's 4.1.7, 4.2.7, 6.1.7 and 6.2.7. Pinned with <c>-a/b</c> of either sign, so the
17+
/// real root among them is positive or negative.
18+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
19+
/// </summary>
20+
[Trait("Area", "Calculus")]
21+
public sealed class AReciprocalOfABinomialInAPowerOfASineIntegralTest
22+
{
23+
private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
24+
{
25+
var integral = integrand.ToEntity().Integrate("x");
26+
Assert.DoesNotContain("integral(", integral.Stringize());
27+
28+
var pinned = integral.Substitute("C", 0);
29+
Entity original = integrand.ToEntity();
30+
foreach (var (name, value) in pins)
31+
{
32+
pinned = pinned.Substitute(name, value);
33+
original = original.Substitute(name, value);
34+
}
35+
var derivative = pinned.Differentiate("x");
36+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
37+
{
38+
var got = derivative.Substitute("x", at).EvalNumerical();
39+
var want = original.Substitute("x", at).EvalNumerical();
40+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
41+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
42+
Assert.True(difference / scale < 1e-8,
43+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
44+
}
45+
}
46+
47+
[Theory]
48+
[InlineData("1/(a+b*sin(x)^5)")]
49+
[InlineData("1/(a-b*sin(x)^6)")]
50+
[InlineData("1/(a+b*sin(x)^8)")]
51+
[InlineData("1/(a-b*cos(x)^5)")]
52+
[InlineData("1/(a+b*cos(x)^6)")]
53+
[InlineData("3/(a+b*sin(2*x+1)^7)")]
54+
[InlineData("1/(a+b*cosh(x)^7)")]
55+
[InlineData("1/(a-b*cosh(x)^5)")]
56+
[InlineData("1/(a+b*sinh(x)^6)")]
57+
[InlineData("1/(a+b*sinh(2*x+1)^5)")]
58+
[InlineData("1/(a-b*cosh(x)^8)")]
59+
public void IsAnsweredAndDifferentiatesBack(string integrand)
60+
{
61+
DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25));
62+
DifferentiatesBack(integrand, ("a", -0.75), ("b", 0.4));
63+
}
64+
}
65+
}

0 commit comments

Comments
 (0)