Opened 14 years ago

Closed 10 years ago

# Make Li symbolic and work with complex input

Reported by: Owned by: William Stein Gary Furnish major sage-5.3 symbolics beginner, Li, log, integral, symbolics, calculus myurko, Benjamin Jones sage-5.3.rc0 Martin Cross Mike Hansen, Karl-Dieter Crisman, Burcin Erocal, Benjamin Jones N/A #11143

Make Li symbolic and work with complex input

Just use mpmath and the ideas from #11143. Probably will have to do a little work to keep the doctests from earlier, maybe deprecate a keyword or two related to precision.

Here is some example code from M. Yurko that explains how to do this. I think something based on this should be put into the Li function itself.

O.K. I defined li(x) as follows:

if z in RR and z >= 2: #check if real number greater than 2
return Li(z) +
in SAGE def
elif z == 1:
return -infinity
else: #mode for complex and below 2 from incomplete gamma
z = CDF(z)
return -gamma_inc(0,-log(z)) + (log(log(z))-log(1/log(z)))/2-
log(-log(z))

The first part uses SAGE's built in Li(x) but adjusts for the offset.
The second part should be self explanatory. The third part uses a
formula involving the incomplete gamma function which I found on the
Wolfram Functions website. On testing different values with an
external calculator,  the third statement appears to only be valid for
negative reals and complex numbers. This leaves the interval [0,2)
undefined. Please note that I have no background in complex analysis
and that my above statements about domain are only based upon
experimentation.


---

Apply trac_3401.v2.patch to the Sage library.

### comment:1 Changed 14 years ago by William Stein

No version

def li(z): #def log integral for real and complex variables
if z in RR and z >= 2: #check if real number greater than 2
return Li(z) +
in SAGE def
elif z == 0:
return 0
elif z > 1 and z < 2:
return Ei(log(z))
elif z == 1:
return -infinity
elif z > 0 and z < 1:
return
else: #mode for complex and below 2 from incomplete gamma
z = CDF(z)
return -gamma_inc(0,-log(z)) + (log(log(z))-log(1/log(z)))/2-
log(-log(z))



### comment:2 Changed 14 years ago by Alex Ghitza

Type: defect → enhancement

### comment:3 Changed 13 years ago by myurko

Status: new → needs_review

Sorry in advance to the reviewer and release manager, but I couldn't figure out how to flatten the patches without applying them.

### comment:5 Changed 13 years ago by Mike Hansen

Authors: → Michael Yurko → Mike Hansen

### comment:6 Changed 13 years ago by Mike Hansen

I've added a patch which folds the above patches together and deprecates the eps_rel and err_bound parameters so that code that uses them won't "break", but will get a deprecation warning.

I'm okay with myurko's changes so if someone could sign off on the deprecation warning, that'd be great.

### comment:7 Changed 13 years ago by Karl-Dieter Crisman

Status: needs_review → needs_work

Overall looks good, but there should be at least one doctest for the new DeprecationWarnings? (I think this was agreed upon somewhere on sage-devel), and there should also be documentation that this actually fulfills the ticket - namely, to extend Li to complex input! It certainly does, but I have no idea if the output is correct (I assume it is):

sage: Li(1+i)
-0.431252110896297 + 2.05958421419258*I
sage: Li(2+i)
0.366095261900308 + 1.22470693841030*I
sage: Li(2+2*i)
0.875423840014232 + 1.96947430597102*I
sage: Li(-2-2*i)
-0.333825651054542 - 3.94714365810975*I
sage: Li(-8)
-1.74509249432858 + 5.26897573517771*I
sage: Li(-10)
-2.04384864349662 + 5.69678038115052*I
sage: Li(-100)
-15.9214591889007 + 17.3366538615045*I


Something like that should be added.

### comment:8 Changed 13 years ago by Mike Hansen

Report Upstream: → N/A needs_work → needs_review

I've put up a new patch which address the above concerns.

### comment:9 Changed 13 years ago by Karl-Dieter Crisman

Reviewers: Mike Hansen → Mike Hansen, Karl-Dieter Crisman needs_review → positive_review

