Skip to content

Compute limits at infinity by Gruntz's algorithm (#231, #353) #694

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 2 commits into from
Aug 4, 2026
Merged

Compute limits at infinity by Gruntz's algorithm (#231, #353) #694

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
90a42d2
Compute limits at infinity by Gruntz's algorithm
Rafael-SOWNet Aug 4, 2026
6fc4949
Merge origin/master into feat/gruntz
Rafael-SOWNet Aug 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
353 changes: 353 additions & 0 deletions Sources/AngouriMath/Functions/Continuous/Limits/Gruntz/AsymptoticSeries.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,353 @@
//
// Copyright (c) 2019-2022 Angouri.
// AngouriMath is licensed under MIT.
// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
// Website: https://am.angouri.org.
//

using AngouriMath.Core.Multithreading;
using PeterO.Numbers;
using System;
using System.Collections.Generic;
using System.Linq;
using static AngouriMath.Entity;
using static AngouriMath.Entity.Number;

namespace AngouriMath.Functions.Algebra
{
/// <summary>
/// A power series in <c>w</c> around <c>w -> 0+</c>, truncated at a known exponent, whose
/// coefficients are arbitrary expressions not containing <c>w</c>.
/// </summary>
/// <remarks>
/// <para>
/// This is what Gruntz's algorithm expands in, and it is not a Taylor series: the
/// exponents may be negative and fractional, and the coefficients are symbolic rather
/// than numeric. What is wanted from it in the end is only the leading term, but the
/// leading terms of a sum can cancel -- which is the whole point of the algorithm, since
/// e^(x + e^-x) - e^x cancels to every order when its two parts are expanded separately
/// -- so the series has to be carried far enough to see past the cancellation, and how
/// far that is is not known in advance. The caller raises the order until a leading term
/// survives.
/// </para>
/// <para>
/// The value of <c>log(w)</c> is supplied from outside rather than left as a logarithm of
/// a symbol. Gruntz chooses w as an exponential, so its logarithm is an ordinary
/// expression in x, and expanding a logarithm needs it: without it the constant term of
/// <c>log(c * w^e)</c> cannot be told from zero.
/// </para>
/// </remarks>
internal sealed class AsymptoticSeries
{
/// <summary>Coefficients by exponent, ascending. None of them contain w.</summary>
internal SortedDictionary<ERational, Entity> Terms { get; }

/// <summary>Nothing is known at or beyond this exponent.</summary>
internal ERational Order { get; }

private AsymptoticSeries(SortedDictionary<ERational, Entity> terms, ERational order)
=> (Terms, Order) = (terms, order);

/// <summary>
/// The order of a series that is not truncated at all. A constant, a power of w, and
/// anything built from those by the four operations is known exactly, and saying so
/// matters: a series carrying a finite order is one whose terms beyond it have been
/// thrown away, and multiplying by it would throw away the other factor's terms too.
/// </summary>
[ConstantField] internal static readonly ERational Exact = ERational.FromInt32(1 << 24);

internal static AsymptoticSeries Constant(Entity value)
{
var terms = new SortedDictionary<ERational, Entity>();
if (!IsKnownZero(value))
terms[ERational.Zero] = value;
return new AsymptoticSeries(terms, Exact);
}

/// <summary>The one term <c>coefficient * w^power</c>.</summary>
internal static AsymptoticSeries Monomial(Entity coefficient, ERational power)
{
var terms = new SortedDictionary<ERational, Entity>();
if (!IsKnownZero(coefficient))
terms[power] = coefficient;
return new AsymptoticSeries(terms, Exact);
}

/// <summary>
/// The leading coefficient and its exponent, or <see langword="null"/> if every term
/// known so far cancelled and the series is indistinguishable from zero at this order.
/// </summary>
internal (Entity Coefficient, ERational Power)? LeadingTerm()
{
foreach (var term in Terms)
if (!IsKnownZero(term.Value))
return (term.Value, term.Key);
return null;
}

/// <summary>
/// Whether the coefficient is zero on the evidence available. A coefficient that
/// cannot be decided is kept, so that a leading term is never claimed for something
/// that might have cancelled.
/// </summary>
private static bool IsKnownZero(Entity coefficient)
{
if (coefficient == Integer.Zero)
return true;
var simplified = coefficient.InnerSimplified;
if (simplified == Integer.Zero)
return true;
return simplified.Evaled is Complex { IsZero: true };
}

private static SortedDictionary<ERational, Entity> Fresh() => new();

internal AsymptoticSeries WithOrder(ERational order)
{
if (order.CompareTo(Order) >= 0)
return this;
var terms = Fresh();
foreach (var term in Terms)
if (term.Key.CompareTo(order) < 0)
terms[term.Key] = term.Value;
return new AsymptoticSeries(terms, order);
}

internal AsymptoticSeries Add(AsymptoticSeries other)
{
var order = Order.CompareTo(other.Order) <= 0 ? Order : other.Order;
var terms = Fresh();
foreach (var term in Terms.Concat(other.Terms))
{
if (term.Key.CompareTo(order) >= 0)
continue;
terms[term.Key] = terms.TryGetValue(term.Key, out var already)
? (already + term.Value).InnerSimplified
: term.Value;
}
return new AsymptoticSeries(terms, order);
}

internal AsymptoticSeries Negate()
{
var terms = Fresh();
foreach (var term in Terms)
terms[term.Key] = (-term.Value).InnerSimplified;
return new AsymptoticSeries(terms, Order);
}

internal AsymptoticSeries Multiply(AsymptoticSeries other)
{
// Everything below the leading exponent of one factor is unknown in the product,
// so the order of the product is each order raised by the other's leading power.
var mine = LeadingTerm()?.Power ?? Order;
var theirs = other.LeadingTerm()?.Power ?? other.Order;
var order = Min(Order.Add(theirs), other.Order.Add(mine));
var terms = Fresh();
foreach (var left in Terms)
foreach (var right in other.Terms)
{
var power = left.Key.Add(right.Key);
if (power.CompareTo(order) >= 0)
continue;
var product = (left.Value * right.Value).InnerSimplified;
terms[power] = terms.TryGetValue(power, out var already)
? (already + product).InnerSimplified
: product;
}
return new AsymptoticSeries(terms, order);
}

private static ERational Min(ERational a, ERational b) => a.CompareTo(b) <= 0 ? a : b;

/// <summary>
/// The series divided by its own leading term, which leaves a series starting at 1.
/// </summary>
private AsymptoticSeries? Normalised(out Entity coefficient, out ERational power)
{
coefficient = Integer.Zero;
power = ERational.Zero;
if (LeadingTerm() is not var (c, e))
return null;
(coefficient, power) = (c, e);
var terms = Fresh();
foreach (var term in Terms)
if (term.Key.CompareTo(e) >= 0)
terms[term.Key.Subtract(e)] = (term.Value / c).InnerSimplified;
return new AsymptoticSeries(terms, Order.Subtract(e));
}

/// <summary>
/// The reciprocal, by the geometric series on what is left once the leading term is
/// taken out. Declines if the leading term cannot be found, since dividing by a
/// series that might be zero is not something to guess at.
/// </summary>
internal AsymptoticSeries? Invert(ERational requested)
{
if (Normalised(out var coefficient, out var power) is not { } unit)
return null;
// 1/(1 + d) = 1 - d + d^2 - ..., which terminates because d starts above w^0.
var rest = unit.Add(Constant(-1));
var sum = Constant(1);
var term = Constant(1);
for (var i = 0; i < MaxExpansionTerms; i++)
{
term = term.Multiply(rest.Negate());
if (term.LeadingTerm() is not var (_, at) || at.CompareTo(requested) >= 0)
break;
sum = sum.Add(term);
}
return sum.WithOrder(requested).Multiply(Monomial((1 / coefficient).InnerSimplified, power.Negate()));
}

/// <summary>How many terms of a geometric, exponential or logarithmic series to take.</summary>
private const int MaxExpansionTerms = 24;

/// <summary>
/// The exponential of the series. The argument must not run off to infinity -- a
/// leading exponent below zero would mean the exponential belongs to a faster
/// comparability class than w, which the algorithm has already ruled out by putting
/// it in the set of most rapidly varying subexpressions instead.
/// </summary>
internal AsymptoticSeries? Exponentiate(ERational requested)
{
var constant = Terms.TryGetValue(ERational.Zero, out var atZero) ? atZero : Integer.Zero;
var rest = Add(Constant(-constant));
if (rest.LeadingTerm() is var (_, least) && least.CompareTo(ERational.Zero) <= 0)
return null;
// exp(constant) is left unexpanded on purpose: it is a coefficient, and opening
// it up would reveal a function of a faster class than w.
var scale = MathS.Pow(MathS.e, constant).InnerSimplified;
var sum = Constant(1);
var term = Constant(1);
for (var i = 1; i < MaxExpansionTerms; i++)
{
term = term.Multiply(rest).Multiply(Constant(Rational.Create(1, i)));
if (term.LeadingTerm() is not var (_, at) || at.CompareTo(requested) >= 0)
break;
sum = sum.Add(term);
}
return sum.WithOrder(requested).Multiply(Constant(scale));
}

/// <summary>
/// The logarithm of the series, given the value of <c>log(w)</c>. The leading term
/// has to be found, since <c>log</c> of something that might be zero is not something
/// to guess at.
/// </summary>
internal AsymptoticSeries? Logarithm(Entity logarithmOfW, ERational requested)
{
if (Normalised(out var coefficient, out var power) is not { } unit)
return null;
var head = (MathS.Ln(coefficient) + Rational.Create(power) * logarithmOfW).InnerSimplified;
// log(1 + d) = d - d^2/2 + d^3/3 - ...
var rest = unit.Add(Constant(-1));
var sum = Constant(head);
var term = Constant(1);
for (var i = 1; i < MaxExpansionTerms; i++)
{
term = term.Multiply(rest);
if (term.LeadingTerm() is not var (_, at) || at.CompareTo(requested) >= 0)
break;
sum = sum.Add(term.Multiply(Constant(Rational.Create(i % 2 == 1 ? 1 : -1, i))));
}
return sum.WithOrder(requested);
}

/// <summary>
/// The series of <paramref name="expr"/> in <paramref name="w"/> about zero, carried
/// to <paramref name="order"/>, or <see langword="null"/> if some part of it is not
/// something this can expand.
/// </summary>
internal static AsymptoticSeries? Expand(Entity expr, Variable w, Entity logarithmOfW, ERational order)
{
MultithreadingFunctional.ExitIfCancelled();
if (!expr.ContainsNode(w))
return Constant(expr);
if (expr == w)
return Monomial(1, ERational.One);
switch (expr)
{
case Sumf(var left, var right):
return Both(left, right, (a, b) => a.Add(b));
case Minusf(var left, var right):
return Both(left, right, (a, b) => a.Add(b.Negate()));
case Mulf(var left, var right):
return Both(left, right, (a, b) => a.Multiply(b));
case Divf(var left, var right):
return Expand(right, w, logarithmOfW, order) is { } divisor
&& divisor.Invert(order) is { } inverted
&& Expand(left, w, logarithmOfW, inverted.Order) is { } dividend
? dividend.Multiply(inverted)
: null;
case Powf(var @base, var power) when @base == MathS.e:
return Expand(power, w, logarithmOfW, order)?.Exponentiate(order);
case Powf(var @base, var power) when !power.ContainsNode(w):
return RaiseToConstant(@base, power, w, logarithmOfW, order);
case Powf(var @base, var power):
// b^p with w in the exponent is exp(p log b), which the two cases above
// between them can expand.
return Expand(MathS.Pow(MathS.e, power * MathS.Ln(@base)), w, logarithmOfW, order);
case Logf(var @base, var antilog) when !@base.ContainsNode(w):
return Expand(antilog, w, logarithmOfW, order)?.Logarithm(logarithmOfW, order) is { } logarithm
? logarithm.Multiply(Constant((1 / MathS.Ln(@base)).InnerSimplified))
: null;
default:
return null;
}

AsymptoticSeries? Both(Entity left, Entity right, Func<AsymptoticSeries, AsymptoticSeries, AsymptoticSeries> combine)
=> Expand(left, w, logarithmOfW, order) is { } a && Expand(right, w, logarithmOfW, order) is { } b
? combine(a, b)
: null;
}

private static AsymptoticSeries? RaiseToConstant(
Entity @base, Entity power, Variable w, Entity logarithmOfW, ERational order)
{
var requested = order;
if (Expand(@base, w, logarithmOfW, order) is not { } expanded)
return null;
if (power.Evaled is not Rational exponent)
return null;
var wanted = exponent.ERational;
// A whole power is a product, which keeps the coefficients as they are.
if (wanted.IsInteger() && wanted.Abs().CompareTo(ERational.FromInt32(MaxExpansionTerms)) <= 0)
{
var times = wanted.ToEIntegerIfExact().ToInt32Checked();
if (times == 0)
return Constant(1);
if (times < 0)
return expanded.Invert(order) is { } inverted ? Repeated(inverted, -times) : null;
return Repeated(expanded, times);
}
// Otherwise take the leading term out and use the binomial series on the rest.
if (expanded.Normalised(out var coefficient, out var leading) is not { } unit)
return null;
var rest = unit.Add(Constant(-1));
var sum = Constant(1);
var term = Constant(1);
for (var i = 1; i < MaxExpansionTerms; i++)
{
// The binomial coefficient built up one factor at a time, so that it works
// for a fractional exponent as much as a whole one.
var factor = (Rational.Create(wanted) - (i - 1)) / i;
term = term.Multiply(rest).Multiply(Constant(factor.InnerSimplified));
if (term.LeadingTerm() is not var (_, at) || at.CompareTo(requested) >= 0)
break;
sum = sum.Add(term);
}
sum = sum.WithOrder(requested);
var scaled = MathS.Pow(coefficient, power).InnerSimplified;
return sum.Multiply(Monomial(scaled, wanted.Multiply(leading)));
}

private static AsymptoticSeries Repeated(AsymptoticSeries series, int times)
{
var result = Constant(1);
for (var i = 0; i < times; i++)
result = result.Multiply(series);
return result;
}
}
}
Loading
Loading