Invalid limit of floating point matrix power

[asmeurer]

This is from https://stackoverflow.com/questions/39736887/use-python-to-calculate-pn-while-n-approach-oo-and-p-is-a-matrix

In [1]: p = Matrix([[1./2,1./4,1./4],[1./2,0,1./2],[1./4,0,3./4]])

In [2]: (p**k)
Out[2]:
⎡                                    k                                      k                        k                                           k
⎢  0.2165423646591⋅-0.154508497187474  + 0.419821271704536⋅0.404508497187474  + 0.363636363636364⋅1.0      - 0.350372906022699⋅-0.154508497187474  + 0.25946
⎢
⎢                                      k                                      k                        k                                       k
⎢- 0.507064433090879⋅-0.154508497187474  + 0.143428069454515⋅0.404508497187474  + 0.363636363636364⋅1.0    0.820447487227238⋅-0.154508497187474  + 0.0886434
⎢
⎢                                       k                                      k                        k                                        k
⎣- 0.0598508375909206⋅-0.154508497187474  - 0.303785526045443⋅0.404508497187474  + 0.363636363636364⋅1.0    0.0968406894772594⋅-0.154508497187474  - 0.18774

                            k                         k                                       k                                      k
3815113608⋅0.404508497187474  + 0.0909090909090909⋅1.0    0.133830541363598⋅-0.154508497187474  - 0.679285086818144⋅0.404508497187474  + 0.545454545454546⋅1

                           k                         k                                         k                                      k
218636708⋅0.404508497187474  + 0.0909090909090909⋅1.0     - 0.31338305413636⋅-0.154508497187474  - 0.232071491318186⋅0.404508497187474  + 0.545454545454546⋅

                           k                         k                                          k                                      k
978038635⋅0.404508497187474  + 0.0909090909090909⋅1.0    - 0.0369898518863388⋅-0.154508497187474  + 0.491535306431793⋅0.404508497187474  + 0.545454545454546

  k  ⎤
.0   ⎥
     ⎥
   k ⎥
1.0  ⎥
     ⎥
    k⎥
⋅1.0 ⎦

In [3]: (p**k).applyfunc(lambda i: limit(i, k, oo))
Out[3]:
⎡∞  ∞  ∞⎤
⎢       ⎥
⎢∞  ∞  ∞⎥
⎢       ⎥
⎣∞  ∞  ∞⎦

I think the problem is the 1.0**k terms, which confuse limit.

[Upabjojr]

In [2]: ;var k
Out[2]: k

In [3]: limit(1.0**k, k, oo)
Out[3]: oo

Wrong, should be 1.

[Upabjojr]

In sympy.series.gruntz.py, function mrv(1.0**k, k) goes into infinite loop because of roundoff errors: log(1.0) ==> 2.36364252615315e-125.

Suggest to replace floating point numbers with rationals during the limit calculation.

[Upabjojr]

https://github.com/sympy/sympy/pull/11681

The problem is not just 1.0**k for k ===> oo.

[asmeurer]

I remember discussing this before. This expression is tricky because any deviation in the base changes the value of the limit.