Looks good - sometimes slower, sometimes faster, but it's much better to have the complex functionality than worry about the rest. I removed an auxiliary function which only existed to allow the previous implementation. Positive review!

### comment:10 Changed 13 years ago by Burcin Erocal

Reviewers: Mike Hansen, Karl-Dieter Crisman → Mike Hansen, Karl-Dieter Crisman, Burcin Erocal positive_review → needs_work

Sorry to come in this late to the discussion, but this needs more work.

The prec argument to symbolic functions is deprecated, adding it to Li now doesn't make sense.

sage: gamma(10,prec=100)
.../_home_burcin__sage_init_sage_0.py:1: DeprecationWarning: The prec keyword argument is deprecated. Explicitly set the precision of the input, for example gamma(RealField(300)(1)), or use the prec argument to .n() for exact inputs, e.g., gamma(1).n(300), instead.
# -*- coding: utf-8 -*-
362880.00000000000000000000000


You can get the precision from the argument provided by the user. If the user needs a higher precision, they should explicitly convert the argument to a higher precision, for example by using RealFiel(300)(val).

We should also start converting these to proper symbolic functions that remain symbolic on exact input, but that can be left to another ticket.

### comment:11 Changed 13 years ago by myurko

What is wrong with the prec argument? By default it is left as None and will get the precision from the argument as you said.

### comment:13 Changed 11 years ago by Karl-Dieter Crisman

Possibly related to #11143.

### comment:14 Changed 11 years ago by Benjamin Jones

In #11143 a fully symbolic function is defined for li(x) called exp_integral_li, the non-offset logarithmic integral. It would be very easy to add Li by simply returning exp_integral_li(x) - exp_integral_li(2). On the other hand, adding a symbolic version of Li would be equally easy by copying the definition of exp_integral_li and making one simple change.

If so, what would a good name for the offset Li be? Maybe exp_integral_li_offset?

Another thought is that in #11143, we could add an optional parameter offset to the init method for exp_integral_li which would return Li instead of li. The eval and evalf methods could be bound in the init call to return the right values and the derivative is obviously the same for both.

Either of these solutions could be put into #11143 without much effort and that would take care of the issue in this ticket because evaluation at complex inputs is handled by mpmath for all the functions defined there.

### comment:15 Changed 11 years ago by Karl-Dieter Crisman

My understanding is that the offset Li is the same as li. But maybe I've missed something while looking into this - I'm not a special functions expert.

I think that as long as we have both of these, and not named super-crazily - such as just being named Li and li - this would be fine. I think the parameter is not needed.

### comment:16 Changed 11 years ago by Leif Leonhardy

Just for the record:

sage: import mpmath
sage: mpmath.li(1+i)
mpc(real='0.61391166922119556', imag='2.0595842141925775')
sage: mpmath.li(1+i, offset=True)
mpc(real='-0.43125211089629728', imag='2.0595842141925775')


But maybe I've missed something (tl;dr).

### comment:17 Changed 11 years ago by Karl-Dieter Crisman

I don't think you're missing anything. In #11143 I think Benjamin is using mpmath as much as possible (though we should be checking timings...). In principle, the hope is that #11143 will render this ticket obsolete, but I like to keep things complete for trollers :)

### comment:18 Changed 11 years ago by Karl-Dieter Crisman

Description: modified (diff) extend li to work with complex input → Make Li symbolic and work with complex input

I'm changing this to make the (offset) Li symbolic and to work with complex input. Simply using the ideas of #11143 should be sufficient.

### comment:19 Changed 11 years ago by Karl-Dieter Crisman

Description: modified (diff)

Symbolic Li

### comment:20 follow-up:  21 Changed 10 years ago by Martin George Cross

Dependencies: → #11143

I have created a symbolic Li patch on top of #11143 on sage-5.2.rc1 . This is my first go at a patch so no doubt will need a good scrubbing...

Please note doing this is a hobby for me and I have little or no time weekdays to do anything so my responses are likely to be slow.

### comment:21 in reply to:  20 Changed 10 years ago by Karl-Dieter Crisman

