Implement an algorithm for limits

[WhiteBlackGoose]

It's a great-looking approach, located here.

[Rafael-SOWNet]

An algorithm for limits landed in #694: Gruntz's method — the MRV (most-rapidly-varying) set approach from his thesis, with an asymptotic series engine underneath it.

I could not open the ResearchGate link to check it is the same paper you meant, so worth confirming. What is implemented is the one that computes a limit by finding the most rapidly varying subexpression, rewriting the whole expression in terms of it, and expanding as a series — not a pattern table.

What it buys over the rules that were already there is the cases nothing else has any reading of, in particular differences whose terms cancel to every order:

lim x->+oo (e^(x + e^(-x)) - e^x)                  // 1
lim x->+oo (sqrt(x^2 + 3x) - sqrt(x^2 + 1))        // 3/2
lim x->+oo (x^20 / e^x)                            // 0, beyond l'Hopital's step bound

It runs last of the readings at infinity, after the cheaper ones, since it is the most expensive.

Leaving this open in case the paper you linked is a different algorithm and you want that one too — tell me which and I will read it. If Gruntz is what you meant, it can be closed.

[Rafael-SOWNet]

This looks done. Measured on master at a45a7256.

The paper linked in the body is the Gruntz algorithm, and it is implemented — Sources/AngouriMath/Functions/Continuous/Limits/Gruntz/ holds Gruntz.cs and AsymptoticSeries.cs. It answers the cases that distinguish a real algorithm from a table of special cases:

limit at +oo	answer
e^x / x^10	+oo
ln(x) / x	0
(1 + 1/x)^x	e
e^x - x	+oo

The last is #231's one unticked case, which also now works.

I have not audited it against the paper section by section — that is a larger claim than I can make from behaviour — so treat this as "the algorithm named here exists and answers the hard cases" rather than "the implementation is complete and faithful". But the issue as filed asks for the algorithm to be implemented, and it is.

[Happypig375]

Decide whether to close after reading the paper and comparing with our implementation

[Rafael-SOWNet]

I read the paper and compared it with our implementation, as you asked.

Which paper. The link is Gonnet and Gruntz, Limit Computation in Computer Algebra, ETH Zürich technical report 187, November 1992. ResearchGate refuses the PDF (HTTP 403), and ETH's Research Collection search needs JavaScript, so I read Bruno Salvy's summary of Gruntz's 29 March 1993 talk at the INRIA Algorithms seminar (gruntz.pdf), which states the algorithm step by step. The implementation cites Gruntz's 1996 thesis, which is the same algorithm worked out in full.

Step by step, it is what Limits/Gruntz/Gruntz.cs does:

The paper	Ours
mrv(f + g), mrv(f g) are the larger of the parts'	Mrv → MrvMax, for sums, differences, products and quotients
mrv(f^c) = mrv(f), mrv(log f) = mrv(f), mrv(x) = x	the matching arms of Mrv
mrv(e^f) contains e^f exactly when f → ∞	decided by LimitInf of the exponent, for either sign of infinity
compare f and g by `lim log	f
when one of them is x, replace x by exp(x) throughout	the lifting in MrvLeadTerm
rewrite the set in terms of one member ω	Rewrite: the member containing none of the others, inverted if needed so that ω → 0
Puiseux expansion in ω, other terms as parameters; recurse on the leading coefficient when ω's exponent is 0	AsymptoticSeries.Expand (rational exponents); LimitInf's exponent-0 branch
the problem reduces to deciding whether an exp-log function is zero	the zero test in AsymptoticSeries; outside the exp-log class it declines rather than guesses

But it did not answer the paper's own example. The paper introduces the algorithm with e^x (exp(1/x + e^(-x)) - exp(1/x)) as x → ∞, where expanding in powers of 1/x gives O(x^-k) for every k. On master and on 2.5.0 it came back as written, while e^(x + e^(-x)) - e^x, x^x / e^(x ln x) and e^(e^x) / e^(e^x - e^(-e^x)) were all answered. GRUNTZ_DEBUG=1 shows the algorithm doing the paper's steps correctly: the set {e^(-x), e^x} and the rewrite (e^(w + 1/x) - e^(1/x))/w are right, and then the series engine declines.

The cause is that a constant term was taken off by subtracting it. Exponentiate subtracts the constant term 1/x, leaving 1/x + -(1/x) at w^0. Normalising a series divides by the leading coefficient and subtracts 1, leaving x/x + -1. The inner simplification writes out neither as zero, and the zero test cannot decide a coefficient that contains x. So the leading exponent read 0 where it is not, and Exponentiate declined. The same x/x + -1 is what makes the geometric series swell in the trace.

Fixed in #1627. The constant term is now removed exactly, the normalised leading coefficient is 1 exactly, and the zero test collects like terms and reads 0 provided not x = 0 as zero. The example is now 1, as are two variants with a different scale and parameter. For the whole computation, GRUNTZ_DEBUG printed 7981 lines before the fix and prints 97 now. The 643 limit tests are unchanged, the suite passes 14153/0, and the gate passes.

So the implementation is the paper's algorithm. With that fix it also answers the example the paper is built around, and I'll close this once the fix merges, unless you'd rather keep it open for something the full report covers and the summary does not.

[Rafael-SOWNet]

Closing, as I said above: the implementation is the paper's algorithm, and with #1627 (371b38e) it answers the example the paper opens with. If you want the full technical report read against it, reopen and I'll find a copy that isn't behind ResearchGate.