Skip to content

Commit c32fdaa

Browse files
Decline four slow routes ahead of the answering one by value, for the hyperbolic integral functions and the tangent (#2021)
Four routes spent gigabytes declining before the one that answers, on Rubi's 8.5 and 4.3.3.1: - SolveBySplittingSum expanded Shi(u) or Chi(u) times (e^u ± e^-u)/x^n into exponential terms, each non-elementary on its own, and declined after 10-13 GB; SolveIntegratingByParts then answered the whole in 0.2-0.9 GB. IsAHyperbolicIntegralBesideItsExponentialsOverAPolynomial sends such an integrand to it first, the same call at the same depth. - SolveByScalingTheVariable ran three Simplify calls before finding that x = c t leaves c inside e^(b x) or Shi(u). TheScaleStaysInsideAFunction evaluates the scaled integrand at two scales and two points first; it can only decline. - SolveBySubstitution simplified each quotient by a special function's derivative, a gigabyte each for x^2 Chi(a + b x)^2. A special function is transcendental over the elementary ones, so with it held a quotient taking two values at two points keeps x through any rewriting: VariesWhereASpecialFunctionIsHeld refuses it before the simplification. - SolveByTangentSubstitution's cot^k leaves 1^k above the bar and u^k below, outside the lowest-terms class that skips Simplify; it is in it now. Answers are identical. Chain against master f7f2a7a: pocket of 3,919 rows (4.5.x, 4.3.3.x and 8.x every 2nd row, the rest of 4.3 and 4.4 every 6th, every 30th elsewhere), verdicts identical, allocation over the rows solved on both 460.9 -> 435.0 GB, 8.5 alone 31.0 -> 11.8 GB; samples unchanged; no row's allocation grew past 1.3x. By the probe against d3a6bce: 8.5:62 12.59 -> 0.17 GB, 8.5:154 12.68 -> 0.19 GB, 4.3.3.1:313 5.01 -> 0.20 GB. 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 6e16aae commit c32fdaa

2 files changed

Lines changed: 266 additions & 2 deletions

File tree

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

Lines changed: 182 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,14 @@ internal static class IndefiniteIntegralSolver
7272
// parts on what it hands down was a two-minute search.
7373
if (integrateByParts && IsAnInverseFunctionTimesARationalFunction(expr, x) && SolveIntegratingByParts(expr, x) is { } byParts)
7474
return byParts;
75+
// And a hyperbolic sine or cosine integral of a linear u beside a sum of e^u and e^(-u),
76+
// which is how a hyperbolic function of u beside it is written, over a polynomial: parts
77+
// against the whole pair the integral's derivative, sinh(u)/u, with the whole hyperbolic
78+
// function, where expanded each exponential leaves `Ei(u) e^(-u)/u` or the like, which
79+
// has no closed form, and `cosh(b x) Shi(b x)/x^2` spent thirteen gigabytes declining its
80+
// two terms before parts answered it whole in a quarter of one.
81+
if (integrateByParts && IsAHyperbolicIntegralBesideItsExponentialsOverAPolynomial(expr, x) && SolveIntegratingByParts(expr, x) is { } byPartsWhole)
82+
return byPartsWhole;
7583
var splitted = TreeAnalyzer.GatherLinearChildrenOverSumAndExpand(expr, e => e.ContainsNode(x));
7684
if (splitted is null || splitted.Count < 2) return null; // nothing to do, let other solvers do the work
7785
// Each expanded term with a factor written on both sides of its bar cancelled:
@@ -428,6 +436,44 @@ private static Entity WithTheFactorsEveryTermSharesOut(Entity sum, Entity.Variab
428436
/// <summary>The largest term a sum is split over as written, before the expanding gather.</summary>
429437
private const int LargestTermTakenAsWritten = 60;
430438

439+
/// <summary>
440+
/// Whether <paramref name="expr"/> is <c>Shi(u)</c> or <c>Chi(u)</c> for a linear <c>u</c>
441+
/// times sums of constant multiples of <c>e^u</c> and <c>e^(-u)</c>, over a polynomial in
442+
/// <paramref name="x"/> that is not a constant, with nothing else in <paramref name="x"/>.
443+
/// </summary>
444+
private static bool IsAHyperbolicIntegralBesideItsExponentialsOverAPolynomial(Entity expr, Entity.Variable x)
445+
{
446+
if (expr is Sumf or Minusf)
447+
return false;
448+
var (above, below) = Functions.SingleQuotient.Of(expr);
449+
if (!below.ContainsNode(x) || !TreeAnalyzer.TryGetPolynomial(below, x, out _))
450+
return false;
451+
var factors = Mulf.LinearChildren(above).Where(factor => factor.ContainsNode(x)).ToList();
452+
if (factors.Count(factor => factor is Shif or Chif) != 1)
453+
return false;
454+
var argument = factors.First(factor => factor is Shif or Chif).DirectChildren.First();
455+
if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out var offset))
456+
return false;
457+
bool AnExponentialOfTheArgument(Entity node)
458+
=> node is Powf(var @base, var exponent) && @base == MathS.e
459+
&& TreeAnalyzer.TryGetPolyLinear(exponent, x, out var rate, out var shift)
460+
&& (VanishesIdentically(rate - slope) && VanishesIdentically(shift - offset)
461+
|| VanishesIdentically(rate + slope) && VanishesIdentically(shift + offset));
462+
bool ATermInTheExponentials(Entity term)
463+
=> Mulf.LinearChildren(term).Count(factor => factor.ContainsNode(x)) == 1
464+
&& Mulf.LinearChildren(term).Where(factor => factor.ContainsNode(x)).All(AnExponentialOfTheArgument);
465+
var sums = 0;
466+
foreach (var factor in factors)
467+
{
468+
if (factor is Shif or Chif)
469+
continue;
470+
if (factor is not (Sumf or Minusf) || !Sumf.LinearChildren(factor).All(term => !term.ContainsNode(x) || ATermInTheExponentials(term)))
471+
return false;
472+
sums++;
473+
}
474+
return sums > 0;
475+
}
476+
431477
/// <summary>
432478
/// Whether the integrand is one logarithm or inverse trigonometric function of
433479
/// <paramref name="x"/>, or a whole power of one, times a rational function of
@@ -13830,7 +13876,11 @@ static bool IsTwice(Entity argument, Entity.Variable x)
1383013876
: Functions.SingleQuotient.Combine(inU / (1 + MathS.Sqr(uSub))) is var overOneBar
1383113877
&& IsAQuotientOfDistinctSumsOfSymbolsInLowestTerms(overOneBar, uSub)
1383213878
? WithNeighboursToOnePowerGathered(overOneBar)
13833-
: overOneBar.Simplify();
13879+
// And so beside a power of u below the bar, what a power of the cotangent leaves
13880+
// with its ones above: see AQuotientOfDistinctSumsOverAPowerOfTheVariable.
13881+
: AQuotientOfDistinctSumsOverAPowerOfTheVariable(overOneBar, uSub) is { } overThePower
13882+
? overThePower
13883+
: overOneBar.Simplify();
1383413884
if (integrand is Providedf(var withoutCondition, _))
1383513885
integrand = withoutCondition;
1383613886
// Asked as the same question: the rewrite renames the variable and adds no step of
@@ -13922,6 +13972,57 @@ private static bool IsAQuotientOfDistinctSumsOfSymbolsInLowestTerms(Entity expr,
1392213972
return holdsASymbol && ReferenceEquals(WithNeighboursToOnePowerGathered(top), top);
1392313973
}
1392413974

