Opened 12 years ago

Closed 7 years ago

#3401 closed enhancement (fixed)

Make Li symbolic and work with complex input

Reported by: was Owned by: gfurnish
Priority: major Milestone: sage-5.3
Component: symbolics Keywords: beginner, Li, log, integral, symbolics, calculus
Cc: myurko, benjaminfjones Merged in: sage-5.3.rc0
Authors: Martin Cross Reviewers: Mike Hansen, Karl-Dieter Crisman, Burcin Erocal, Benjamin Jones
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: #11143 Stopgaps:

Description (last modified by benjaminfjones)

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:

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) +
1.045163780117492784844588889194613136522615578151 #adjust for offset
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.

Attachments (8)

Li(x).patch (3.3 KB) - added by myurko 10 years ago.
Lix_2.patch (651 bytes) - added by myurko 10 years ago.
Lix_3.patch (599 bytes) - added by myurko 10 years ago.
trac_3401.patch (5.1 KB) - added by mhansen 10 years ago.
trac_3401-reviewed.patch (5.5 KB) - added by kcrisman 10 years ago.
trac_3401_Li.patch (9.8 KB) - added by martinx 7 years ago.
Symbolic Li
trac_3401_review.patch (2.2 KB) - added by benjaminfjones 7 years ago.
reviewer patch
trac_3401.v2.patch (21.7 KB) - added by martinx 7 years ago.
Revised Li including reviewer patch

Download all attachments as: .zip

Change History (49)

comment:1 Changed 12 years ago by was

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) +
1.045163780117492784844588889194613136522615578151 #adjust for offset
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 11 years ago by AlexGhitza

  • Type changed from defect to enhancement

Changed 10 years ago by myurko

Changed 10 years ago by myurko

Changed 10 years ago by myurko

comment:3 Changed 10 years ago by myurko

  • Status changed from new to 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:4 Changed 10 years ago by myurko

  • Cc myurko added

comment:5 Changed 10 years ago by mhansen

  • Authors set to Michael Yurko
  • Reviewers set to Mike Hansen

comment:6 Changed 10 years ago by mhansen

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 10 years ago by kcrisman

  • Status changed from needs_review to 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.

Changed 10 years ago by mhansen

comment:8 Changed 10 years ago by mhansen

  • Report Upstream set to N/A
  • Status changed from needs_work to needs_review

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

Changed 10 years ago by kcrisman

comment:9 Changed 10 years ago by kcrisman

  • Reviewers changed from Mike Hansen to Mike Hansen, Karl-Dieter Crisman
  • Status changed from needs_review to 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 10 years ago by burcin

  • Reviewers changed from Mike Hansen, Karl-Dieter Crisman to Mike Hansen, Karl-Dieter Crisman, Burcin Erocal
  • Status changed from positive_review to 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 10 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:12 Changed 9 years ago by burcin

  • Keywords beginner added

comment:13 Changed 8 years ago by kcrisman

  • Cc benjaminfjones added

Possibly related to #11143.

comment:14 Changed 8 years ago by benjaminfjones

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 8 years ago by kcrisman

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 8 years ago by leif

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 8 years ago by kcrisman

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 8 years ago by kcrisman

  • Description modified (diff)
  • Summary changed from extend li to work with complex input to 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 8 years ago by kcrisman

  • Description modified (diff)

Changed 7 years ago by martinx

Symbolic Li

comment:20 follow-up: Changed 7 years ago by martinx

  • Dependencies set to #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 7 years ago by kcrisman

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: Changed 7 years ago by benjaminfjones

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.

Changed 7 years ago by benjaminfjones

reviewer patch

comment:23 Changed 7 years ago by benjaminfjones

  • Reviewers changed from Mike Hansen, Karl-Dieter Crisman, Burcin Erocal to Mike Hansen, Karl-Dieter Crisman, Burcin Erocal, Benjamin Jones
  • Status changed from needs_work to needs_review

comment:24 in reply to: ↑ 22 Changed 7 years ago by leif

Replying to benjaminfjones:

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 7 years ago by martinx

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: Changed 7 years ago by benjaminfjones

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 7 years ago by leif

Replying to martinx:

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>
license                lim                    linear_program         list                   list_plot_semilogx
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 7 years ago by leif

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

s/eulerian/Eulerian/

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

comment:29 Changed 7 years ago by leif

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: Changed 7 years ago by leif

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: Changed 7 years ago by leif

Replying to leif:

... 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 7 years ago by martinx

Replying to benjaminfjones:

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 7 years ago by martinx

Replying to leif:

Replying to leif:

... 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 7 years ago by martinx

Revised Li including reviewer patch

comment:34 Changed 7 years ago by martinx

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 7 years ago by benjaminfjones

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: Changed 7 years ago by benjaminfjones

  • Authors changed from Michael Yurko to 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 7 years ago by martinx

Replying to benjaminfjones:

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 7 years ago by benjaminfjones

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.

comment:39 Changed 7 years ago by benjaminfjones

  • Component changed from calculus to symbolics
  • Description modified (diff)
  • Keywords Li log integral symbolics calculus added
  • Status changed from needs_review to positive_review

Ok, I'm giving the most recent patch a positive review. If someone can quickly review the small, most recent fix at #11143, perhaps both of these tickets can be closed in sage-5.3.

comment:40 Changed 7 years ago by jdemeyer

  • Authors changed from martinx to Martin Cross

comment:41 Changed 7 years ago by jdemeyer

  • Merged in set to sage-5.3.rc0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.