Skip to content

Commit 8d4af56

Browse files
Rafael-SOWNetRafael-SOWNet
authored
Apply l'Hopital's rule at infinity too (#596) (#680)
* Apply l'Hopital's rule at infinity too (#596) The rule was only ever reached when the destination was finite, so the ordinary competing-growth limits had nothing to catch them: lim x->+oo ln(x) / x was unevaluated, and NaN once simplified lim x->+oo x / e^x " lim x->+oo ln(x) / sqrt(x) " lim x->+oo e^x / (e^x + 1) " Three things had to be settled before the rule could be let loose there. It has to stop. Differentiating both parts does not always make the quotient simpler: x / sqrt(x^2 + 1) turns into sqrt(x^2 + 1) / x, which turns back, so neither form ever settles. The recursion is bounded, at a depth that still leaves room for the degree of a polynomial, since x^10 / e^x takes ten steps. It has to leave a wrong answer no worse. NaN asserts that a limit does not exist, which is a stronger claim than the rule can make when it has not settled, so the rule's answer is only taken when it is not NaN; otherwise whatever was there before stands, unevaluated included. It has to recognise the quotient. Two shapes were being turned away: - Differentiation does not produce quotients. d/dx ln(x)^2 is 2 * ln(x) * (1 / x), a product with a quotient inside rather than 2 * ln(x) / x, and only the second has a limit here, so the derivative quotient is simplified before it is used. The domain condition that simplification leaves behind is dropped with it: a Providedf is not a continuous node and would be turned away unread, and limits already treat the expression as continuous. - A product with a reciprocal factor is a quotient. That is the shape in #596, where lim x->+oo x^4 * e^(-x) came out as NaN while x^4 / e^x gave 0. Corpus 88/117 -> 91/117 with no wrong answers: lim:oo/oo goes to 4 out of 4 and x / ln(x)^2 joins it. 3909 unit tests and 127 F# tests pass, and the property checker's 1320 checks are clean. Still out of reach and still left unevaluated rather than guessed: e^x - x, sqrt(x^2 + x) - x and e^(x + e^(-x)) - e^x, none of which is a quotient. * Bound the work l'Hopital's rule does, not just its depth The rule as added in this branch could run for more than twenty seconds on lim x -> +oo sqrt(x^2 - x) / x and on x^(3/2) * sqrt(1 + 1/x^2) / x^2, both of which master answers, unevaluated, in under two hundred milliseconds. The bound on the depth was not a bound on the work. Every step asks what the two parts of its quotient tend to before differentiating them, and each of those is a limit the rule may be applied to in turn, so the steps fan out rather than follow one another and sixteen deep is far more than sixteen steps. Two conditions stop the fan-out at its source: the rule declines a quotient it is already partway through, since differentiating sqrt(x^2 - x) / x gives its reciprocal and differentiating that gives the first back; and it declines one that has grown by more than the derivative of a power or a logarithm accounts for, since a bigger quotient makes the two questions each step asks harder than the one being answered. The measured room is wide: a step on the way to an answer adds two nodes and the most any was seen to add is six, against eighteen for the quotient above. A cap on the total number of applications is left behind them as a backstop. Depth is now held for the whole of the rule rather than only for its recursive call, so the two limits each step asks for are counted as well. Both forms now terminate, the second of them in 86 ms, and the two are pinned. Every limit the branch already answered is unchanged and no faster or slower: x^10 / e^x still takes ten steps and 40 ms. 3911 tests, corpus 91/117 with nothing wrong, in error or timing out.
1 parent 7cbdc3e commit 8d4af56

3 files changed

Lines changed: 232 additions & 4 deletions

File tree

‎Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Solvers.Definition.cs‎

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,9 +62,21 @@ private static Entity ExpandLogarithm(Entity expr)
6262

6363
MultithreadingFunctional.ExitIfCancelled();
6464
if (!dest.IsFinite)
65+
{
6566
// just compute limit with no check for left/right equality
6667
// here approach left will be ignored anyways, as dest is infinite number
67-
return expr.ComputeLimitDivideEtImpera(x, dest, ApproachFrom.Left);
68+
var atInfinity = expr.ComputeLimitDivideEtImpera(x, dest, ApproachFrom.Left);
69+
if (atInfinity is { } found && found.Evaled != MathS.NaN)
70+
return found;
71+
// l'Hopital's rule was only ever reached for a finite destination, so the
72+
// textbook oo/oo cases at infinity -- ln(x) / x, x / e^x, ln(x) / sqrt(x) --
73+
// were left with nothing to catch them. The rule is only allowed to improve
74+
// on what is already there: sqrt(x) / sqrt(x + 1) merely turns into its own
75+
// reciprocal, and a NaN from it would claim the limit does not exist.
76+
if (ApplylHopitalRule(expr, x, dest) is { } lhopital && lhopital.Evaled != MathS.NaN)
77+
return lhopital;
78+
return atInfinity;
79+
}
6880
else if (expr.ComputeLimitDivideEtImpera(x, dest, ApproachFrom.Left) is { } fromLeft
6981
&& expr.ComputeLimitDivideEtImpera(x, dest, ApproachFrom.Right) is { } fromRight)
7082
{

‎Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs‎

Lines changed: 123 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77

88
using AngouriMath.Core.Multithreading;
99
using System;
10+
using System.Linq;
1011
using static AngouriMath.Entity;
1112

1213
namespace AngouriMath.Functions.Algebra
@@ -72,19 +73,138 @@ private static bool DivergesInMagnitude(Entity power, Variable x, Entity dest)
7273
private static bool IsFiniteNode(Entity expr)
7374
=> !IsInfiniteNode(expr) && expr != MathS.NaN;
7475

76+
/// <summary>
77+
/// How deep the rule may be applied to one limit. Differentiating both parts does not
78+
/// always make the quotient simpler -- x / sqrt(x^2 + 1) turns into sqrt(x^2 + 1) / x,
79+
/// which turns back, so the two stay indeterminate forever -- and without a bound the
80+
/// recursion below would not end. The bound has to leave room for the degree of a
81+
/// polynomial: x^10 / e^x takes ten steps.
82+
/// </summary>
83+
private const int MaxlHopitalDepth = 16;
84+
85+
[ThreadStatic] private static int lHopitalDepth;
86+
87+
/// <summary>
88+
/// How many times in all the rule may differentiate a quotient while one limit is being
89+
/// computed. Bounding the depth alone does not bound the work: each step asks what the
90+
/// two parts of its quotient tend to, and those are limits the rule may be applied to
91+
/// in turn, so the steps fan out rather than follow one another, and sixteen deep is
92+
/// far more than sixteen of them. This is a backstop rather than the answer to that --
93+
/// the two conditions below are what keep the fan-out from starting.
94+
/// </summary>
95+
private const int MaxlHopitalApplications = 64;
96+
97+
[ThreadStatic] private static int lHopitalApplications;
98+
99+
[ThreadStatic] private static List<Entity>? lHopitalChain;
100+
101+
/// <summary>
102+
/// Whether the rule is already partway through this very quotient. Differentiating both
103+
/// parts of sqrt(x^2 - x) / x gives its reciprocal, and differentiating that gives the
104+
/// first back, so the rule can go round for as long as the bound on its depth allows.
105+
/// Each step of the way costs a simplification and two limits of its own, which is why
106+
/// bounding the depth is not by itself enough to stop that limit running for a minute.
107+
/// </summary>
108+
private static bool AlreadyBeingDifferentiated(Entity quotient)
109+
=> (lHopitalChain ??= new()).Contains(quotient);
110+
111+
/// <summary>
112+
/// Whether differentiating both parts has left a quotient bigger than the one it came
113+
/// from by more than the derivative of a power or a logarithm accounts for. Each step
114+
/// asks what the two parts of its quotient tend to, and a bigger quotient makes those
115+
/// two questions harder than the one being answered, so a step that grows is a step
116+
/// away from an answer rather than towards one. The room left over is measured: the
117+
/// steps that do reach an answer add two nodes, and the most any of them was seen to
118+
/// add is six, while one step of x^(3/2) * sqrt(1 + 1/x^2) / x^2 adds eighteen and
119+
/// leaves the rule tens of seconds of work that ends in nothing.
120+
/// </summary>
121+
private static bool GrewTooMuch(Entity quotient, Entity applied)
122+
=> applied.Nodes.Count() > quotient.Nodes.Count() + 8;
123+
124+
/// <summary>
125+
/// A product that has a reciprocal factor in it, rewritten as a quotient, or
126+
/// <see langword="null"/> if it has none. The rule below only reads quotients, so
127+
/// <c>lim x -&gt; +oo x^4 * e^(-x)</c> came out as NaN while the same limit written
128+
/// <c>x^4 / e^x</c> gave 0 --
129+
/// https://github.com/asc-community/AngouriMath/issues/596.
130+
/// </summary>
131+
private static Entity? AsQuotient(Entity expr)
132+
{
133+
if (expr is not Mulf)
134+
return null;
135+
Entity numerator = 1, denominator = 1;
136+
var reciprocal = false;
137+
foreach (var factor in Mulf.LinearChildren(expr))
138+
switch (factor)
139+
{
140+
case Powf(var @base, Real { IsNegative: true } power):
141+
denominator *= @base.Pow(-power);
142+
reciprocal = true;
143+
break;
144+
case Powf(var @base, Mulf(Real { IsNegative: true } coefficient, var rest)):
145+
denominator *= @base.Pow(-coefficient * rest);
146+
reciprocal = true;
147+
break;
148+
case Divf(var dividend, var divisor):
149+
numerator *= dividend;
150+
denominator *= divisor;
151+
reciprocal = true;
152+
break;
153+
default:
154+
numerator *= factor;
155+
break;
156+
}
157+
return reciprocal ? numerator / denominator : null;
158+
}
159+
75160
private static Entity? ApplylHopitalRule(Entity expr, Variable x, Entity dest)
76161
{
162+
if (lHopitalDepth == 0)
163+
lHopitalApplications = 0;
164+
if (lHopitalDepth >= MaxlHopitalDepth || lHopitalApplications >= MaxlHopitalApplications)
165+
return null;
166+
// Held for the whole of the rule and not just for the recursive call below, since
167+
// asking what the two parts tend to is itself a limit that the rule may be applied
168+
// to, and counting only the recursion left those two free to nest without a bound.
169+
lHopitalDepth++;
170+
try { return ApplylHopitalRuleImpl(expr, x, dest); }
171+
finally { lHopitalDepth--; }
172+
}
173+
174+
private static Entity? ApplylHopitalRuleImpl(Entity expr, Variable x, Entity dest)
175+
{
176+
if (expr is not Divf && AsQuotient(expr) is { } quotient)
177+
expr = quotient;
77178
if (expr is Divf(var num, var den))
78179
if (EvalAssumingContinuous(num.Limit(x, dest)) is var numLimit && EvalAssumingContinuous(den.Limit(x, dest)) is var denLimit)
79180
if (numLimit == 0 && denLimit == 0 ||
80181
IsInfiniteNode(numLimit) && IsInfiniteNode(denLimit))
81182
if (num is not Number && den is not Number)
82183
if (num.ContainsNode(x) && den.ContainsNode(x))
83184
{
84-
var applied = num.Differentiate(x) / den.Differentiate(x);
185+
// Simplified, because the shape of the quotient is what the
186+
// machinery below matches on and differentiation does not
187+
// produce it: d/dx ln(x)^2 is 2 * ln(x) * (1 / x), a product
188+
// with a quotient inside rather than the quotient
189+
// 2 * ln(x) / x, and only the second has a limit here.
190+
var applied = (num.Differentiate(x) / den.Differentiate(x)).Simplify();
191+
// The domain condition simplification leaves behind is what the
192+
// derivative of a logarithm or a root carries, and a Providedf is
193+
// not a continuous node, so the quotient would be turned away
194+
// unread. Limits already treat the expression as continuous, as
195+
// SimplifyAndComputeLimitToInfinity does with the same shape.
196+
while (applied is Providedf(var body, _)) applied = body;
197+
if (AlreadyBeingDifferentiated(applied) || GrewTooMuch(expr, applied))
198+
return null;
199+
lHopitalApplications++;
85200
MultithreadingFunctional.ExitIfCancelled();
86-
if (ComputeLimit(applied, x, dest) is { } resLim)
87-
return resLim;
201+
lHopitalChain!.Add(applied);
202+
try
203+
{
204+
if (ComputeLimit(applied, x, dest) is { } resLim)
205+
return resLim;
206+
}
207+
finally { lHopitalChain.RemoveAt(lHopitalChain.Count - 1); }
88208
}
89209
return null;
90210
}

‎Sources/Tests/UnitTests/Calculus/LimitAtInfinityTest.cs‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
//
2+
// Copyright (c) 2019-2022 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 System.Threading.Tasks;
10+
using AngouriMath;
11+
using AngouriMath.Extensions;
12+
using Xunit;
13+
14+
namespace AngouriMath.Tests.Calculus
15+
{
16+
/// <summary>
17+
/// l'Hopital's rule was only applied when the destination was finite, so the ordinary
18+
/// competing-growth limits at infinity came back unevaluated.
19+
/// </summary>
20+
public sealed class LimitAtInfinityTest
21+
{
22+
private static Entity Limit(string expression, string destination) =>
23+
expression.ToEntity().Limit("x", destination.ToEntity());
24+
25+
[Theory]
26+
[InlineData("ln(x) / x", "+oo", "0")]
27+
[InlineData("x / e ^ x", "+oo", "0")]
28+
[InlineData("x ^ 2 / e ^ x", "+oo", "0")]
29+
[InlineData("ln(x) / sqrt(x)", "+oo", "0")]
30+
[InlineData("e ^ x / (e ^ x + 1)", "+oo", "1")]
31+
[InlineData("ln(x) ^ 2 / x", "+oo", "0")]
32+
[InlineData("x / e ^ (2 * x)", "+oo", "0")]
33+
[InlineData("x / ln(x) ^ 2", "+oo", "+oo")]
34+
[InlineData("x ^ 3 / ln(x)", "+oo", "+oo")]
35+
[InlineData("ln(ln(x)) / ln(x)", "+oo", "0")]
36+
[InlineData("sqrt(x) / sqrt(x + 1)", "+oo", "1")]
37+
[InlineData("x ^ 10 / e ^ x", "+oo", "0")]
38+
public void CompetingGrowthAtInfinity(string expression, string destination, string expected) =>
39+
Assert.Equal(expected.ToEntity(), Limit(expression, destination).Simplify());
40+
41+
// https://github.com/asc-community/AngouriMath/issues/596. The same limits written as
42+
// products with a reciprocal factor. Only the shape differs, and only the quotient
43+
// shape was read.
44+
[Theory]
45+
[InlineData("(x ^ (5 - 1)) * e ^ (-x)", "+oo", "0")]
46+
[InlineData("x ^ 4 * e ^ (-x)", "+oo", "0")]
47+
[InlineData("x * e ^ (-x)", "+oo", "0")]
48+
[InlineData("ln(x) * x ^ (-1)", "+oo", "0")]
49+
[InlineData("e ^ (-2 * x) * x ^ 3", "+oo", "0")]
50+
public void AProductWithAReciprocalFactorIsAQuotient(string expression, string destination, string expected) =>
51+
Assert.Equal(expected.ToEntity(), Limit(expression, destination).Simplify());
52+
53+
// The forms that already had an answer must keep it, in particular the ones that do
54+
// not go through the rule at all.
55+
[Theory]
56+
[InlineData("x ^ 2 / (x ^ 2 + 1)", "+oo", "1")]
57+
[InlineData("(2 * x + 1) / (x - 3)", "+oo", "2")]
58+
[InlineData("1 / x", "+oo", "0")]
59+
[InlineData("x ^ 3", "+oo", "+oo")]
60+
[InlineData("sin(x) / x", "0", "1")]
61+
[InlineData("(1 + x) ^ (1/x)", "0", "e")]
62+
[InlineData("(x ^ 2 - 1) / (x - 1)", "1", "2")]
63+
[InlineData("sqrt(x + 1) / sqrt(x)", "+oo", "1")]
64+
[InlineData("sqrt(x + 5) / sqrt(x + 1)", "+oo", "1")]
65+
public void EstablishedLimitsAreUnaffected(string expression, string destination, string expected) =>
66+
Assert.Equal(expected.ToEntity(), Limit(expression, destination).Simplify());
67+
68+
/// <summary>
69+
/// Differentiating both parts of x / sqrt(x^2 + 1) gives its own reciprocal, so the rule
70+
/// never settles and has to be stopped; x^20 / e^x would settle but only after more
71+
/// steps than the bound allows. Either way the limit is left unevaluated, which is
72+
/// honest. What must not happen is a hang, and the answer must not become NaN either,
73+
/// since NaN asserts that the limit does not exist.
74+
/// </summary>
75+
/// <remarks>
76+
/// The last two are the forms that showed a bound on the number of steps is not by
77+
/// itself a bound on the work. Each step asks what the two parts of its quotient tend
78+
/// to, and those are limits the rule may be applied to in turn, so sixteen steps deep
79+
/// is far more than sixteen steps: both of these ran for over twenty seconds before the
80+
/// rule was stopped from going round a cycle or handing on a quotient bigger than the
81+
/// one it came from.
82+
/// </remarks>
83+
[Theory]
84+
[InlineData("x / sqrt(x ^ 2 + 1)")]
85+
[InlineData("(x + sin(x)) / x")]
86+
[InlineData("x ^ 20 / e ^ x")]
87+
[InlineData("sqrt(x ^ 2 - x) / x")]
88+
[InlineData("x ^ (3/2) * sqrt(1 + 1 / x ^ 2) / x ^ 2")]
89+
public void AFormTheRuleCannotSettleTerminatesWithoutClaimingAnAnswer(string expression)
90+
{
91+
var task = Task.Run(() => Limit(expression, "+oo"));
92+
Assert.True(task.Wait(TimeSpan.FromSeconds(30)), "the limit did not terminate");
93+
Assert.IsType<Entity.Limitf>(task.Result);
94+
}
95+
}
96+
}

0 commit comments

Comments
 (0)