define symbolic functions for exponential integrals

[kcrisman]

We're missing some conversions from Maxima. Like exponential integrals of various kinds.

sage: f(x) = e^(-x) * log(x+1)
sage: uu = integral(f,x,0,oo)
sage: uu
x |--> e*expintegral_e(1, 1)

See this ask.sagemath post for some details.

Current symbol conversion table

From sage.symbolic.pynac.symbol_table['maxima'] as of Sage-4.7

Maxima ---> Sage

%gamma ---> euler_gamma
%pi ---> pi
(1+sqrt(5))/2 ---> golden_ratio
acos ---> arccos
acosh ---> arccosh
acot ---> arccot
acoth ---> arccoth
acsc ---> arccsc
acsch ---> arccsch
asec ---> arcsec
asech ---> arcsech
asin ---> arcsin
asinh ---> arcsinh
atan ---> arctan
atan2 ---> arctan2
atanh ---> arctanh
binomial ---> binomial
brun ---> brun
catalan ---> catalan
ceiling ---> ceil
cos ---> cos
delta ---> dirac_delta
elliptic_e ---> elliptic_e
elliptic_ec ---> elliptic_ec
elliptic_eu ---> elliptic_eu
elliptic_f ---> elliptic_f
elliptic_kc ---> elliptic_kc
elliptic_pi ---> elliptic_pi
exp ---> exp
expintegral_e ---> En
factorial ---> factorial
gamma_incomplete ---> gamma
glaisher ---> glaisher
imagpart ---> imag_part
inf ---> +Infinity
infinity ---> Infinity
khinchin ---> khinchin
kron_delta ---> kronecker_delta
li[2] ---> dilog
log ---> log
log(2) ---> log2
mertens ---> mertens
minf ---> -Infinity
psi[0] ---> psi
realpart ---> real_part
signum ---> sgn
sin ---> sin
twinprime ---> twinprime

Summary of missing conversions

Special functions defined in Maxima

