Skip to content

Commit 2bec4ec

Browse files
Flatten a whole power of a constant times an exponential (#1942)
`(-e^(-b x))^2` is `e^(-2 b x)`, and nothing read it as one: the expansion of `sinh(a + b x)^3 sinh(c + d x)^3` writes eight such products and the integral stopped there, where 2.5.0 answered it by parts. `SolveByFlatteningAPowerOfAnExponential` now distributes a whole power over a constant times an exponential, `(k E)^n = k^n E^n`, which holds for any constant k when n is whole. Rubi 6.7.1 row 195 answers again, and rows 327, 347, 365 and 366, which timed out, answer in 80-210 ms. Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura Co-authored-by: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
1 parent da06bac commit 2bec4ec

2 files changed

Lines changed: 80 additions & 8 deletions

File tree

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

Lines changed: 26 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -22156,17 +22156,35 @@ private static (Entity Constant, Entity Base, Entity ExponentInU, Variable U, En
2215622156
/// exponential of the product: <c>(e^x)^(1/3)</c> is <c>e^(x/3)</c> and <c>(3^(3x))^(1/4)</c>
2215722157
/// is <c>3^(3x/4)</c>, exactly, since the base is positive; and it is only in that
2215822158
/// spelling that the exponential rules read them. Timofeev's
22159-
/// <c>(cos(x/2) + sin(x/2))/(e^x)^(1/3)</c> and <c>cos(3x/2)/(3^(3x))^(1/4)</c>.
22159+
/// <c>(cos(x/2) + sin(x/2))/(e^x)^(1/3)</c> and <c>cos(3x/2)/(3^(3x))^(1/4)</c>. And a whole
22160+
/// power of a constant times one, <c>(-e^(-b x))^2</c>, as the constant's power times the
22161+
/// exponential of the multiple, whatever the constant is.
2216022162
/// </summary>
2216122163
internal static Entity? SolveByFlatteningAPowerOfAnExponential(Entity expr, Entity.Variable x, bool integrateByParts)
2216222164
{
22163-
var flattened = expr.Replace(node =>
22164-
node is Powf(Powf(var @base, var inner), var outer)
22165-
&& !@base.ContainsNode(x) && inner.ContainsNode(x) && !outer.ContainsNode(x)
22166-
&& (@base == MathS.e || @base.Evaled is Number.Real { IsPositive: true })
22167-
&& outer.Evaled is Number.Real
22168-
? MathS.Pow(@base, (outer * inner).InnerSimplified)
22169-
: node);
22165+
bool IsAnExponential(Entity @base, Entity inner) =>
22166+
!@base.ContainsNode(x) && inner.ContainsNode(x)
22167+
&& (@base == MathS.e || @base.Evaled is Number.Real { IsPositive: true });
22168+
// A whole power of a constant times an exponential: (k E)^n is k^n E^n for every k and
22169+
// E when n is whole, so (-e^(-b x))^2 is e^(-2 b x). The product of two of them, which
22170+
// sinh(b x)^3 sinh(d x)^3 expands into, is read by nothing else.
22171+
Entity? Distributed(Entity left, Entity right, Number.Integer whole) =>
22172+
(left, right) switch
22173+
{
22174+
(Powf(var b, var p), var k) when !k.ContainsNode(x) && IsAnExponential(b, p) =>
22175+
MathS.Pow(k, whole) * MathS.Pow(b, (whole * p).InnerSimplified),
22176+
(var k, Powf(var b, var p)) when !k.ContainsNode(x) && IsAnExponential(b, p) =>
22177+
MathS.Pow(k, whole) * MathS.Pow(b, (whole * p).InnerSimplified),
22178+
_ => null
22179+
};
22180+
var flattened = expr.Replace(node => node switch
22181+
{
22182+
Powf(Powf(var @base, var inner), var outer)
22183+
when IsAnExponential(@base, inner) && !outer.ContainsNode(x) && outer.Evaled is Number.Real =>
22184+
MathS.Pow(@base, (outer * inner).InnerSimplified),
22185+
Powf(Mulf(var left, var right), Number.Integer whole) => Distributed(left, right, whole) ?? node,
22186+
_ => node
22187+
});
2217022188
// The same question in another spelling, asked as one of its own, so that the
2217122189
// closed rules that answer only at the top -- the exponential times a
2217222190
// trigonometric -- are consulted for it.
Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
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 whole power of a constant times an exponential, <c>(k e^u)^n = k^n e^(n u)</c>. The terms of
16+
/// <c>sinh(a + b x)^3 sinh(c + d x)^3</c>, written in exponentials, are products of two of them,
17+
/// <c>(-e^(-b x))^2 (-e^(-d x))^2</c>, which nothing else read: 2.5.0 answered the product, and
18+
/// 2.6.0's development had stopped. Rubi's 6.7.1.
19+
/// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
20+
/// </summary>
21+
[Trait("Area", "Calculus")]
22+
public sealed class AWholePowerOfAMultipleOfAnExponentialIntegralTest
23+
{
24+
[Theory]
25+
[InlineData("(-e^(-b*x))^2*(-e^(-d*x))^2")]
26+
[InlineData("(-e^(-b*x))^2")]
27+
[InlineData("(2*e^x)^3")]
28+
[InlineData("sinh(b*x)^3*sinh(d*x)^3")]
29+
[InlineData("sinh(a + b*x)^3*sinh(c + d*x)^3")]
30+
public void AsAnExponentialOfTheMultiple(string integrand)
31+
{
32+
var integral = integrand.ToEntity().Integrate("x");
33+
var text = integral.Stringize();
34+
Assert.DoesNotContain("integral(", text);
35+
Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}");
36+
Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.3);
37+
var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
38+
var original = Pinned(integrand.ToEntity());
39+
var compared = 0;
40+
foreach (var at in new[] { -1.2, -0.7, 0.3, 0.8, 1.3, 2.9 })
41+
{
42+
var want = original.Substitute("x", at).EvalNumerical();
43+
var got = derivative.Substitute("x", at).EvalNumerical();
44+
if (want.IsNaN)
45+
continue;
46+
compared++;
47+
Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
48+
< 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)),
49+
$"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
50+
}
51+
Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}");
52+
}
53+
}
54+
}

0 commit comments

Comments
 (0)