I have created a symbolic Li patch on top of #11143 on sage-5.2.rc1 . This is my first go at a patch so no doubt will need a good scrubbing...

That's okay, we all have to start somewhere!

Please note doing this is a hobby for me and I have little or no time weekdays to do anything so my responses are likely to be slow.

That's also very true for many of us. So we may also be slow to respond.

### comment:22 follow-up:  24 Changed 10 years ago by Benjamin Jones

Hi @martinx, your patch looks very good. I spotted a few whitespace issues to clean up (I'll post a reviewer patch to do that). I'm running full tests now, but I expect a positive review.

I wonder if we really need three aliases at the top level for this function. Having log_integral_eulerian in addition to Li and log_integral_offset seems excessive to me, but if that's a common name for the function I'm OK with it.

reviewer patch

### comment:23 Changed 10 years ago by Benjamin Jones

Reviewers: Mike Hansen, Karl-Dieter Crisman, Burcin Erocal → Mike Hansen, Karl-Dieter Crisman, Burcin Erocal, Benjamin Jones needs_work → needs_review

### comment:24 in reply to:  22 Changed 10 years ago by Leif Leonhardy

I wonder if we really need three aliases at the top level for this function. Having log_integral_eulerian in addition to Li and log_integral_offset seems excessive to me, but if that's a common name for the function I'm OK with it.

There's also "European Li" (for the offset one) IIRC... ;-)

### comment:25 Changed 10 years ago by Martin George Cross

I will go with whatever the sagemath intelligentsia thinks appropriate for the number and name of any aliases.

And I had better read up on coding conventions before my next efforts :)

### comment:26 follow-up:  32 Changed 10 years ago by Benjamin Jones

Don't worry about the coding conventions, some of them are unwritten and some of them are subtle. I found a few doctest errors after running make ptestlong with the patches applied. These should be simple to fix:

Change 9 to 10:

File "/home/jonesbe/sage/sage-5.2/devel/sage-11143/sage/misc/sagedoc.py", line 971:
sage: len(search_src('log', 'derivative', interact=False).splitlines()) < 9
Expected:
True
Got:
False


simple change, Li is now fully symbolic:

File "/home/jonesbe/sage/sage-5.2/devel/sage-11143/sage/functions/transcendental.py", lin
e 195:
sage: Li(100)
Expected:
29.080977804
Got:
-log_integral(2) + log_integral(100)


This one is more mysterious:

File "/home/jonesbe/sage/sage-5.2/devel/sage-11143/sage/symbolic/random_tests.py", line 2
36:
sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_eleme
nt) # random
Exception raised:
Traceback (most recent call last):      File "/home/jonesbe/sage/sage-5.2/local/bin/ncadoctest.py", line 1231, in run_one_test
self.run_one_example(test, example, filename, compileflags)
File "/home/jonesbe/sage/sage-5.2/local/bin/sagedoctest.py", line 38, in run_one_ex
ample
OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
File "/home/jonesbe/sage/sage-5.2/local/bin/ncadoctest.py", line 1172, in run_one_e
xample
compileflags, 1) in test.globs
File "<doctest __main__.example_5[4]>", line 1, in <module>
print "ignore this";  random_expr(Integer(50), nvars=Integer(3), coeff_generator=CDF.random_element) # random###line 236:
sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_element) # random

sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_element) # random
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 258, in random_expr
return random_expr_helper(size, internal, leaves, verbose)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 206, in random_expr_helper
children = [random_expr_helper(n+1, internal, leaves, verbose) for n in nodes_per_child]
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 209, in random_expr_helper
return r[1](*children)
File "function.pyx", line 432, in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:4941)
res = g_function_evalv(self._serial, vec, hold)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/integration/integral.py", line 173, in _eval_
return integrator(*args)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/integration/external.py", line 21, in maxima_integrator
result = maxima.sr_integral(expression, v, a, b)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/interfaces/maxima_lib.py", line 746, in sr_integral
raise error
RuntimeError: ECL says: Error executing code in Maxima: defint: upper limit of integration must be real; found
elliptic_eu(.18648175298340663*I-.7457199773032457,
coth(v3)*(v3*(.12348638361486497*I+.29875723285490263)
+v1*(.12348638361486497*I+.29875723285490263)
-sinh(v3^(.5481180571998028*I-.5534231539946481))))