13975+
/// <summary>
13976+
/// <paramref name="expr"/>, a quotient of sums as <see cref="IsAQuotientOfDistinctSumsOfSymbolsInLowestTerms"/>
13977+
/// reads them beside one positive whole power of <paramref name="u"/> below the bar, with the
13978+
/// ones above it -- <c>1</c> or a power of it -- left out, as the simplifier hands it back;
13979+
/// null where it is not such a quotient.
13980+
/// </summary>
13981+
/// <remarks>
13982+
/// A power of the cotangent is <c>(1/u)^k</c> under the tangent, and over one bar it is
13983+
/// <c>1^k</c> above and <c>u^k</c> below: Rubi's 4.3.3.1
13984+
/// <c>cot(c + d x)^7 (a + b tan(c + d x))^4 (A + B tan(c + d x))</c> is
13985+
/// <c>1^7 (a + b u)^4 (A + B u)/(u^7 (1 + u^2))</c>, which the simplifier spent 4.8 of the
13986+
/// integral's 5 GB on to hand back without the <c>1^7</c>. The sums have constant terms, so
13987+
/// the power of u divides none of them and the quotient is in lowest terms as before. Only
13988+
/// with an exponent no other power below the bar has, since the gathering of neighbours to
13989+
/// one power was measured for sums alone.
13990+
/// https://github.com/asc-community/AngouriMath/issues/233
13991+
/// </remarks>
13992+
private static Entity? AQuotientOfDistinctSumsOverAPowerOfTheVariable(Entity expr, Entity.Variable u)
13993+
{
13994+
if (expr is not Divf(var top, var bottom))
13995+
return null;
13996+
static bool IsOne(Entity factor) => factor == Number.Integer.One || factor is Powf(var one, Number.Integer) && one == Number.Integer.One;
13997+
var above = Mulf.LinearChildren(top).Where(factor => !IsOne(factor)).ToList();
13998+
var below = Mulf.LinearChildren(bottom);
13999+
if (above.Count == 0)
14000+
return null;
14001+
Entity? power = null;
14002+
var exponent = EInteger.One;
14003+
var others = new List<Entity>();
14004+
foreach (var factor in below)
14005+
{
14006+
if (factor == u || factor is Powf(var @base, Number.Integer { EInteger.Sign: > 0 }) && @base == u)
14007+
{
14008+
if (power is not null)
14009+
return null;
14010+
power = factor;
14011+
exponent = factor is Powf(_, Number.Integer { EInteger: var whole }) ? whole : EInteger.One;
14012+
continue;
14013+
}
14014+
others.Add(factor);
14015+
}
14016+
if (power is null || others.Count == 0
14017+
|| others.Any(factor => factor is Powf(_, Number.Integer { EInteger: var whole }) && whole.Equals(exponent)))
14018+
return null;
14019+
static Entity Product(IEnumerable<Entity> factors) => factors.Aggregate((left, right) => left * right);
14020+
var withoutThePower = Product(above) / Product(others);
14021+
return IsAQuotientOfDistinctSumsOfSymbolsInLowestTerms(withoutThePower, u)
14022+
? WithNeighboursToOnePowerGathered(Product(above) / bottom)
14023+
: null;
14024+
}
14025+
1392514026
/// <summary>
1392614027
/// <paramref name="expr"/> with each product of two powers to one whole exponent written
1392714028
/// as that power of their product, from the leaves up: <c>(P^2 Q^2) R^2</c> is
@@ -30724,6 +30825,45 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons
3072430825
return scaled;
3072530826
}
3072630827

