First implementation of the logarithmic integral 'li(z)'
Fix doctests for 'li(z)'
Add 'li(x)' to the doc
Add rewrite from Ei(log(f(x))) to li(f(x))

For example we can now do:

In [8]: f = Function("f")

In [9]: Ei(log(f(z))).rewrite(li)
Out[9]: li(f(z))

In [10]: _.rewrite(Ei)
Out[10]: Ei(log(f(z)))
Add rewrite of 'li(z)' into hypergeometric form

In [3]: li(z).rewrite(hyper)
Out[3]: log(z)*hyper((1, 1), (2, 2), log(z)) - log(1/log(z))/2 + log(log(z))/2 + EulerGamma

In [4]: pprint(_)
                                 ⎛  1   ⎞
                              log⎜──────⎟
        ┌─  ⎛1, 1 │       ⎞      ⎝log(z)⎠   log(log(z))
log(z)⋅ ├─  ⎜     │ log(z)⎟ - ─────────── + ─────────── + EulerGamma
       2╵ 2 ⎝2, 2 │       ⎠        2             2

In [5]: hyperexpand(_)
Out[5]: ((-log(log(z)) + Ei(log(z)))/log(z) - EulerGamma/log(z))*log(z) - log(1/log(z))/2 + log(log(z))/2 + EulerGamma

In [6]: pprint(_)
                                                     ⎛  1   ⎞
                                                  log⎜──────⎟
⎛-log(log(z)) + Ei(log(z))   EulerGamma⎞             ⎝log(z)⎠   log(log(z))
⎜───────────────────────── - ──────────⎟⋅log(z) - ─────────── + ─────────── + EulerGamma
⎝          log(z)              log(z)  ⎠               2             2

In [7]: simplify(_)
Out[7]: -log(1/log(z))/2 - log(log(z))/2 + Ei(log(z))

In [8]: pprint(_)
     ⎛  1   ⎞
  log⎜──────⎟
     ⎝log(z)⎠   log(log(z))
- ─────────── - ─────────── + Ei(log(z))
       2             2

In [8]: _.rewrite(li)
Out[8]: -log(1/log(z))/2 - log(log(z))/2 + li(z)

In [9]: pprint(_)
     ⎛  1   ⎞
  log⎜──────⎟
     ⎝log(z)⎠   log(log(z))
- ─────────── - ─────────── + li(z)
       2             2

This is (maybe modulo a mess with branch cuts) a pure 'li(z)'.
Add rewrite of 'li(z)' in terms of Meijer G functions

In [3]: li(z).rewrite(meijerg)
Out[3]: -log(1/log(z))/2 - log(-log(z)) + log(log(z))/2 - meijerg(((), (1,)), ((0, 0), ()), -log(z))

In [4]: pprint(_)
     ⎛  1   ⎞
  log⎜──────⎟
     ⎝log(z)⎠                  log(log(z))   ╭─╮2, 0 ⎛      1 │        ⎞
- ─────────── - log(-log(z)) + ─────────── - │╶┐     ⎜        │ -log(z)⎟
       2                            2        ╰─╯1, 2 ⎝0, 0    │        ⎠

In [5]: hyperexpand(_)
Out[5]: ((-log(log(z)*exp_polar(2*I*pi)) + Ei(log(z)) + 2*I*pi)/log(z) - EulerGamma/log(z))*log(z) + log(log(z)*exp_polar(I*pi)) - log(1/log(z))/2 - log(-log(z)) + log(log(z))/2 + EulerGamma

In [6]: pprint(_)
                                                                                        ⎛  1   ⎞
⎛     ⎛        2⋅ⅈ⋅π⎞                                  ⎞                             log⎜──────⎟
⎜- log⎝log(z)⋅ℯ     ⎠ + Ei(log(z)) + 2⋅ⅈ⋅π   EulerGamma⎟             ⎛        ⅈ⋅π⎞      ⎝log(z)⎠                  log(log(z))
⎜───────────────────────────────────────── - ──────────⎟⋅log(z) + log⎝log(z)⋅ℯ   ⎠ - ─────────── - log(-log(z)) + ─────────── + EulerGamma
⎝                  log(z)                      log(z)  ⎠                                  2                            2

In [7]: simplify(_)
Out[7]: -log(1/log(z))/2 - log(-log(z)) + log(log(z))/2 + Ei(log(z)) + I*pi

In [8]: pprint(_)
     ⎛  1   ⎞
  log⎜──────⎟
     ⎝log(z)⎠                  log(log(z))
- ─────────── - log(-log(z)) + ─────────── + Ei(log(z)) + ⅈ⋅π
       2                            2

In [9]: _.rewrite(li)
Out[9]: -log(1/log(z))/2 - log(-log(z)) + log(log(z))/2 + li(z) + I*pi

In [10]: pprint(_)
     ⎛  1   ⎞
  log⎜──────⎟
     ⎝log(z)⎠                  log(log(z))
- ─────────── - log(-log(z)) + ─────────── + li(z) + ⅈ⋅π
       2                            2

We have to see how we can get rid of those nested 'log's in front ...
Helper function '_eis(z)' to make 'Ei(z)' and 'li(z)' tractable for the Gruntz algorithm.
Use '_eis(z)' for computing limits of 'li(z)'

For example we can compute that:

In [11]: gruntz(1/li(z), z, oo)
Out[11]: 0
Use '_eis(z)' for computing limits of 'Ei(z)'

For example we can compute that:

In [4]: gruntz(Ei(z), z, oo)
Out[4]: oo

In [5]: gruntz(Ei(z-Ei(z)), z, oo)
Out[5]: 0

This one is obvious (IF we would use the relation between 'Ei' <-> 'li')

In [6]: gruntz(Ei(log(z))/(2*li(z)), z, oo)
Out[6]: 1/2
Fix args test failure
Tests for the 'li(z)' function
Tests for the 'Ei(z)' function rewriting
Tests for the '_eis(z)' helper function
Add a limit of 'Ei' which was wrong before adding '_eis'
Implementation of 'Li(z)'
Test args of 'Li(z)'
Improve doc of 'Li(z)' and add differentiation and numerical evaluation
Fix review comment
Import mpmath from sympy
Add tests for the 'Li(z)' function