### comment:27 Changed 10 years ago by Leif Leonhardy

I will go with whatever the sagemath intelligentsia thinks appropriate for the number and name of any aliases.

Well, with names put into the global namespace, the most important thing is that tab completion is likely to suggest you what you're looking for, i.e., the prefix of each name matters. So in this case I think a single instance of log_int* (in addition to [Ll]i*, maybe more) would be sufficent.

With Sage 5.3.beta0, I currently get:

sage: log<TAB>
log             log_b           log_gamma       log_text        logstr
log2            log_dvi         log_html        lognormvariate
sage: li<TAB>
lie                    limit                  linear_relation        list_composition       list_plot_semilogy
lie_console            line                   linear_transformation  list_plot
lift                   line2d                 lisp                   list_plot3d
lift_to_sl2z           line3d                 lisp_console           list_plot_loglog
sage: Li<TAB>
Li                         LinearCode                 LinearCodeFromVectorSpace
LiE                        LinearCodeFromCheckMatrix  Lisp
sage: euler<TAB>
euler_gamma             euler_phi               eulers_method_2x2
euler_number            eulers_method           eulers_method_2x2_plot


(Of course also the docstring for e.g. li, perhaps that of Ei, too, should refer to Li and vice versa.)

### comment:28 Changed 10 years ago by Leif Leonhardy

(Also on top of the reviewer patch: ;-) )

s/eulerian/Eulerian/

s/for x > 1/for x \geq 2/

### comment:29 Changed 10 years ago by Leif Leonhardy

Nice theorem:

\exists x : \pi(x) > \operatorname{Li}(z)

(s/z/x/)

s/

However it is a theorem that there are very large, (e.g., around 10^{316}) values of x


/

However, it is a theorem that there are very large values of x (e.g., around 10^{316})


/

### comment:30 follow-up:  31 Changed 10 years ago by Leif Leonhardy

More nitpicking: s/"polylog()"/polylog()/ and/or make (it) a cross-reference (:func:polylog -- not sure whether it has to be :class:sage.functions.log.Function_polylog).

### comment:31 in reply to:  30 ; follow-up:  33 Changed 10 years ago by Leif Leonhardy

... not sure whether it has to be :class:sage.functions.log.Function_polylog).

... or rather :class:polylog <sage.functions.log.Function_polylog> or something like that.

### comment:32 in reply to:  26 Changed 10 years ago by Martin George Cross

Don't worry about the coding conventions, some of them are unwritten and some of them are subtle.

and I am a slow learner anyway :)

I found a few doctest errors after running make ptestlong with the patches applied. These should be simple to fix:

Change 9 to 10:

Agreed.

simple change, Li is now fully symbolic:

File "/home/jonesbe/sage/sage-5.2/devel/sage-11143/sage/functions/transcendental.py", lin
e 195:
sage: Li(100)
Expected:
29.080977804
Got:
-log_integral(2) + log_integral(100)


The function can be removed since it was just a convenience one to support the original li and Li.

This one is more mysterious:

File "/home/jonesbe/sage/sage-5.2/devel/sage-11143/sage/symbolic/random_tests.py", line 2
36:
sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_eleme
nt) # random
Exception raised:
Traceback (most recent call last):      File "/home/jonesbe/sage/sage-5.2/local/bin/ncadoctest.py", line 1231, in run_one_test
self.run_one_example(test, example, filename, compileflags)
File "/home/jonesbe/sage/sage-5.2/local/bin/sagedoctest.py", line 38, in run_one_ex
ample
OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
File "/home/jonesbe/sage/sage-5.2/local/bin/ncadoctest.py", line 1172, in run_one_e
xample
compileflags, 1) in test.globs
File "<doctest __main__.example_5[4]>", line 1, in <module>
print "ignore this";  random_expr(Integer(50), nvars=Integer(3), coeff_generator=CDF.random_element) # random###line 236:
sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_element) # random