30828+
/// <summary>
30829+
/// Whether <paramref name="expr"/> under <c>x = c t</c>, for <c>c</c> the
30830+
/// <paramref name="scale"/>, is shown by its values not to be a function of the scale times
30831+
/// one of <c>t</c>: where a function that is not algebraic has an argument holding the
30832+
/// variable, the quotient of the scaled integrand's values at two points of <c>t</c> is taken
30833+
/// at two values of the scale, and the two quotients differ. False wherever a value cannot be
30834+
/// had, the values agree, or no such function is written.
30835+
/// </summary>
30836+
private static bool TheScaleStaysInsideAFunction(Entity expr, Entity.Variable x, Entity.Variable scale)
30837+
{
30838+
if (!expr.Nodes.Any(node => node switch
30839+
{
30840+
Sumf or Minusf or Mulf or Divf or Number or Variable => false,
30841+
Powf(_, var exponent) => exponent.ContainsNode(x),
30842+
_ => node.DirectChildren.Any(child => child.ContainsNode(x)),
30843+
})
30844+
|| PowerOfXAsWritten(expr, x) > LargestPowerEvaluatedForARepeat)
30845+
return false;
30846+
// In decimals: downcast, every value is folded exactly.
30847+
using var _ = MathS.Settings.DowncastingEnabled.Set(false);
30848+
Number.Complex? At(string scaleValue, string point)
30849+
{
30850+
var c = Number.Real.Create(EDecimal.FromString(scaleValue));
30851+
var value = (expr.Substitute(x, c * Number.Real.Create(EDecimal.FromString(point))).Substitute(scale, c) * c).Evaled;
30852+
return value is Number.Complex { IsNaN: false } number && number.IsFinite && !number.IsZero ? number : null;
30853+
}
30854+
if (At("1", "0.83") is not { } oneFirst || At("1", "1.29") is not { } oneSecond
30855+
|| At("1.37", "0.83") is not { } scaledFirst || At("1.37", "1.29") is not { } scaledSecond)
30856+
return false;
30857+
var first = ToComplex(scaledFirst) / ToComplex(oneFirst);
30858+
var second = ToComplex(scaledSecond) / ToComplex(oneSecond);
30859+
var magnitude = System.Math.Max(1.0, System.Numerics.Complex.Abs(first));
30860+
var difference = System.Numerics.Complex.Abs(first - second);
30861+
return !double.IsNaN(difference) && !double.IsInfinity(difference) && !double.IsNaN(magnitude) && !double.IsInfinity(magnitude)
30862+
&& difference > 1e-6 * magnitude;
30863+
static System.Numerics.Complex ToComplex(Number.Complex number)
30864+
=> new(number.RealPart.EDecimal.ToDouble(), number.ImaginaryPart.EDecimal.ToDouble());
30865+
}
30866+
3072730867
/// <summary>
3072830868
/// An integrand carrying one symbolic parameter, scaled by it — <c>x = c t</c> — so that
3072930869
/// what is left to integrate has the variable alone in it.
@@ -30775,6 +30915,14 @@ bool IsAConstantTimesAPowerOfOneMinusUSquared(Entity polynomial, out Entity cons
3077530915
}
3077630916
if (scale is null)
3077730917
return null;
30918+
// A function that is not algebraic, of an argument holding the variable, takes the scale
30919+
// inside it under x = c t unless the argument divides it out -- `e^(b x)` is `e^(b^2 t)`
30920+
// and `Shi(x)` is `Shi(b t)` -- and then no power of it comes out, which two values of the
30921+
// scale say before the three simplifications below: they spent a gigabyte and more on
30922+
// `(b x e^(b x) + e^(b x)) (b x Shi(b x) - cosh(b x))`, what parts leaves of `x Shi(b x)^2`, to
30923+
// find it declined. Only a decline is decided so; agreeing values leave it to the check below.
30924+
if (TheScaleStaysInsideAFunction(expr, x, scale))
30925+
return null;
3077830926