(http://maxima.sourceforge.net/docs/manual/en/maxima_16.html#SEC56)

bessel_j (index, expr)         Bessel function, 1st kind
bessel_y (index, expr)         Bessel function, 2nd kind
bessel_i (index, expr)         Modified Bessel function, 1st kind
bessel_k (index, expr)         Modified Bessel function, 2nd kind

• Notes: bessel_I, bessel_J, etc. are functions in Sage for numerical evaluation. There is also the Bessel class, but no conversions from Maxima's bessel_i etc. to Sage.

hankel_1 (v,z)                 Hankel function of the 1st kind
hankel_2 (v,z)                 Hankel function of the 2nd kind
struve_h (v,z)                 Struve H function
struve_l (v,z)                 Struve L function

• Notes: None of these functions are currently exposed at the top level in Sage. Evaluation is possible using mpmath.

assoc_legendre_p[v,u] (z)      Legendre function of degree v and order u 
assoc_legendre_q[v,u] (z)      Legendre function, 2nd kind

• Notes: In Sage we have legendre_P(n, x) and legendre_Q(n, x) both described as Legendre functions. It's not clear to me how there are related to Maxima's versions since the number of arguments differs.

%f[p,q] ([], [], expr)         Generalized Hypergeometric function
hypergeometric(l1, l2, z)      Hypergeometric function
slommel
%m[u,k] (z)                    Whittaker function, 1st kind
%w[u,k] (z)                    Whittaker function, 2nd kind

• Notes: hypergeometric(l1, l2, z) needs a conversion to Sage's hypergeometric_U. The others can be evaluated using mpmath. slommel is presumably mpmath's lommels1() or lommels2() (or both?). This isn't well documented in Maxima.

expintegral_e (v,z)            Exponential integral E
expintegral_e1 (z)             Exponential integral E1
expintegral_ei (z)             Exponential integral Ei
expintegral_li (z)             Logarithmic integral Li
expintegral_si (z)             Exponential integral Si
expintegral_ci (z)             Exponential integral Ci
expintegral_shi (z)            Exponential integral Shi
expintegral_chi (z)            Exponential integral Chi
erfc (z)                       Complement of the erf function

• Notes: The exponential integral functions expintegral_e1 and expintegral_ei (z) are called exponential_integral_1 and Ei resp. in Sage. They both need conversions. The rest need BuiltinFunction classes defined for them with evaluation handled by mpmath and the symbol table conversion added. Also, erfc is called error_fcn, so also needs a conversion.

kelliptic (z)                  Complete elliptic integral of the first 
                               kind (K)
parabolic_cylinder_d (v,z)     Parabolic cylinder D function

• Notes: kelliptic(z) needs a conversion to elliptic_kc in Sage and parabolic_cylinder_d (v,z) does not seem to be exposed at top level. It can be evaluated by mpmath.

---

Apply to the Sage library:

1. attachment: trac_11143-v2.5-rebased.4.patch
2. attachment: trac-11143-ref.2.patch
3. attachment: trac_11143-pynac-serials.patch

Depends on #13109

Component: symbolics

Keywords: ei Ei special function maxima sd32 sd40.5

Author: Benjamin Jones, Volker Braun

Reviewer: Burcin Erocal, Karl-Dieter Crisman, William Stein

Merged: sage-5.3.beta2

Issue created by migration from https://trac.sagemath.org/ticket/11143

[benjaminfjones]

comment:1

As far as I can tell, the general exponential_e function isn't available directly in Sage or in PARI (which is used to evaluate the exponential_integral_1 function in Sage).

Also, it's possible to get maxima to rewrite the exponential integrals in terms of gamma functions like so:

sage: maxima.eval('expintrep:gamma_incomplete')
'gamma_incomplete'
sage: maxima.integrate(exp(-x)*log(x+1), x, 0, oo)
%e*gamma_incomplete(0,1)
sage: N(e*gamma(0,1), digits=18)
0.596347362323194074

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma. Anyway, it should be possible to have the maxima interface (with the help of maxima itself) rewrite any exponential integral that Sage doesn't have in terms of gamma functions.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

[kcrisman]

comment:3

Replying to @benjaminfjones:

As far as I can tell, the general exponential_e function isn't available directly in Sage or in PARI (which is used to evaluate the exponential_integral_1 function in Sage).

That's unfortunate. However, mpmath seems to have it. So we could create a symbolic version and have the _eval_ method call mpmath, which we seem to be moving to.

Also, it's possible to get maxima to rewrite the exponential integrals in terms of gamma functions like so:

sage: maxima.eval('expintrep:gamma_incomplete')
'gamma_incomplete'
sage: maxima.integrate(exp(-x)*log(x+1), x, 0, oo)
%e*gamma_incomplete(0,1)
sage: N(e*gamma(0,1), digits=18)
0.596347362323194074

Interesting.

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma.

That should be ok; the whole point of the table is to convert into the Sage equivalent.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

No, that is an automatic thing that happens. It is possible to be designated an 'owner' of a ticket in a given component, which basically means you automatically get updates. If you want to 'own' it, please do! We really have plenty of special functions in Sage, but they are not always well exposed at the top level.

Incidentally, once you comment on a ticket, I believe the default is to copy you in on all replies. So you didn't have to cc: yourself specially :)

[burcin]

comment:4

Replying to @benjaminfjones:

But as you see, gamma_incomplete isn't defined in Sage either, but the table sage.symbolic.pynac.symbol_table['maxima'] lists the Sage equivalent gamma. Anyway, it should be possible to have the maxima interface (with the help of maxima itself) rewrite any exponential integral that Sage doesn't have in terms of gamma functions.

Incomplete gamma is defined in Sage. You can access it directly though incomplete_gamma() or gamma_inc(). The top level function gamma() behaves like incomplete gamma if you give it two arguments. IIRC, this is similar to maple.

By the way, the owner on the ticket is @burcin, does that mean they are working on it currently?

I am not working on it. The ticket status assigned is supposed to indicate that the owner is working on the problem, but we don't use that much either.

[burcin]

comment:5

I guess the point of this ticket is to define symbolic function in Sage to represent exponential integrals, etc. The symbolic function class handles adding stuff to the conversion table automatically.

Can we replace this ticket with several beginner tickets? One for each function we are missing.

[kcrisman]

comment:6

Replying to @burcin:

I guess the point of this ticket is to define symbolic function in Sage to represent exponential integrals, etc. The symbolic function class handles adding stuff to the conversion table automatically.

Can we replace this ticket with several beginner tickets? One for each function we are missing.

Well, that would be nice. But we could also presumably do it directly, if that would solve this problem for now. Well, either is fine as long as it were to happen. If you do split this, be sure to give a good template (I mean a link to the template, not write it yourself).

[benjaminfjones]

comment:7

I'm attempting to write a template for the expintegral_e function (denoted E_n(z) in A&S). As I'm looking through the code, I see several models used for the functions and classes in sage/functions/special.py and sage/functions/transcendental.py

• Functions like Function_exp_integral (also called Ei) are defined as classes that inherit from BuiltInFunction and call the mpmath implementation when evaluated. The function DickmanRho also does this and includes other nice methods for approximating values and power series.
• Functions like EllipticE inherit from MaximaFunction which handles evaluation, etc. through Maxima. It seems there is an advantage to the mpmath implementation because presumably the interface is faster and the precision is arbitrary (whereas Maxima is limited to 53 bits).
• Functions like Li and error_fcn are simply wrapper functions that try to evaluate the input symbolically or numerically depending on context.
• In sage/functions/trig.py there is a mixture of classes that derive from GinacFunction (and include information in their __init__ methods about conversions to other systems like Maxima or Mathematica) and also functions that derive from BuiltinFunction. It's not clear to me why some functions are Ginac and some are Builtin.

Questions:

• For the purposes of this ticket, what is recommended? @kcrisman 's comment leads me to believe that inheriting from the BuiltinFunction class and using mpmath for evaluation is preferable.
• Where should the various exponential integral special functions that we are missing go? sage/functions/special.py, sage/functions/transcendental.py, or somewhere else?

[kcrisman]

comment:8

Those are good questions.

I think the most important thing is to make sure that whatever is implemented has

• good numerical evaluation (perhaps via mpmath)
• translates back and forth to Maxima properly (for integration and limits)

My sense is that BuiltinFunction would be best. GinacFunction probably is only good for things that in fact evaluate in Ginac. This explains trig.py. This Ginac page shows that the ones which are GinacFunctions are exactly the ones which Ginac in fact has. I don't think it has most of these other functions.

As for MaximaFunction, it seems to inherit from BuiltinFunction and lives in the special functions file. This dates from the days when Sage had very few options for symbolic stuff and evaluation, and so it just does a few extra things. If we moved to mpmath for a given function, we would probably use BuiltinFunction and then add evaluation options for Maxima and add to the conversion dictionary as needed.

In fact, it wouldn't be a bad idea to have two different eval procedures if possible...

---

As for where such things go, probably it would make sense to separate a lot of these special functions out into a separate file. The distinction between transcendental and special is not totally obvious, for instance :)

