From 9b66e6dbdbed655714feb67cb9cf24521118a52a Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Fri, 23 Oct 2026 05:26:40 +0000 Subject: [PATCH] Integrate a power beside a sine or a cosine as incomplete gamma functions A power of x or of a linear that is not a polynomial, beside sin(u)^k cos(u)^j, is integrated through Euler's formula: sin(u)^k cos(u)^j is (2 i)^(-k) 2^(-j) sum C(k, r) (-1)^(k - r) C(j, t) e^(i (2 r - k + 2 t - j) u), each e^(i q u) beside the power is a gamma function of the gamma rule for a power beside an exponential, and the term with q = 0 is the power's own integral. Where any term is declined -- a whole power, whose integral is elementary and answered by parts, or a half-whole order, which is Fresnel's -- so is the integrand. The gamma functions' arguments are imaginary on the real line, off the cut. Rubi's suite on the harness that reads GAMMA, against master 23c8eca1: 48,700 -> 48,850 of 49,052, 0 wrong, no solved row lost, no row to a timeout, no answer's allocation grown. The gains are 57 in 4.1.12, 26 in 4.1.10, 25 in 4.2.10, 18 in 4.2.12, 12 in 5.1.4, 9 in 4.7.3, 2 in 8.4 and one Moses problem. Four definite integrals agree with quadrature to 25 digits, real. Part of #2075. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 5 +- .../Integration/IntegralPatterns.cs | 118 ++++++++++++++++++ ...BesideAnExponentialOfAPowerIntegralTest.cs | 15 +++ 3 files changed, 137 insertions(+), 1 deletion(-) diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 7ecac3822..9dca6c2ea 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -148,7 +148,8 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa - **Integrals in it.** `x^m F^(a + b x^n)` and `(c + d x)^m F^(a + b x)` are `-F^a x^(m + 1) z^(-s) gamma(s, z)/n`, with `z = -b ln(F) x^n` and `s = (m + 1)/n`, where `s` is neither whole nor half a whole number, or is past 30, as `x^31 e^(-x)` is, whose sum of 32 terms is `-gamma(32, x)` now. A definite integral of such an answer that is - undefined at a bound is left as written. + undefined at a bound is left as written. And beside a power of a sine and a cosine, which Euler's formula writes as + those exponentials, the gamma functions of an imaginary argument. ([#2075](https://github.com/asc-community/AngouriMath/issues/2075)) | Input | Was (2.6.0) | Now | @@ -164,6 +165,8 @@ as `a - b = 0`, as a statement's own solve reads it: `"x + 1 = 3".Solve("x")` wa | `"x^(1/3)*gamma(2, a + b*x)".ToEntity().Integrate("x")` | a `FunctionArgumentCountException` | `-b * e ^ (-a) * x ^ (7/3) * (b * x) ^ (-7/3) * gamma(7/3, b * x) + -(1 + a) * e ^ (-a) * x ^ (4/3) * (b * x) ^ (-4/3) * gamma(4/3, b * x) + C` | | `"(c + d*x)^m*e^(a + b*x)".ToEntity().Integrate("x")` | `integral((c + d * x) ^ m * e ^ (a + b * x), x)` | `e ^ (a - b * c / d) * (c + d * x) ^ m * (-b * (c + d * x) / d) ^ (-m) * gamma(m + 1, -b * (c + d * x) / d) / b + C` | | `"x^31*e^(-x)".ToEntity().Integrate("x")` | a sum of 32 terms, 6,172 characters | `-gamma(32, x) + C` | +| `"x^m*sin(b*x)".ToEntity().Integrate("x")` | `integral(x ^ m * sin(b * x), x)` | `gamma(m + 1, -i b x)` and `gamma(m + 1, i b x)`, each beside `x ^ (m + 1) (∓i b x) ^ (-(m + 1))` | +| `"sin(c + d*x^3)".ToEntity().Integrate("x")` | `integral(sin(c + d * x ^ 3), x)` | `gamma(1/3, -i d x ^ 3)` and `gamma(1/3, i d x ^ 3)`, beside `e ^ (±i c)` | | `"integral(e^(x^3), x, 0, 1)".ToEntity().Simplify()` | `integral(e ^ x ^ 3, x, 0, 1)` | `(gamma(1/3, -1) - (-2/3)!) / 3`, which is `1.3419044179774…` | ### Numerical equality is to half the evaluation's digits diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index 52254ec9d..e9165010c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -497,6 +497,120 @@ private static bool TakeMomentFactor(Entity factor, bool inverted, Entity.Variab * MathS.Pow(under, -exponentOfPower) * MathS.IncompleteGamma((exponentOfPower + 1).InnerSimplified, under) / (logarithm * rate); } + /// + /// int P sin(u)^k cos(u)^j dx for a power P of x, or of a linear, which is not a polynomial -- or + /// none, beside a function of a power of x, sin(a + b x^3): by + /// Euler's formula sin(u)^k cos(u)^j is a sum of e^(i q u), each beside P a gamma function of + /// , and the term with q = 0 is P's own integral. Where any + /// term is declined -- a whole power of x, whose integral is elementary and answered by parts -- so is the + /// integrand. The gamma functions' arguments are imaginary on the real line, off the cut. Rubi's 4.1.10 and 4.1.12. + /// https://github.com/asc-community/AngouriMath/issues/2075 + /// + private static Entity? APowerBesideATrigonometricFunction(Entity expr, Entity.Variable x) + { + // Asked of every integrand that reaches the table, so a type test comes first, and nothing here captures a + // local before it has found a sine or a cosine. + if (expr is not (Entity.Mulf or Entity.Divf or Entity.Sinf or Entity.Cosf or Entity.Powf(Entity.Sinf or Entity.Cosf, _))) + return null; + var (numerator, denominator) = expr is Entity.Divf(var n, var d) ? (n, d) : (expr, (Entity)Entity.Number.Integer.One); + Entity coefficient = Entity.Number.Integer.One; + Entity? argument = null; + var sines = 0; + var cosines = 0; + Entity? above = null; + Entity? below = null; + foreach (var factor in Entity.Mulf.LinearChildren(numerator)) + { + if (!factor.ContainsNode(x)) + coefficient *= factor; + else if (TheTrigonometricPower(factor) is var (isSine, trigonometricArgument, times) + && (argument is null || argument == trigonometricArgument)) + { + argument = trigonometricArgument; + if (isSine) sines += times; else cosines += times; + } + else + above = above is null ? factor : above * factor; + } + if (argument is null || sines + cosines > 8) + return null; + foreach (var factor in Entity.Mulf.LinearChildren(denominator)) + if (!factor.ContainsNode(x)) + coefficient /= factor; + else + below = below is null ? factor : below * factor; + + // sin(u)^k cos(u)^j = (2 i)^(-k) 2^(-j) sum C(k, r) (-1)^(k - r) C(j, t) e^(i (2 r - k + 2 t - j) u). + var weights = new SortedDictionary(); + for (var r = 0; r <= sines; r++) + for (var t = 0; t <= cosines; t++) + { + var q = 2 * r - sines + 2 * t - cosines; + var weight = Binomial(sines, r).EInteger.Multiply(Binomial(cosines, t).EInteger); + if ((sines - r) % 2 != 0) + weight = weight.Negate(); + weights[q] = weights.TryGetValue(q, out var already) ? already.Add(weight) : weight; + } + Entity sum = Entity.Number.Integer.Zero; + // The highest frequency first: where the power is whole, the gamma rule declines it at once. + foreach (var weight in weights.Reverse()) + { + if (weight.Value.IsZero) + continue; + Entity? term; + if (weight.Key == 0) + term = TheIntegralOfAPower(above, below, x); + else + { + var exponential = MathS.Pow(MathS.e, (MathS.i * weight.Key * argument).Expand()); + var integrand = above is null ? (below is null ? exponential : exponential / below) : below is null ? above * exponential : above * exponential / below; + term = APowerBesideAnExponentialOfAPower(integrand, x); + } + if (term is null) + return null; + sum += Entity.Number.Integer.Create(weight.Value) * term; + } + return coefficient / (MathS.Pow(2 * MathS.i, sines) * MathS.Pow(2, cosines)) * sum; + } + + /// as sin(u)^k or cos(u)^k for a whole k >= 1. + private static (bool IsSine, Entity Argument, int Times)? TheTrigonometricPower(Entity factor) => factor switch + { + Entity.Sinf(var u) => (true, u, 1), + Entity.Cosf(var u) => (false, u, 1), + Entity.Powf(Entity.Sinf(var u), Entity.Number.Integer { EInteger: var k }) when k.Sign > 0 && k.CompareTo(EInteger.FromInt32(8)) <= 0 + => (true, u, k.ToInt32Checked()), + Entity.Powf(Entity.Cosf(var u), Entity.Number.Integer { EInteger: var k }) when k.Sign > 0 && k.CompareTo(EInteger.FromInt32(8)) <= 0 + => (false, u, k.ToInt32Checked()), + _ => null + }; + + /// + /// The integral of a power (d x)^m above or below the bar, or of a power (c + d x)^m, for any m but -1, + /// or . + /// + private static Entity? TheIntegralOfAPower(Entity? above, Entity? below, Entity.Variable x) + { + if (above is null && below is null) + return x; + if (above is not null && below is not null) + return null; + var power = above ?? below!; + Entity @base, m; + if (power == x) + (@base, m) = (x, Entity.Number.Integer.One); + else if (power is Entity.Powf(var written, var exponent) && !exponent.ContainsNode(x)) + (@base, m) = (written, exponent); + else + return null; + if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _)) + return null; + var mPlusOne = (above is null ? -m + 1 : m + 1).InnerSimplified; + if (TreeAnalyzer.IsZero(mPlusOne)) + return null; + return MathS.Pow(@base, mPlusOne) / (slope * mPlusOne); + } + /// /// The gamma function's order , as the number it is where its symbols cancel: (-1 + n + 1)/n /// is 1, which inner simplification leaves as written, and (p x)^(n - 1) e^(x^n), whose integral is @@ -802,6 +916,10 @@ _ when IncompleteGammaBesideAPower(expr, x) is { } gamma => gamma, // neither elementary nor an error function. https://github.com/asc-community/AngouriMath/issues/2075 _ when APowerBesideAnExponentialOfAPower(expr, x) is { } gammaForm => gammaForm, + // And beside a power of a sine and a cosine, as the exponentials Euler's formula writes them. + // https://github.com/asc-community/AngouriMath/issues/2075 + _ when APowerBesideATrigonometricFunction(expr, x) is { } trigonometricGammaForm => trigonometricGammaForm, + Entity.Absf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => // ∫ |ax + b| dx = sgn(ax + b) * (ax + b)^2 / (2a) MathS.Signum(arg) * MathS.Pow(arg, 2) / (2 * a), diff --git a/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs index cae26c093..a318cef0a 100644 --- a/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/APowerBesideAnExponentialOfAPowerIntegralTest.cs @@ -97,6 +97,19 @@ public void ADefiniteIntegralFromZero(string integral, double value) public void APositiveRateIsTakenOutOfThePower() => Assert.Equal("-gamma(4/3, x)", DifferentiatesBack("x^(1/3)*e^(-x)").Substitute("C", 0).InnerSimplified.Stringize()); + /// + /// Beside a power of a sine and a cosine, by Euler's formula: each e^(i q u) beside the power is a gamma function + /// of an imaginary argument. Rubi's 4.1.10 and 4.1.12. + /// + [Theory] + [InlineData("(c + d*x)^m*sin(a + b*x)")] + [InlineData("(c + d*x)^m*cos(a + b*x)^2")] + [InlineData("x^m*sin(b*x)")] + [InlineData("x^(1/3)*sin(x)*cos(x)")] + [InlineData("x*(a + b*sin(c + d*x^3))")] + [InlineData("sin(c + d*x^3)^2")] + public void APowerBesideASineOrACosine(string integrand) => DifferentiatesBack(integrand); + /// /// Not where the integral is elementary, an exponential integral or an error function: those keep the form /// they had. @@ -106,6 +119,8 @@ public void APositiveRateIsTakenOutOfThePower() [InlineData("sqrt(x)*e^(-x)")] [InlineData("e^(b*x^2)")] [InlineData("e^(b*x)/x")] + [InlineData("x^2*sin(x)")] + [InlineData("sin(x)/x")] public void NotWhereTheIntegralIsElementaryOrAnErrorFunction(string integrand) => Assert.DoesNotContain("gamma(", DifferentiatesBack(integrand).Stringize()); }