3077930927
// dx = c dt.
3078030928
var t = Variable.CreateUnique(expr, "u_scale");
@@ -38427,7 +38575,7 @@ private static Dictionary<Variable, double> PinsFor(IEnumerable<Variable> symbol
3842738575
/// or a numeric power, where a symbol is not pinned, and where a value is not real and
3842838576
/// finite. A root of a negative number is not real here, odd or even.
3842938577
/// </summary>
38430-
private static (double Value, double Scale)? AtRealPoint(Entity expr, Dictionary<Variable, double> pins, Entity.Variable x, double at)
38578+
private static (double Value, double Scale)? AtRealPoint(Entity expr, Dictionary<Variable, double> pins, Entity.Variable x, double at, bool exponentials = false)
3843138579
{
3843238580
var scale = 0.0;
3843338581
double Of(Entity node)
@@ -38448,6 +38596,7 @@ double Of(Entity node)
3844838596
=> System.Math.Pow(Of(@base), whole.EInteger.ToInt32Checked()),
3844938597
Powf(var @base, Number.Rational fraction) => Of(@base) is var inner && inner >= 0
3845038598
? System.Math.Pow(inner, fraction.ERational.ToDouble()) : double.NaN,
38599+
Powf(var @base, var exponent) when exponentials && @base == MathS.e => System.Math.Exp(Of(exponent)),
3845138600
_ => double.NaN,
3845238601
};
3845338602
if (System.Math.Abs(value) > scale)
@@ -38516,6 +38665,29 @@ bool Constant(System.Func<int, double> of)
3851638665
return false;
3851738666
}
3851838667