sage: print "ignore this";  random_expr(50, nvars=3, coeff_generator=CDF.random_element) # random
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 258, in random_expr
return random_expr_helper(size, internal, leaves, verbose)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 206, in random_expr_helper
children = [random_expr_helper(n+1, internal, leaves, verbose) for n in nodes_per_child]
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/random_tests.py", line 209, in random_expr_helper
return r[1](*children)
File "function.pyx", line 432, in sage.symbolic.function.Function.__call__ (sage/symbolic/function.cpp:4941)
res = g_function_evalv(self._serial, vec, hold)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/integration/integral.py", line 173, in _eval_
return integrator(*args)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/symbolic/integration/external.py", line 21, in maxima_integrator
result = maxima.sr_integral(expression, v, a, b)
File "/home/jonesbe/sage/sage-5.2/local/lib/python/site-packages/sage/interfaces/maxima_lib.py", line 746, in sr_integral
raise error
RuntimeError: ECL says: Error executing code in Maxima: defint: upper limit of integration must be real; found
elliptic_eu(.18648175298340663*I-.7457199773032457,
coth(v3)*(v3*(.12348638361486497*I+.29875723285490263)
+v1*(.12348638361486497*I+.29875723285490263)
-sinh(v3^(.5481180571998028*I-.5534231539946481))))


The function randomly selects from a list of all Pynac functions, to which is now added Li, and so now the functions chosen have changed. The docs state that it will often raise an error because it tries to create an erroneous expression. In this case it is trying to pass a complex expression to Maxima. Trial and error and got a return result setting the seed to 1.

When I have worked out how to merge patches I'll post a revised patch for review that addresses all comments made so far.

### comment:33 in reply to:  31 Changed 10 years ago by Martin George Cross

... not sure whether it has to be :class:sage.functions.log.Function_polylog).

... or rather :class:polylog <sage.functions.log.Function_polylog> or something like that.

Both seem to be allowed: http://www.sagemath.org/doc/developer/sage_manuals.html#chapter-sage-manuals-links

### Changed 10 years ago by Martin George Cross

Revised Li including reviewer patch

### comment:34 Changed 10 years ago by Martin George Cross

I think v2 addresses all comments made so far and includes reviewer patch. ptestlong passes I think; I had some issues with sagedoc.py but that now passes all tests. Too tired to run ptestlong yet again...........

### comment:35 Changed 10 years ago by Benjamin Jones

There are lots of lines touched in sage/functions/transcendental.py that don't seem to show any difference. What's going on there? Have invisible non-whitespace characters been deleted?

I looked over the command line and HTML docs for Li and they look good. Running tests now.

### comment:36 follow-up:  37 Changed 10 years ago by Benjamin Jones

Authors: Michael Yurko → martinx

I see, it's just whitespace at the beginnings of lines. I guess that doesn't get highlighted by trac when you view the patch. Anyway, I think it's a good thing to clean up trailing whitespace, but you'll have to watch out that touching so many lines of the file doesn't cause conflicts with other patches that modify sage/functions/transcendental.py. In this case I think it's probably OK, I don't know of any other currently pending positively reviewed patches that touch that file.

If leif is happy with the patch, I'll give it a positive review.

### comment:37 in reply to:  36 Changed 10 years ago by Martin George Cross

I see, it's just whitespace at the beginnings of lines. I guess that doesn't get highlighted by trac when you view the patch. Anyway, I think it's a good thing to clean up trailing whitespace, but you'll have to watch out that touching so many lines of the file doesn't cause conflicts with other patches that modify sage/functions/transcendental.py. In this case I think it's probably OK, I don't know of any other currently pending positively reviewed patches that touch that file.

I got carried away with strip trailing whitespace command in Geany, in response to the previous review comments. Will try to be more restrained next time ;-) Martin

### comment:38 Changed 10 years ago by Benjamin Jones

Patches apply cleanly on top of those at #11143 on top of sage-5.3.beta0. All long tests pass. I think this is ready to go.