Assorted special functions update

Over the week-and-a-half since the last blog update, I've gotten a bunch more work done on special functions in mpmath.

Function plots in the documentation

I have started adding graphics of special functions to the documentation. So far, I've done most of the Bessel-type functions and orthogonal polynomials. Many more to come!

Example screenshot (see the Bessel functions page):

New inhomogeneous Bessel functions

There exists a large number of lesser-known special functions which are essentially variations of Bessel functions. These include functions which solve the generalized (inhomogeneous) Bessel differential equation



with some specific right-hand side g(z). New additions to mpmath in this category are the Anger function (angerj()), Weber function (webere()), and Lommel functions (lommels1(), lommels2()). See commits here and here.

More information about Anger-Weber functions and Lommel functions can be found in the DLMF.

In the near future, I will probably further improve the implementations of the main Bessel functions. The Bessel functions are mostly implemented as generic hypergeometric functions, but the standard cases can be tuned a great deal with special-purpose code.

Airy functions and related functions

The Airy functions Ai and Bi have been present for quite some time in mpmath. In a recent commit, I have rewritten them for improved rigor and better performance at high precision. There are also some new features, such as the ability to evaluate derivatives or iterated integrals of arbitrary order.

Derivatives:

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> airyai(1.5, derivative=5)
0.211387453153454489799743
>>> diff(airyai, 1.5, 5)
0.211387453153454489799743
>>> airyai(1.5, derivative=100)
-6.480220187791312407132043e+49
>>> airybi(0, derivative=1000)
3.754976097101270163249629e+854
>>> airybi(0, derivative=1001)
0.0
>>> airybi(0, derivative=1002)
3.756228172694934424506624e+856

Integrals:

>>> airyai(5, derivative=-1)
0.3332875903059178794866562
>>> quad(airyai, [0,5])
0.3332875903059178794866562
>>> airyai(-100000, derivative=-1)
-0.6665753658794626398413214

Also, functions for computing the zeros of Ai and Bi (and the first derivatives) have been added:

>>> airyaizero(1)
-2.338107410459767038489197
>>> airyaizero(2)
-4.087949444130970616636989
>>> airybizero(1)
-1.17371322270912792491998
>>> airybizero(1, derivative=1)
-2.294439682614123246622459
>>> airybizero(1, derivative=1, complex=True)
(0.2149470745374305676088329 + 1.100600143302797880647194j)
>>> airybizero(10000)
-1304.584974702601410702964
>>> airybizero(10000, complex=True)
(652.3059222438076432024695 + 1129.846189716375208308414j)

I have also implemented two new functions related to Airy functions: the Scorer functions Gi and Hi. These are available as scorergi() and scorerhi() respectively.

Here are two plots of the Gi-function, which can also be seen in the documentation:

Interval gamma functions

The interval arithmetic context now implements gamma, rgamma (reciprocal gamma function), factorial as well as loggamma for real as well as complex arguments (commit). For example:

>>> iv.dps = 10
>>> iv.gamma('50.3')
[1.96282982095908e+63, 1.96282982457481e+63]
>>> iv.gamma(iv.mpc('2.7','5.9'))
([0.00269836072064322, 0.00269836072271801] + 
 [0.0120124287790304, 0.0120124287810768]*j)

As a "practical" example, consider evaluating the Riemann-Siegel theta function which involves computing the difference of two log-gamma functions. For input with a large real part, the imaginary part in the result suffers from massive cancellation and may end up with the wrong sign:

>>> mp.dps = 15
>>> mp.siegeltheta(10**50 + 0.25j)
(5.61456887916465e+51 - 0.143091235731175j)
>>> mp.dps = 10; nprint(mp.siegeltheta(10**50 + 0.25j).imag)
-0.143091
>>> mp.dps = 100; nprint(mp.siegeltheta(10**50 + 0.25j).imag)
14.1614

With interval arithmetic, the sign uncertainty is reflected in the output:

>>> iv.dps = 15
>>> iv.siegeltheta(10**50 + 0.25j)
([5.6145688791646467648e+51, 5.6145688791646474294e+51] +
[-5.0706024009129187319e+30, 5.070602400912917606e+30]*j)
>>> iv.dps = 50
>>> iv.siegeltheta(10**50 + 0.25j)
([5614568879164647368060513633451316140100495086670736.0, 5614568879164647368060
513633451316140100495086670744.0] + 
[14.1613521236438249782320715810808676610440
8814838558203, 14.16147419395632497823207158108086766104408814838560341]*j)

As another example, consider evaluating the gamma function of a huge argument. The digits in the answer may be "wrong" because the input is converted from decimal to binary, and the gamma function is sensitive to the input being perturbed:

>>> mp.dps = 15
>>> mp.gamma('123456789012345.1')
6.11544992055093e+1686076589184486
>>> mp.dps = 30
>>> mp.gamma('123456789012345.1')
7.49032018540342193592769680745e+1686076589184486
>>> mp.dps = 60
>>> mp.gamma('123456789012345.1')
7.49032018540342058679709881225047421518964527875047787194339e+1686076589184486

With interval arithmetic, the uncertainty in the input is propagated correctly:

>>> iv.dps = 15
>>> iv.nprint(iv.gamma('123456789012345.1'), mode='diff')
[6.11545e+1686076589184486, 1.01533e+1686076589184487]
>>> iv.dps = 30
>>> iv.nprint(iv.gamma('123456789012345.1'), 20, mode='diff')
7.4903201854034[185631, 219359]e+1686076589184486
>>> iv.dps = 60
>>> iv.nprint(iv.gamma('123456789012345.1'), 50, mode='diff')
7.49032018540342058679709881225047421518964527[60898, 87505]e+1686076589184486

Rewritten Lambert W function

Lastly, the Lambert W function has received a much-needed rewrite (commit) mainly to improve evaluation very close to the branch cut along the negative axis and particularly near the branch point at -1/e for the k = -1, 0, 1 branches.

With the previous implementation, results were frequently inaccurate or ended up on the wrong branch in this region. Here are some hard cases that now work perfectly:

>>> mp.dps = 1000
>>> x = -1/e + mpf('1e-900')
>>> y = -1/e - mpf('1e-900')
>>> z = -1/e + mpf('1e-900')*1j
>>> w = -1/e - mpf('1e-900')*1j
>>> mp.dps = 25
>>> lambertw(x,0); lambertw(y,0); lambertw(z,0); lambertw(w,0)
-1.0
(-1.0 + 2.331643981597124203363536e-450j)
(-1.0 + 1.648721270700128146848651e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
>>> lambertw(x,1); lambertw(y,1); lambertw(z,1); lambertw(w,1)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-3.088843015613043855957087 + 7.461489285654254556906117j)
(-1.0 + 1.648721270700128146848651e-450j)
>>> lambertw(x,-1); lambertw(y,-1); lambertw(z,-1); lambertw(w,-1)
-1.0
(-1.0 - 2.331643981597124203363536e-450j)
(-1.0 - 1.648721270700128146848651e-450j)
(-3.088843015613043855957087 - 7.461489285654254556906117j)

To finish this post, I present the following Mathematica bug:

remote1:frejohl:[~]$ math
Mathematica 7.0 for Linux x86 (32-bit)
Copyright 1988-2008 Wolfram Research, Inc.

In[1]:= Im[LambertW[0,-1/E-10^(-900)]]

Out[1]= 0