Assorted special functions update
June 15, 2010
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
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.
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.
>>> 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
>>> 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:= Im[LambertW[0,-1/E-10^(-900)]] Out= 0