38668+
/// <summary>
38669+
/// <see cref="VariesWithTheVariableWhereTheCandidateIsHeld"/> for a quotient by the derivative
38670+
/// of a special function, with that function written as <paramref name="uSub"/>: in doubles,
38671+
/// where it is a rational function of x, the symbols and exponentials, of any size.
38672+
/// </summary>
38673+
/// <remarks>
38674+
/// A special function is transcendental over the elementary functions -- it is not one of
38675+
/// them, and what is algebraic over them is -- so a quotient elementary in x and rational in
38676+
/// the function is a function of the function alone only where, with it held, no x is left:
38677+
/// two values at two points then say that no rewriting that keeps the value takes the x out.
38678+
/// Real values at two positive points only, as there.
38679+
/// </remarks>
38680+
private static bool VariesWhereASpecialFunctionIsHeld(Entity expr, Variable uSub, Entity.Variable x)
38681+
{
38682+
if (HoldsTheImaginaryUnit(expr) || !expr.ContainsNode(uSub)
38683+
|| !expr.Nodes.All(node => node is Sumf or Minusf or Mulf or Divf or Number or Variable || node is Powf(_, Number) || node is Powf(var @base, _) && @base == MathS.e))
38684+
return false;
38685+
var pins = PinsFor(expr.Vars, x);
38686+
return AtRealPoint(expr, pins, x, 0.29, exponentials: true) is (var first, var scaleFirst)
38687+
&& AtRealPoint(expr, pins, x, 0.53, exponentials: true) is (var second, var scaleSecond)
38688+
&& System.Math.Abs(first - second) > 1e-6 * System.Math.Max(1.0, System.Math.Max(scaleFirst, scaleSecond));
38689+
}
38690+
3851938691
/// <summary>
3852038692
/// Whether <paramref name="expr"/>, written in <paramref name="x"/> and a candidate's variable
3852138693
/// <paramref name="uSub"/>, takes two real values at two points of x, <paramref name="uSub"/>
@@ -42280,6 +42452,14 @@ private static bool IsAFunctionOfTheTangentAlone(Entity expr, Entity.Variable u)
4228042452
// simplification each time it was offered.
4228142453
if (KeepsTheVariableThroughEveryRewriting(inTermsOfU, u, uSub, x, source, duDx))
4228242454
continue;
42455+
// And under a special function, where with it held the quotient takes two values at two
42456+
// points of x, refused before it is simplified, for the reason in
42457+
// VariesWhereASpecialFunctionIsHeld: `x^2 Chi(a + b x)^2` offered `Chi(a + b x)` to three
42458+
// remainders of its parts, and simplifying each quotient to refuse it was a gigabyte.
42459+
if (IsASpecialFunction(u) && inTermsOfU.ContainsNode(x)
42460+
&& !inTermsOfU.Nodes.Any(node => IsASpecialFunction(node) && node.ContainsNode(x))
42461+
&& VariesWhereASpecialFunctionIsHeld(inTermsOfU, uSub, x))
42462+
continue;
4228342463
integrandInU = !inTermsOfU.ContainsNode(x)
4228442464
&& TreeAnalyzer.TryGetPolyLinear(u, x, out var slope, out var offset) && !slope.ContainsNode(x) && !offset.ContainsNode(x)
4228542465
&& IsAFunctionOfTheTangentAlone(inTermsOfU, uSub) && !HoldsTheImaginaryUnit(inTermsOfU)

0 commit comments

Comments
 (0)