sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.9k stars 4.42k forks source link

Gruntz sign function, imaginary number #23835

Open egolf-cs opened 2 years ago

egolf-cs commented 2 years ago

G = (-2*sqrt(5)*(-1 + sqrt(5))*(1/2 - sqrt(5)/2)**k/5 - 2*sqrt(5)*(1/2 + sqrt(5)/2)**k*(1 + sqrt(5))/5)/(-2*sqrt(5)*(-1 + sqrt(5))*(1/2 - sqrt(5)/2)**k/5 + (-1 + sqrt(5))*(1/2 - sqrt(5)/2)**k*(sqrt(5)/5 + 1) - 2*sqrt(5)*(1/2 + sqrt(5)/2)**k*(1 + sqrt(5))/5 + (-1 + sqrt(5)/5)*(1/2 + sqrt(5)/2)**k*(1 + sqrt(5)))

I'm trying to compute gruntz(G, sympify("k"), oo)---error + trace below. The result should be 1/GoldenRatio. What would need to change for this limit to go through? Is the issue that log(-1/2 + sqrt(5)/2) + I*pi is a complex number? I noticed a comment in the code:

    # For limits of complex functions, the algorithm would have to be
    # improved, or just find limits of Re and Im components separately.

Would finding these limits separately make the limit go through and how much work would it be to treat these components separately? Thanks!

Traceback (most recent call last):
  File "file.py", line 23, in <module>
    gruntz(sympify(Gexpr), sympify("k"), sympify(oo))
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/series/gruntz.py", line 709, in gruntz
    r = limitinf(e0, z)
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/core/cache.py", line 70, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/series/gruntz.py", line 452, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/core/cache.py", line 70, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/series/gruntz.py", line 545, in mrv_leadterm
    f, logw = rewrite(exps, Omega, x, w)
  File "/home/user/mambaforge/envs/sage/lib/python3.9/site-packages/sympy/series/gruntz.py", line 630, in rewrite
    raise NotImplementedError('Result depends on the sign of %s' % sig)
NotImplementedError: Result depends on the sign of sign(log(-1/2 + sqrt(5)/2) + I*pi)
jksuom commented 2 years ago

gruntz is basically designed to work with real exp-log functions, i.e., functions obtained by composing exp, log and rational functions with real coefficients, and taking real values in a semi-infinite interval (a, oo) on the real line. The implementation also works with complex multiples of such functions and some special functions like gamma that have a suitable asymptotic expansion.

The message Result depends on the sign of ,,, usually means that the expression is not in the domain of gruntz. In this case, it cannot handle the subexpression (1/2 - sqrt(5)/2)**k because it is not real as a function of k in any semi-infinite interval. It can be expressed as exp(log(1/2 - sqrt(5)/2)*k) but the logarithm is not real.

Would finding these limits separately make the limit go through

The real and imaginary parts of exp(log(1/2 - sqrt(5)/2)*k) are not exp-log functions even if they are real. A basic property of exp-log functions is that they (or actually their germs) form a field, in particular they can be divided. That is possible because an exp-log function has no zeros in some interval (a, oo) so division does not introduce singularities. The real and imaginary parts of exp(log(1/2 - sqrt(5)/2)*k) change sign repeatedly in each interval (a, oo) so I don't think gruntz could work on them.

egolf-cs commented 2 years ago

Thank you for your detailed response. I don’t want to go off topic too much, but are you aware of any extensions to Gruntz’s work which handle strictly more general limits? It seems like his algorithm is the state of the art, but it’s also 20+ years old.

oscarbenjamin commented 3 months ago

After a little rearrangement the limit here becomes trivial:

In [27]: G = ratsimp(G)

In [28]: G
Out[28]: 
                             k              k                      
                  10⋅(1 - √5)  + 10⋅(1 + √5)                      1
─────────────────────────────────────────────────────────────── + ─
               k              k              k                k   3
- 9⋅√5⋅(1 - √5)  + 15⋅(1 - √5)  + 15⋅(1 + √5)  + 9⋅√5⋅(1 + √5)     

In [29]: G2 = G.xreplace({1 - sqrt(5): x, 1 + sqrt(5): y}).xreplace({x: (1 - sqrt(5))/(1 + sqrt(5)),
    ...: y:1})

In [30]: G2
Out[30]: 
                        k                      
                ⎛1 - √5⎞                       
             10⋅⎜──────⎟  + 10                 
                ⎝1 + √5⎠                      1
─────────────────────────────────────────── + ─
               k              k               3
       ⎛1 - √5⎞       ⎛1 - √5⎞                 
- 9⋅√5⋅⎜──────⎟  + 15⋅⎜──────⎟  + 15 + 9⋅√5    
       ⎝1 + √5⎠       ⎝1 + √5⎠                 

In [31]: G2.subs(k, oo)
Out[31]: 
   10       1
───────── + ─
15 + 9⋅√5   3

In [32]: simplify(_)
Out[32]: 
  1   √5
- ─ + ──
  2   2 

It doesn't seem that any current limit finding routines can handle this even after rearrangement though. There should probably be better handling of rational functions like this somehow.

arnabnandikgp commented 3 months ago

ohk, I din't look at this workaround previously. Maybe need to look at how such cases can be handled.