Gruntz Demo

[anutosh491]

This is a draft PR trying to replicate the gruntz algorithm and a few test cases for the same.

[anutosh491]

cc @certik

[anutosh491]

We can demonstrate the sample case that sympy use in their documentation which is limit(sin(x)/x, x, 0) or something even simpler limit(sin(x), x, 0) . The breakdown for the following would be this

In [1]: limit(sin(x), x, 0)
limitinf(sin(1/_x), _x) = 0
+-mrv_leadterm(sin(1/x), x) = (1, 1)
| +-mrv(sin(1/_x), _x) = set([_x])
| +-rewrite(sin(1/_x), set([exp(_x)]), _x, _w) = (sin(_w), _x)
+-sign(1, _x) = 1
+-limitinf(1, _x) = 0

So limitinf technically uses mrv_leadterm & sign.

i) mrv_leadterm basically returns the the coeff & exponent of the most rapidly varying leadterm of the expression. In this case sin(x) (or rather sin(_w)) itself is the most rapidly varying terms from our input expression sin(x). Hence we return (1, 1) which basically comes from 1*x **1 which is the first term of the series expansion of sin(x)

ii) We then use sign on the generated exponent and the input variable x . So we call sign(1, x) and the logic here is pretty simple

        e >  0 for x sufficiently large ...  1
        e == 0 for x sufficiently large ...  0
        e <  0 for x sufficiently large ... -1

Now based on the value returned through sign, we return the output

    sig = sign(e0, x)
   if sig == 1:
       return S.Zero  # e0>0: lim f = 0
   elif sig == -1:  # e0<0: lim f = +-oo (the sign depends on the sign of c0)
       if c0.match(I*Wild("a", exclude=[I])):
           return c0*oo
       s = sign(c0, x)
       # the leading term shouldn't be 0:
       if s == 0:
           raise ValueError("Leading term should not be 0")
       return s*oo
   elif sig == 0:
       if c0 == old:
           c0 = c0.cancel()
       return limitinf(c0, x)  # e0=0: lim f = lim c0
   else:
       raise ValueError("{} could not be evaluated".format(sig))

[anutosh491]

Internally mrv_leadterm depends on mrv , rewrite and leadterm

mrv_leadterm(sin(1/x), x) 
1) exps = mrv(sin(1/x), x) = sin(1/x)
2) 
    w = Dummy("w", positive=True)
    f, logw = rewrite(exps, Omega, x, w)

This gives us f which would be sin(_w) 3) Then we calculate the leadterm for sin(_w) which is x -> (1, 1)

[anutosh491]

A slightly more involved example would be this

In [1]: limit(sin(x)/x, x, 0)
limitinf(_x*sin(1/_x), _x) = 1
+-mrv_leadterm(_x*sin(1/_x), _x) = (1, 0)
| +-mrv(_x*sin(1/_x), _x) = set([_x])
| | +-mrv(_x, _x) = set([_x])
| | +-mrv(sin(1/_x), _x) = set([_x])
| |   +-mrv(1/_x, _x) = set([_x])
| |     +-mrv(_x, _x) = set([_x])
| +-mrv_leadterm(exp(_x)*sin(exp(-_x)), _x, set([exp(_x)])) = (1, 0)
|   +-rewrite(exp(_x)*sin(exp(-_x)), set([exp(_x)]), _x, _w) = (1/_w*sin(_w), -_x)
|     +-sign(_x, _x) = 1
|     +-mrv_leadterm(1, _x) = (1, 0)
+-sign(0, _x) = 0
+-limitinf(1, _x) = 1

[anutosh491]

The current implementation can be simplified quite a bit

1. limitinf technically depends on

• mrv_leadterm : To calculate leadterm of most rapidly varying subexpression
• sign : To compute limitinf based on exponent of calculated leadterm

2. mrv_leadterm technically depends on

• rewrite : To basically rewrite the expression in terms of the most rapidly varying (mrv) expression
• " Any general leadterm computing function": Like we have the leadterm from sympy which is actually the function that is responsible for generating the leading term's coeff and exponent

3. rewrite technically depends on

• mrv : To calculate the most rapidly varying term.

So we have something like this

I shall try simplifying our algorithm in the next commit.

[certik]

Perfect. Yes, let's do the sin(x)/x, and you can take the running example in sympy and just simplify it to the bare minimum. Then we'll port the simplified code to LPython.

[anutosh491]

The latest commit fixes some import errors and any issues with the code and now we can get results for our simplified gruntz algorithm through sympy. I've added a test for the same saying

# tests
x = Symbol('x')
ans = gruntz(sin(x)/x, x, 0)
print(ans)

I could maybe add another file that is relevant only for the sin(x)/x case. We could keep this general file for a broader use case whenever required.

[anutosh491]

The recent commits simplifies the previous version of gruntz by quite a bit. For our use case i.e. sin(x)/x

1. We don't have to worry about ideas related to tractable/intractable functions That would only come into picture we would be dealing with inverse trig or inverse hyperbolic functions etc

2. We don't need to worry about the logw parameter so that has been removed. That would only come into picture if we would be dealing with nested logs, exponentials for example log(log(1/x)) or exp(exp(log(1/x))

[anutosh491]

Also talking about the rewrite function. The rewrite functions is cruicial and can't be simplified any further but it's actual requirements comes into picture when we have exp-log functions. For our use case we can simplify it if we want

So the mrv term would in almost all cases be a part of one of these classes const, log, linear (x) or exp . Now if we consider our test case or any case where we know the mrv is linear, we essentially know what the rewrite would be.

>>> x = Symbol('x')
>>> rewrite(x, x, w)
(1/_w, -x)
>>> rewrite(log(x), x, w)
(log(1/_w), -x)
>>> rewrite(x**2, x, w)
((1/_w)**2, -x)
>>> rewrite(exp(log(x)), x, w)
(1/_w, -x)
>>> rewrite(sin(x), x, w)
(sin(1/_w), -x)
>>> rewrite(sin(x)/x, x, w)
(sin(1/_w)/1/_w, -x)

All of these have x as their mrv and so the rewrite would just be substituing 1/x in place of x in the expression.

[anutosh491]

For the signinf function, we would need the mathematical function sign

def signinf(e, x):
    r"""
    Determine sign of the expression at the infinity.

    Returns
    =======

    {1, 0, -1}
        One or minus one, if `e > 0` or `e < 0` for `x` sufficiently
        large and zero if `e` is *constantly* zero for `x\to\infty`.

    """

    if not e.has(x):
        return sign(e)
    if e == x or (e.is_Pow and signinf(e.base, x) == 1):
        return S(1)     

I've raised a PR to symengine to add functionality for sign ( using basic_sign) here https://github.com/symengine/symengine/pull/2016

[anutosh491]

The limitinf function would require support for the subs function

def gruntz(e, z, z0, dir="+"):
    if str(dir) == "-":
        e0 = e.subs(z, z0 - 1/z)
    elif str(dir) == "+":
        e0 = e.subs(z, z0 + 1/z)
    else:
        raise NotImplementedError("dir must be '+' or '-'")

    r = limitinf(e0, z)
    return r 

And technically the rewrite function would also be using subs or xreplace function. Hence I've raised a PR adding support for the subs attrbute here https://github.com/lcompilers/lpython/pull/2695

[anutosh491]

Listing all blockers below.