---

By the way, we certainly want fewer of the wrapper functions! But that will take a lot of tedious work.

[kcrisman]

comment:9

By the way, adding this will really be a great help. Sage has so many of the things almost anybody needs, but if you have to use mpmath or GSL or R or something else directly, it sort of makes Sage a moot point. The idea is having a one-stop shop.

[benjaminfjones]

comment:10

I've uploaded a first shot at a template for the functions we want to add in this ticket. The patch exposes the general complex exponential integral function En also called Function_exp_integral_En by adding a class to sage/functions/special.py (I can change where it goes later on if needed / desired). Numerical evaluation is handled by mpmath and symbolics are handled by Sage (e.g. the derivative) and Maxima (e.g. the antiderivative).

One of the docstring examples shows that the integral of e^(-x) * log(x+1) from the ticket description is now evaluated properly.

Any comments or suggestions?

[benjaminfjones]

Description changed:

--- 
+++ 
@@ -8,7 +8,120 @@
 ```
 See [this ask.sagemath post](http://ask.sagemath.org/question/488/calculating-integral) for some details.
 
+## Current symbol conversion table
+From `sage.symbolic.pynac.symbol_table['maxima']` as of Sage-4.7
+
 ```
-sage: sage.symbolic.pynac.symbol_table['maxima']
-{'elliptic_e': elliptic_e, 'imagpart': imag_part, 'acsch': arccsch, 'glaisher': glaisher, 'asinh': arcsinh, 'minf': -Infinity, 'elliptic_f': elliptic_f, '(1+sqrt(5))/2': golden_ratio, 'inf': +Infinity, 'log(2)': log2, 'kron_delta': kronecker_delta, 'asin': arcsin, 'log': log, 'atanh': arctanh, 'brun': brun, '%pi': pi, 'acosh': arccosh, 'sin': sin, 'mertens': mertens, 'ceiling': ceil, 'infinity': Infinity, 'elliptic_ec': elliptic_ec, 'atan': arctan, 'factorial': factorial, 'twinprime': twinprime, 'khinchin': khinchin, 'catalan': catalan, 'signum': sgn, 'binomial': binomial, 'delta': dirac_delta, 'asec': arcsec, 'elliptic_kc': elliptic_kc, '%gamma': euler_gamma, 'realpart': real_part, 'elliptic_eu': elliptic_eu, 'cos': cos, 'acoth': arccoth, 'gamma_incomplete': gamma, 'li[2]': dilog, 'atan2': arctan2, 'exp': exp, 'psi[0]': psi, 'asech': arcsech, 'acos': arccos, 'acot': arccot, 'acsc': arccsc, 'elliptic_pi': elliptic_pi}
+Maxima ---> Sage
+
+%gamma ---> euler_gamma
+%pi ---> pi
+(1+sqrt(5))/2 ---> golden_ratio
+acos ---> arccos
+acosh ---> arccosh
+acot ---> arccot
+acoth ---> arccoth
+acsc ---> arccsc
+acsch ---> arccsch
+asec ---> arcsec
+asech ---> arcsech
+asin ---> arcsin
+asinh ---> arcsinh
+atan ---> arctan
+atan2 ---> arctan2
+atanh ---> arctanh
+binomial ---> binomial
+brun ---> brun
+catalan ---> catalan
+ceiling ---> ceil
+cos ---> cos
+delta ---> dirac_delta
+elliptic_e ---> elliptic_e
+elliptic_ec ---> elliptic_ec
+elliptic_eu ---> elliptic_eu
+elliptic_f ---> elliptic_f
+elliptic_kc ---> elliptic_kc
+elliptic_pi ---> elliptic_pi
+exp ---> exp
+expintegral_e ---> En
+factorial ---> factorial
+gamma_incomplete ---> gamma
+glaisher ---> glaisher
+imagpart ---> imag_part
+inf ---> +Infinity
+infinity ---> Infinity
+khinchin ---> khinchin
+kron_delta ---> kronecker_delta
+li[2] ---> dilog
+log ---> log
+log(2) ---> log2
+mertens ---> mertens
+minf ---> -Infinity
+psi[0] ---> psi
+realpart ---> real_part
+signum ---> sgn
+sin ---> sin
+twinprime ---> twinprime
 ```
+
+# Summary of missing conversions
+
+## Special functions defined in Maxima
+(http://maxima.sourceforge.net/docs/manual/en/maxima_16.html#SEC56)
+
+```
+bessel_j (index, expr)         Bessel function, 1st kind
+bessel_y (index, expr)         Bessel function, 2nd kind
+bessel_i (index, expr)         Modified Bessel function, 1st kind
+bessel_k (index, expr)         Modified Bessel function, 2nd kind
+```
+
+* Notes: bessel_I, bessel_J, etc. are functions in Sage for numerical evaluation. There is also the `Bessel` class, but no conversions from Maxima's bessel_i etc. to Sage.
+
+```
+hankel_1 (v,z)                 Hankel function of the 1st kind
+hankel_2 (v,z)                 Hankel function of the 2nd kind
+struve_h (v,z)                 Struve H function
+struve_l (v,z)                 Struve L function
+```
+
+* Notes: None of these functions are currently exposed at the top level in Sage. Evaluation is possible using mpmath.
+
+```
+assoc_legendre_p[v,u] (z)      Legendre function of degree v and order u 
+assoc_legendre_q[v,u] (z)      Legendre function, 2nd kind
+```
+
+* Notes: In Sage we have `legendre_P(n, x)` and `legendre_Q(n, x)` both described as Legendre functions. It's not clear to me how there are related to Maxima's versions since the number of arguments differs.
+
+```
+%f[p,q] ([], [], expr)         Generalized Hypergeometric function
+hypergeometric(l1, l2, z)      Hypergeometric function
+slommel
+%m[u,k] (z)                    Whittaker function, 1st kind
+%w[u,k] (z)                    Whittaker function, 2nd kind
+```
+
+* Notes: `hypergeometric(l1, l2, z)` needs a conversion to Sage's `hypergeometric_U`. The others can be evaluated using mpmath. `slommel` is presumably mpmath's `lommels1()` or `lommels2()` (or both?). This isn't well documented in Maxima.
+
+```
+expintegral_e (v,z)            Exponential integral E
+expintegral_e1 (z)             Exponential integral E1
+expintegral_ei (z)             Exponential integral Ei
+expintegral_li (z)             Logarithmic integral Li
+expintegral_si (z)             Exponential integral Si
+expintegral_ci (z)             Exponential integral Ci
+expintegral_shi (z)            Exponential integral Shi
+expintegral_chi (z)            Exponential integral Chi
+erfc (z)                       Complement of the erf function
+```
+
+* Notes: The exponential integral functions `expintegral_e1` and `expintegral_ei (z)` are called `exponential_integral_1` and `Ei` resp. in Sage. They both need conversions. The rest need `BuiltinFunction` classes defined for them with evaluation handled by mpmath and the symbol table conversion added. Also, `erfc` is called `error_fcn`, so also needs a conversion.
+
+```
+kelliptic (z)                  Complete elliptic integral of the first 
+                               kind (K)
+parabolic_cylinder_d (v,z)     Parabolic cylinder D function
+```
+
+* Notes: `kelliptic(z)` needs a conversion to `elliptic_kc` in Sage and `parabolic_cylinder_d (v,z)` does not seem to be exposed at top level. It can be evaluated by mpmath.

[benjaminfjones]

Changed keywords from ei Ei to ei Ei special function maxima