Skip to content

Commit 0d72aa4

Browse files
Integrate a product of powers of a sine and a cosine beside a non-polynomial by linearising it (#1924)
* Integrate a product of powers of a sine and a cosine beside a non-polynomial by linearising it Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Read the dictionaries' pairs without deconstruction, which netstandard2.0 lacks Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura * Record the linearised powers of a sine and a cosine in 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 6796bc8 commit 0d72aa4

4 files changed

Lines changed: 171 additions & 0 deletions

File tree

‎BREAKING-CHANGES.md‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -281,6 +281,19 @@ divided by Smith's method now, which scales by the larger part and squares nothi
281281
|---|---|---|
282282
| `(1.98e119 + i)/(3.39e125 + i)` with the downcasting off | `0` | `5.8407…e-7` |
283283

284+
### A product of powers of a sine and a cosine beside a non-polynomial is linearised
285+
286+
**Answers where there were none.** `x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))` was declined, and its square ran out of
287+
time (Rubi's 5.1.4). Under `x = sin(u)/c` they become `sin(u)^p cos(u)^q` beside a power of `a + b u`. A product of
288+
powers of a sine and a cosine beside such a factor is now written as a sum of sines and cosines of multiples of the
289+
angle, each of which the integrator answers
290+
([#718](https://github.com/asc-community/AngouriMath/issues/718)).
291+
292+
| Input | Was (2.5.0) | Now |
293+
|---|---|---|
294+
| `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))".ToEntity().Integrate("x")` | `integral(...)` | 508 characters |
295+
| `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))^2".ToEntity().Integrate("x")` | `integral(...)` | 1,470 characters |
296+
284297
### A piecewise is defined where the case it takes is
285298

286299
**Domain conditions that were false.** The domain condition of a piecewise asked every case's predicate

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

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26347,6 +26347,104 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e
2634726347
return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
2634826348
}
2634926349

26350+
/// <summary>
26351+
/// A product of whole powers of a sine and a cosine of one linear times something else in x
26352+
/// that is no polynomial, <c>sin(u)^3 cos(u)^6/(a + b u)</c>, linearised: the product is a sum
26353+
/// of <c>cos(k u)</c> and <c>sin(k u)</c> by the product-to-sum formulas, and each term beside
26354+
/// the rest its own question.
26355+
/// </summary>
26356+
/// <remarks>
26357+
/// <c>sin(3 u)/(a + b u)</c> is answered in sines and cosine integrals, and the product of
26358+
/// powers was not: Rubi's 5.1.4 <c>x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))</c>, which under
26359+
/// <c>u = asin(c x)</c> is that product, was declined. Last but one, so that everything
26360+
/// answering the product as written goes first.
26361+
/// https://github.com/asc-community/AngouriMath/issues/718
26362+
/// </remarks>
26363+
internal static Entity? SolveAProductOfPowersOfASineAndACosineByLinearising(Entity expr, Entity.Variable x)
26364+
{
26365+
if (!Integration.AnsweringTheQuestionAsked)
26366+
return null;
26367+
static IEnumerable<(Entity Factor, bool Below)> Factors(Entity e, bool below) => e switch
26368+
{
26369+
Mulf(var l, var r) => Factors(l, below).Concat(Factors(r, below)),
26370+
Divf(var l, var r) => Factors(l, below).Concat(Factors(r, !below)),
26371+
_ => new[] { (e, below) }
26372+
};
26373+
Entity? argument = null;
26374+
int sines = 0, cosines = 0;
26375+
Entity rest = Number.Integer.One;
26376+
foreach (var (factor, below) in Factors(expr, false))
26377+
{
26378+
var (function, power) = factor is Powf(var b0, Number.Integer { EInteger: var p0 }) && b0 is Sinf or Cosf && p0.Sign > 0 && p0.CanFitInInt32()
26379+
? (b0, p0.ToInt32Unchecked()) : (factor, 1);
26380+
if (!below && function is Sinf or Cosf && function.ContainsNode(x))
26381+
{
26382+
var y0 = function is Sinf(var ys) ? ys : ((Cosf)function).Argument;
26383+
if (argument is not null && argument != y0)
26384+
return null;
26385+
argument = y0;
26386+
if (function is Sinf) sines += power; else cosines += power;
26387+
}
26388+
else
26389+
rest = below ? rest / factor : rest * factor;
26390+
}
26391+
if (argument is null || sines + cosines < 2 || sines + cosines > 24
26392+
|| !TreeAnalyzer.TryGetPolyLinear(argument, x, out _, out _)
26393+
|| !rest.ContainsNode(x) || rest.Nodes.Any(node => node is Sinf or Cosf or Tanf or Cotanf or Secantf or Cosecantf)
26394+
|| TreeAnalyzer.TryGetPolynomial(rest, x, out var asPolynomial) && asPolynomial.Keys.All(key => key.Sign >= 0))
26395+
return null;
26396+
// cos(k y) and sin(k y) coefficients, multiplied by sin y and cos y one at a time.
26397+
var cosine = new Dictionary<int, ERational> { [0] = ERational.One };
26398+
var sine = new Dictionary<int, ERational>();
26399+
static void Add(Dictionary<int, ERational> into, int k, ERational value)
26400+
=> into[k] = (into.TryGetValue(k, out var old) ? old + value : value).ToLowestTerms();
26401+
var half = ERational.Create(1, 2);
26402+
for (var step = 0; step < sines + cosines; step++)
26403+
{
26404+
var bySine = step < sines;
26405+
var (nextCosine, nextSine) = (new Dictionary<int, ERational>(), new Dictionary<int, ERational>());
26406+
foreach (var (k, value) in cosine.Select(pair => (pair.Key, pair.Value)))
26407+
if (bySine)
26408+
{ // cos(k y) sin y = (sin((k + 1) y) - sin((k - 1) y))/2
26409+
Add(nextSine, k + 1, value * half);
26410+
if (k - 1 > 0) Add(nextSine, k - 1, -value * half);
26411+
else if (k - 1 < 0) Add(nextSine, 1 - k, value * half);
26412+
}
26413+
else
26414+
{ // cos(k y) cos y = (cos((k + 1) y) + cos((k - 1) y))/2
26415+
Add(nextCosine, k + 1, value * half);
26416+
Add(nextCosine, System.Math.Abs(k - 1), value * half);
26417+
}
26418+
foreach (var (k, value) in sine.Select(pair => (pair.Key, pair.Value)))
26419+
if (bySine)
26420+
{ // sin(k y) sin y = (cos((k - 1) y) - cos((k + 1) y))/2
26421+
Add(nextCosine, System.Math.Abs(k - 1), value * half);
26422+
Add(nextCosine, k + 1, -value * half);
26423+
}
26424+
else
26425+
{ // sin(k y) cos y = (sin((k + 1) y) + sin((k - 1) y))/2
26426+
Add(nextSine, k + 1, value * half);
26427+
if (k - 1 > 0) Add(nextSine, k - 1, value * half);
26428+
else if (k - 1 < 0) Add(nextSine, 1 - k, -value * half);
26429+
}
26430+
(cosine, sine) = (nextCosine, nextSine);
26431+
}
26432+
Entity answer = Number.Integer.Zero;
26433+
foreach (var (terms, isSine) in new[] { (cosine, false), (sine, true) })
26434+
foreach (var (k, value) in terms.OrderBy(term => term.Key).Select(pair => (pair.Key, pair.Value)))
26435+
{
26436+
if (value.IsZero)
26437+
continue;
26438+
var multiple = k == 1 ? argument : Number.Integer.Create(k) * argument;
26439+
Entity piece = k == 0 ? rest : (isSine ? MathS.Sin(multiple) : MathS.Cos(multiple)) * rest;
26440+
if (Integration.ComputeAsAQuestionOfItsOwn(piece, x, integrateByParts: true) is not { } integrated)
26441+
return null;
26442+
answer += Number.Rational.Create(value.Numerator, value.Denominator) * integrated;
26443+
}
26444+
answer = answer.InnerSimplified;
26445+
return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
26446+
}
26447+
2635026448
/// <summary>
2635126449
/// A rational function of a sine, or of a cosine, over a quadratic in it with a symbol
2635226450
/// 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: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1364,6 +1364,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
13641364
if ((answer = IndefiniteIntegralSolver.SolvePowersOfOneLinearAsOnePower(expr, x, integrateByParts)) is { }) return answer;
13651365
// And beside a polynomial too, where everything before has declined: `x ln(e (f (a + b x)^p (c + d x)^q)^r)`.
13661366
if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionTimesALogarithmOfPowersOfLinears(expr, x, integrateByParts, besideAPolynomialToo: true)) is { }) return answer;
1367+
// A product of powers of a sine and a cosine beside something that is no polynomial,
1368+
// linearised, where everything answering it as written has declined: `sin(u)^3 cos(u)^6/(a + b u)`.
1369+
if ((answer = IndefiniteIntegralSolver.SolveAProductOfPowersOfASineAndACosineByLinearising(expr, x)) is { }) return answer;
13671370
// Last of the rewrites, because it is the only one that fires on an integrand nothing
13681371
// is wrong with -- it clears a parameter rather than a shape -- so everything that
13691372
// can answer the problem as written gets to try first.
Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
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 product of powers of a sine and a cosine beside something that is no polynomial, linearised:
16+
/// Rubi's 5.1.4 <c>x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))</c>, which was declined.
17+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
18+
/// </summary>
19+
[Trait("Area", "Calculus")]
20+
public sealed class AProductOfPowersOfASineAndACosineLinearisedIntegralTest
21+
{
22+
private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
23+
{
24+
var integral = integrand.ToEntity().Integrate("x");
25+
Assert.DoesNotContain("integral(", integral.Stringize());
26+
27+
var pinned = integral.Substitute("C", 0);
28+
Entity original = integrand.ToEntity();
29+
foreach (var (name, value) in pins)
30+
{
31+
pinned = pinned.Substitute(name, value);
32+
original = original.Substitute(name, value);
33+
}
34+
var derivative = pinned.Differentiate("x");
35+
foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
36+
{
37+
var got = derivative.Substitute("x", at).EvalNumerical();
38+
var want = original.Substitute("x", at).EvalNumerical();
39+
var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
40+
var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
41+
Assert.True(difference / scale < 1e-8,
42+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
43+
}
44+
}
45+
46+
[Theory]
47+
[InlineData("x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))")]
48+
[InlineData("x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))^2")]
49+
[InlineData("sin(3+2*x)^3*cos(3+2*x)^6/(a+b*x)")]
50+
public void IsAnsweredAndDifferentiatesBack(string integrand)
51+
{
52+
// |c x| below 1 at every sampled x, where the inverse sine is real.
53+
DifferentiatesBack(integrand, ("a", 1.3), ("b", 0.7), ("c", 0.3));
54+
DifferentiatesBack(integrand, ("a", -2), ("b", 1.5), ("c", -0.2));
55+
}
56+
}
57+
}

0 commit comments

Comments
 (0)