The generalized hypergeometric function can itself be generalized endlessly. I have recently implemented three such extensions in mpmath: bilateral series, two-dimensional hypergeometric series, and q-analog (or "basic") hypergeometric series.

Bilateral series

The bilateral hypergeometric series is the simplest extension, and consists of taking the usual hypergeometric series and extending the range of summation from [0,∞) to (-∞,∞):



This series only converges when |z| = 1 and A = B, but interpreted as a sum of two ordinary hypergeometric series, it can be assigned a value for arbitrary z through the analytic continuation (or Borel regularization) of the ordinary hypergeometric series, which mpmath implements. Anyway, the convergent case is the interesting one. Here one obtains, for instance, Dougall's identity for 2H2:

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> a,b,c,d = 0.5, 1.5, 2.25, 3.25
>>> bihyper([a,b],[c,d],1)
-14.49118026212345786148847
>>> gammaprod([c,d,1-a,1-b,c+d-a-b-1],[c-a,d-a,c-b,d-b])
-14.49118026212345786148847

As an example of regularization, the divergent 1H0 series can be expressed as the sum of one 2F0 function and one 1F1 function:

>>> a = mpf(0.25)
>>> z = mpf(0.75)
>>> bihyper([a], [], z)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)
>>> hyper([a,1],[],z) + (hyper([1],[1-a],-1/z)-1)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)
>>> hyper([a,1],[],z) + hyper([1],[2-a],-1/z)/z/(a-1)
(0.2454393389657273841385582 + 0.2454393389657273841385582j)

Two-dimensional series

The most common hypergeometric series of two variables, i.e. a twodimensional series whose summand is a hypergeometric expression with respect to both indices separately, is the Appell F1 function, previously available in mpmath as appellf1:



However, much more general functions are possible. There are three other Appell functions: F2, F3, F4. The Horn functions are 34 distinct functions of order two, containing the Appell functions as special cases. The Kampé de Fériet function provides a generalization of Appell functions to arbitrary orders.

The new hyper2d function in mpmath can evaluate all these named functions, and more general functions still. The trick for speed is to write the series as a nested series, where the inner series is a generalized hypergeometric series that can be evaluated efficiently with hyper, and where the outer series has a rational recurrence formula. This rewriting also permits evaluating the analytic continuation with respect to the inner variable (as implemented by hyper).

The user specifies the format of the series in quasi-symbolic form, and the rewriting to nested form is done automatically by mpmath. For example, the Appell F1 function can be computed as

hyper2d({'m+n':[a], 'm':[b1], 'n':[b2]}, {'m+n':[c]}, x, y)

and indeed, this is essentially what appellf1 now does internally. The Appell F2-F4 functions have also been added explicitly as appellf2, appellf3, appellf4.

Hypergeometric functions of two (or more) variables have numerous applications, such as solving high-order algebraic equations, expressing various derivatives and integrals in closed form, and solving differential equations, but I have not yet found any simple examples that make good demonstrations except for F1. (I have mostly found examples of that take half a page to write down.) Any such examples for the documentation would be a welcome contribution!

Some trivial examples from the documentation are:

>>> x, y = mpf(0.25), mpf(0.5)
>>> hyper2d({'m':1,'n':1}, {}, x,y)
2.666666666666666666666667
>>> 1/(1-x)/(1-y)
2.666666666666666666666667
>>> hyper2d({'m':[1,2],'n':[3,4]}, {'m':[5],'n':[6]}, x,y)
4.164358531238938319669856
>>> hyp2f1(1,2,5,x)*hyp2f1(3,4,6,y)
4.164358531238938319669856
>>> hyper2d({'m':1,'n':1},{'m+n':1},x,y)
2.013417124712514809623881
>>> (exp(x)*x-exp(y)*y)/(x-y)
2.013417124712514809623881

An example of a Horn function, H3:

>>> x, y = 0.0625, 0.125
>>> a,b,c = 0.5,0.75,0.625
>>> hyper2d({'2m+n':a,'n':b},{'m+n':c},x,y)
1.190003093972956004227425
>>> nsum(lambda m,n: rf(a,2*m+n)*rf(b,n)/rf(c,m+n)*\
...     x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1.190003093972956004227425

This also demonstrates the recently added generic support for multidimensional infinite series in mpmath. But of course, nsum is much slower than hyper2d.

Hypergeometric q-series

Before introducing the q-analog of the hypergeometric series, I should introduce the q-Pochhammer symbol,


This itself is a new function in mpmath, implemented as qp(a,q,n) (with two- and one-argument forms qp(a,q) and qp(q) also permitted) and is the basis for more general computation involving q-analogs. The q-factorial and q-gamma function have also been added (as qfac and qgamma), but are not yet documented.

The q-analogs have important applications in number theory. As a very neat example, numerically computing the Taylor series of 1/(q, q)_∞ with mpmath gives

>>> taylor(lambda q: 1/qp(q), 0, 10)
[1.0, 1.0, 2.0, 3.0, 5.0, 7.0, 11.0, 15.0, 22.0, 30.0, 42.0]

These are the values of the partition function P(n) for n = 0,1,2,..., i.e. the number of ways of writing n as a sum of positive integers.

Replacing the rising factorials (Pochhammer symbols) in the generalized hypergeometric series with their q-analogs gives the hypergeometric q-series or basic hypergeometric series



This function is implemented as qhyper. Like hyper, it supports arbitrary combinations of real and complex arguments (assuming |q| 1). Some examples from the documentation:

>>> qhyper([0.5], [2.25], 0.25, 4)
-0.1975849091263356009534385
>>> qhyper([0.5], [2.25], 0.25-0.25j, 4)
(2.806330244925716649839237 + 3.568997623337943121769938j)
>>> qhyper([1+j], [2,3+0.5j], 0.25, 3+4j)
(9.112885171773400017270226 - 1.272756997166375050700388j)

Like hyper, it automatically ensures accurate evaluation for alternating series:

>>> q = 0.998046875
>>> mp.dps=5; qhyper([2],[0.5], q, -0.5)
6.6738e-69
>>> mp.dps=15; qhyper([2],[0.5], q, -0.5)
6.67376764851253e-69
>>> mp.dps=25; qhyper([2],[0.5], q, -0.5)
6.673767648512527695718826e-69
>>> mp.dps=100; qhyper([2],[0.5], q, -0.5)
6.673767648512527695718826106778799352769798151218768443717704076836963752188876
561888441662933081804e-69

With the q-analog of the generalized hypergeometric function implemented, it becomes possible to compute q-exponentials, q-sines, q-orthogonal polynomials, q-Bessel functions, and pretty much anything else. If there is interest, such function could be added explicitly to mpmath.

For more information and examples of the functions discussed in this post, see the sections on hypergeometric functions and q-functions (a little terse at the moment) in the mpmath documentation.

This post demonstrates how the knapsack problem can be used to construct a digital signature. As with the knapsack systems themselves, it has been broken. However, it is an elegant system, and also very fast. It was designed by Shamir in 1978; our description is based on an article by Andrew Odlyzko. It’s quite [...]

Juan Arias de Reyna, who wrote the code used in mpmath for evaluating the Riemann zeta function near the critical strip, has contributed a new implementation of the zetazero function which computes the nth zero on the critical line.

The old version (written by myself) used a lookup table for initial approximations, so it was limited to computing the first few zeros. Juan's code calculates the position of arbitrary zeros using Gram's law, with a sophisticated algorithm (about 300 lines of code) to find the correct zero when Gram's law (which is actually just a heuristic) fails.

A bunch of zeros, from small to large:

from mpmath import *
from timeit import default_timer as clock
mp.dps = 20
for n in range(14):
    t1 = clock()
    v1 = zetazero(10**n)
    t2 = clock()
    v2 = zetazero(10**n+1)
    t3 = clock()
    print "10<sup>%i</sup>  "%n, v1, "%.2f" % (t2-t1)
    print "10<sup>%i</sup>+1"%n, v2, "%.2f" % (t3-t2)


n      value                          time (s)
100    (0.5 + 14.13472514173469379j)  0.12
100+1  (0.5 + 21.022039638771554993j) 0.05
101    (0.5 + 49.773832477672302182j) 0.05
101+1  (0.5 + 52.970321477714460644j) 0.04
102    (0.5 + 236.5242296658162058j)  0.32
102+1  (0.5 + 237.769820480925204j)   0.30
103    (0.5 + 1419.4224809459956865j) 0.93
103+1  (0.5 + 1420.416526323751136j)  0.84
104    (0.5 + 9877.7826540055011428j) 3.75
104+1  (0.5 + 9878.6547723856922882j) 3.71
105    (0.5 + 74920.827498994186794j) 2.53
105+1  (0.5 + 74921.929793958414308j) 2.58
106    (0.5 + 600269.67701244495552j) 2.43
106+1  (0.5 + 600270.30109071169866j) 2.89
107    (0.5 + 4992381.014003178666j)  1.63
107+1  (0.5 + 4992381.2627112065366j) 2.30
108    (0.5 + 42653549.760951553903j) 2.36
108+1  (0.5 + 42653550.046758478876j) 2.93
109    (0.5 + 371870203.83702805273j) 6.98
109+1  (0.5 + 371870204.36631304458j) 6.72
1010   (0.5 + 3293531632.3971367042j) 11.48
1010+1 (0.5 + 3293531632.6869557853j) 15.92
1011   (0.5 + 29538618431.613072811j) 53.67
1011+1 (0.5 + 29538618432.07777426j)  43.00
1012   (0.5 + 267653395648.62594824j) 162.01
1012+1 (0.5 + 267653395648.84752313j) 174.67
1013   (0.5 + 2445999556030.2468814j) 1262.55
1013+1 (0.5 + 2445999556030.6222451j) 1456.38

The function correctly separates close zeros, and can optionally return information about the separation. For example:

>>> mp.dps = 20; mp.pretty = True
>>> zetazero(542964976,info=True)
((0.5 + 209039046.578535272j), [542964969, 542964978], 6, '(013111110)')

The extra information is explained in the docstring for zetazero:

[This means that the] zero is between Gram points 542964969 and 542964978, it is the 6-th zero between them. Finally (01311110) is the pattern of zeros in this interval. The numbers indicate the number of zeros in each Gram interval, (Rosser Blocks between parenthesis). In this case only one Rosser Block of length nine.

Juan reports having computed the 10¹⁵th zero, which is

208514052006405.46942460229754774510611

and verifying that it satisfies the Riemann hypothesis. And fortunately, the values of large zeros computed with mpmath seem to agree with those reported by Andrew Odlyzko and Xavier Gourdon.

The new zetazero function in mpmath is indeed able to separate the closest known pair of zeros, found in Xavier Gourdon's exhaustive search up to 10¹³:

>>> mp.dps = 25
>>> v1 = zetazero(8637740722917, info=True)
>>> v2 = zetazero(8637740722918, info=True)
>>> v1
((0.5 + 2124447368584.392964661515j), [8637740722909L, 8637740722925L], 7L,
'(1)(1)(1)(1)(1)(1)(02)(1)(02)(1)(1)(20)(1)')
>>> v2
((0.5 + 2124447368584.392981706039j), [8637740722909L, 8637740722925L], 8L,
'(1)(1)(1)(1)(1)(1)(02)(1)(02)(1)(1)(20)(1)')
>>> nprint(abs(v1[0]-v2[0]))
1.70445e-5

This computation takes over 2 hours on my laptop.

Of course, zeros can be computed to any desired precision:

>>> mp.dps = 1000
>>> print zetazero(1)
(0.5 + 14.1347251417346937904572519835624702707842571156992431756855674601499634
29809256764949010393171561012779202971548797436766142691469882254582505363239447
13778041338123720597054962195586586020055556672583601077370020541098266150754278
05174425913062544819786510723049387256297383215774203952157256748093321400349904
68034346267314420920377385487141378317356396995365428113079680531491688529067820
82298049264338666734623320078758761792005604868054356801444424651065597568665903
22868651054485944432062407272703209427452221304874872092412385141835146054279015
24478338354254533440044879368067616973008190007313938549837362150130451672696838
92003917628512321285422052396913342583227533516406016976352756375896953767492033
61272092599917304270756830879511844534891800863008264831251691127106829105237596
17977431815170713545316775495153828937849036474709727019948485532209253574357909
22612524773659551801697523346121397731600535412592674745572587780147260983080897
860071253208750939599796666067537838121489190886j)

The only real limitation of the code, as far as I can tell, is that it's impractical for computing a large number of consecutive zeros. For that one needs to use something like the Odlyzko–Schönhage algorithm for multi-evaluation of the zeta function.

A nice exercise would be to parallelize the zeta code using multiprocessing. The speedup for large zeros should be essentially linear with the number of processors.

I have a fairly large backlog of changes in mpmath to blog about. For today, I'll just briefly mention the most recent one, which is a set of optimizations to elementary functions. The commit is here. I intend to provide much faster Cython versions of the elementary functions some time in the future, but since there was room for optimizing the Python versions, I decided to invest a little time in that.

Most importantly, the asymptotic performance of trigonometric function has been improved greatly, due to a more optimized complexity reduction strategy. The faster strategy was previously used only for exp, but cos/sin are virtually identical in theory so fixing this was mostly a coding problem. In fact, the high-precision code for exp/cosh/sinh/cos/sin is largely shared now, which is nice because only one implementation needs to be optimized.

The following image shows performance for computing cos(3.7). The red graph is the old implementation; blue is new. The top subplot shows evaluations/second at low precision, i.e. higher is better, and the bottom subplot shows seconds/evaluation at higher precision, i.e. lower is better.



Already at a few thousand bits, the new code is more than twice as fast, and then it just gets better and better. There is also a little less overhead at low precision now, so the new code is uniformly faster.

Also exp, cosh and sinh have been optimized slightly. Most importantly, there is a little less overhead at low precision, but asymptotic speed has also improved a bit.

Performance for computing cosh(3.7):



Performance for computing exp(3.7):


I should note that the performance depends on some tuning parameters which naturally are system-specific. The results above are on a 32-bit system with gmpy as the backend. I hope to later implement optimized tuning parameters for other systems as well.

Elementary function performance is of course important globally, but it's especially important for some specific functions. For example, the speed of the Riemann zeta function depends almost proportionally on the total speed of evaluating one exp, one log, one sine, and one cosine, since the main part of the work consists of summing the truncated L-series



With the improved elementary functions, and some overhead removals to the zetasum code in the same commit, the Riemann zeta function is up to 2x faster now, and about 50% faster on the critical line.

If you have a lengthy Sage Notebook and you want to quickly jump to a certain header or paragraph as a reference, use HTML anchors. They work like follows:

1. An <a name="anchorname"></a> tag must be inserted at the position you want to jump to.
2. Reference it via <a href="#anchorname">some text</a>.

(The ideas in this post are based on a worksheet by Jason Grout. Thanks Jason!)

Sage 2D plots use the excellent library matplotlib as a backend, which gives us beautiful publication-quality graphics. Sage only exposes a small subset of matplotlib's capabilities, and this subset is generally sufficient in day-to-day work. However, especially if you are preparing a graph for a class or a presentation, you might need to tweak things to get it to look "just right".

I recently had to produce some graphs of trigonometric functions for my calculus class, and the default tick locations and labels are not really well-suited. Here is one example (code and graph):

var('x')
p = plot(2*sin(x+pi/3)+2, (x, -pi/3, 5*pi/3))
p.save('2sinx_pi_3_2_original.png') 

There is too much numerical information displayed, and a lot of it is just distracting for this function. I would much rather just see some multiples of (say) marked along the x-axis, and a few significant numbers marked along the y-axis. For this we need to access matplotlib directly. Here is the new code, and the resulting graph:

var('x')
p = plot(2*sin(x+pi/3)+2, (x, -pi/3, 5*pi/3))
fig = p.matplotlib(axes_labels=['$x$', '$y$'])
from matplotlib.backends.backend_agg import FigureCanvasAgg
fig.set_canvas(FigureCanvasAgg(fig))
axes = fig.get_axes()[0]
axes.minorticks_off()
axes.set_xticks([float((j-1/3)*pi) for j in range(0,3)])
axes.set_xticklabels(['$-\pi/3$', '$2\pi/3$', '$5\pi/3$'])
axes.set_yticks([float(2), float(4)])
fig.suptitle('$2\,\sin(x+\pi/3)+2$', color='navy')
fig.savefig('2sinx_pi_3_2_matplotlib.png')

One of the reasons for the popularity of modular forms in number theory is that they are highly computable objects. There is a variety of ways in which you can get your hands on explicit examples (see for instance William Stein's Modular forms: a computational approach or Lloyd Kilford's Modular forms: a classical and computational introduction). And, of course, you can (and probably should) use Sage to generate and play with such examples.

One of the reasons for the lesser popularity of various generalisations of modular forms is that they currently lack good computational infrastructure. The good news is that this is a fairly active field of research. Thinking specifically of Siegel modular forms, there are rumours of an upcoming paper (and code) by Raum-Ryan-Skoruppa-Tornaría. While we wait for this, there is a recent preprint by Martin Raum with the descriptive title Efficiently generated spaces of classical Siegel modular forms and the Böcherer conjecture. More precisely, its algorithmic focus is on genus 2, level 1, scalar-valued Siegel modular forms of high weight and with rational coefficients, and he conjectures the existence of a simple, nice-looking basis for these vector spaces (which we are encouraged to think of as an analogue of the Victor Miller basis for classical modular forms). He checked the conjecture computationally up to weight 172 (which is quite impressive) by using Sage.

This is excellent news, since the Victor Miller basis is a particularly efficient method for computing with classical modular forms, and it is wonderful to have this in the Siegel context. He shows off some of the power of this approach with a few interesting applications.

I've been amusing myself all day with doing some complex plotting in Sage. The default color scheme has 0=black, infinity=white, and red=real. A lot of the time, this doesn't work too well, but if you plot f/abs(f) the color will just depend on the complex argument. For example, if f = (z^5+1)/(z^5-1) we get:



Or a more baroque example, f = z*log(z)*exp((z^5+1)/(z^5-1)^(1/2)):




Sometimes the default coloring is helpful, for example when showing the essential singularity at 0 of exp(1/z):

The newest release of Sage is now out. See this announcement for full details, including a list of tickets closed.

There's one improvement I'm particularly happy about: Dima Pasechnik upgraded the version of GAP included in Sage to the latest upstream, namely 4.4.12. The previous version (4.4.10) was seriously outdated and had several known issues.

Thanks to Dima for the work in upgrading the spkg, and to David Joyner, William Stein, Minh Van Nguyen, and Robert Bradshaw for reviewing the various tickets involved in this upgrade.

Pre-requisites. I don’t know how much of it its really needed, but I install the whole Sun’s Java 6 packages, $ sudo apt-get install sun-java6-jdk sun-java-6-jre sun-java6-pluging sun-java6-source sun-java6-bin For sure you’ll need the jdk and plugin one… If your internet connection is slow you might try installing these two first and try the above command. Configuration [...]

The gamma function is probably the most important nonelementary special function. The gamma function or ratios of gamma functions appear in everything from the functional equation of the Riemann zeta function to the normalization factors of hypergeometric functions. Even expressions that are usually rational functions (such as binomial coefficients) are often most conveniently expressed (or in software, implemented) in terms of the gamma function since this automatically provides correct asymptotics and extension to noninteger parameters. Needless to say, a well-implemented gamma function is fundamental for a special functions library.

The gamma function implementation used in mpmath until now is one of the oldest parts of the code (I wrote it around three years ago), and although it has done its job, it has some flaws. It mainly uses Spouge's approximation, which is fairly efficient but nonoptimal for asymptotically large arguments. The implementation also has some subtle precision issues.

The new gamma function added in this commit is the result of over a year of intermittent work. I did most of the work during last summer but didn't have time to finish it until recently. The improvements are essentially the following:

• Removal of various overheads


• Special-casing half-integers


• Using Taylor series for small real arguments


• Using Stirling's series for complex and/or large arguments


• Optimizing for computing loggamma


• Handling values near poles and near the real axis extremely carefully

The difference in performance is quite dramatic. The smallest speedup occurs for small complex arguments, but it turns out that even there, the Stirling series is faster than Spouge's approximation (at least with the removal of overheads in the new implementation) despite the need to perform argument transformations. For real arguments, and large complex arguments, the new implementation is consistently faster by a large factor, especially for loggamma.

I'll show some benchmarks. The times are in milliseconds; the first value in each cell is the time for the old implementation and the second is the time for the new one, followed by the speedup. The results are on my laptop, with gmpy (I'll try with Sage later):

Firstly, results for gamma(z):

z	dps=15	dps=50	dps=250	dps=1000
3.0	0.0086 0.0077 (1.13x)	0.0089 0.0076 (1.16x)	0.0087 0.0075 (1.16x)	0.0085 0.0077 (1.12x)
5.27	0.1200 0.0342 (3.51x)	0.2513 0.0545 (4.61x)	2.5964 0.4220 (6.15x)	80.7089 11.3011 (7.14x)
2.5	0.1255 0.0158 (7.94x)	0.2398 0.0138 (17.41x)	2.5519 0.0144 (177.72x)	79.0651 0.0149 (5319.87x)
-3.7	0.2518 0.0457 (5.51x)	0.4006 0.0615 (6.52x)	3.1073 0.4283 (7.26x)	87.4200 11.3981 (7.67x)
(-3.0+4.0j)	0.6205 0.3942 (1.57x)	1.0095 0.6093 (1.66x)	9.0896 5.9099 (1.54x)	253.9219 191.3358 (1.33x)
(3.25-4.125j)	0.3548 0.2864 (1.24x)	0.7456 0.4791 (1.56x)	8.6651 5.6367 (1.54x)	252.7598 189.2187 (1.34x)
(15.0+20.0j)	0.3980 0.1973 (2.02x)	0.7526 0.4496 (1.67x)	8.6583 5.4604 (1.59x)	250.8646 188.5849 (1.33x)
52.1	0.1661 0.0586 (2.83x)	0.3050 0.0894 (3.41x)	2.9862 0.5685 (5.25x)	84.8334 12.4169 (6.83x)
123.25	0.1643 0.0638 (2.57x)	0.3071 0.0963 (3.19x)	2.8908 0.9194 (3.14x)	83.8168 14.6588 (5.72x)
-200.7	0.3011 0.1485 (2.03x)	0.4897 0.1914 (2.56x)	3.5145 1.1359 (3.09x)	89.8617 17.5500 (5.12x)
(100.25+100.0j)	0.4238 0.1857 (2.28x)	0.7641 0.2757 (2.77x)	8.6807 3.1403 (2.76x)	249.6731 171.4369 (1.46x)
12345678.9	0.3032 0.0701 (4.32x)	0.4382 0.0857 (5.11x)	3.4078 0.3257 (10.46x)	88.4082 5.2828 (16.74x)
(0.0+12345678.9j)	0.7814 0.1409 (5.54x)	1.1671 0.1802 (6.48x)	10.3985 0.6169 (16.86x)	265.6192 8.3569 (31.78x)
1.0e+20	0.8833 0.0798 (11.06x)	1.3019 0.0899 (14.48x)	6.4010 0.2939 (21.78x)	107.0924 3.5141 (30.47x)
(0.0+1.0e+20j)	3.5025 0.1515 (23.11x)	3.6163 0.1834 (19.72x)	19.5676 0.5308 (36.86x)	322.9033 6.2148 (51.96x)

Secondly, results for loggamma(z):

z	dps=15	dps=50	dps=250	dps=1000
3.0	0.0700 0.0115 (6.07x)	0.0677 0.0116 (5.84x)	0.0683 0.0119 (5.73x)	0.0683 0.0121 (5.65x)
5.27	0.2207 0.0564 (3.91x)	0.3548 0.0796 (4.45x)	2.8925 0.5221 (5.54x)	88.3556 12.6452 (6.99x)
2.5	0.2210 0.0319 (6.93x)	0.3565 0.0374 (9.53x)	2.8595 0.1011 (28.28x)	81.8333 1.3546 (60.41x)
-3.7	0.4392 0.0782 (5.62x)	0.6418 0.0987 (6.50x)	3.5767 0.5429 (6.59x)	88.8484 12.7363 (6.98x)
(-3.0+4.0j)	1.1235 0.4942 (2.27x)	1.5685 0.7224 (2.17x)	10.0879 6.0459 (1.67x)	267.2893 198.7613 (1.34x)
(3.25-4.125j)	0.8789 0.2704 (3.25x)	1.2326 0.4752 (2.59x)	9.7501 5.4572 (1.79x)	264.7649 190.4007 (1.39x)
(15.0+20.0j)	0.8983 0.1248 (7.20x)	1.2923 0.4221 (3.06x)	9.7917 5.2766 (1.86x)	265.2711 187.9775 (1.41x)
52.1	0.2604 0.0812 (3.21x)	0.4186 0.1140 (3.67x)	3.2493 0.6890 (4.72x)	87.0111 13.9280 (6.25x)
123.25	0.2629 0.0478 (5.50x)	0.4227 0.0692 (6.11x)	3.1110 1.0216 (3.05x)	85.9028 15.9917 (5.37x)
-200.7	0.5439 0.1631 (3.34x)	0.7465 0.2221 (3.36x)	3.9305 1.2537 (3.14x)	92.5380 18.8309 (4.91x)
(100.25+100.0j)	0.9103 0.1179 (7.72x)	1.3332 0.1849 (7.21x)	10.1655 3.1083 (3.27x)	266.3565 170.3099 (1.56x)
12345678.9	0.4026 0.0491 (8.21x)	0.5439 0.0525 (10.36x)	3.7183 0.2051 (18.13x)	89.9402 4.0044 (22.46x)
(0.0+12345678.9j)	1.2628 0.0765 (16.52x)	1.7085 0.0858 (19.91x)	11.5023 0.2749 (41.84x)	275.6174 4.4548 (61.87x)
1.0e+20	0.9820 0.0495 (19.84x)	1.3882 0.0540 (25.72x)	6.6983 0.1481 (45.21x)	110.1757 2.1954 (50.19x)
(0.0+1.0e+20j)	4.0385 0.0788 (51.24x)	4.3758 0.0833 (52.54x)	20.9478 0.1936 (108.20x)	335.9407 2.3911 (140.50x)

All times are from a warm cache. The new implementation has slightly longer precomputation time (for the Taylor series coefficients) at high precision, and the Taylor series therefore isn't used above a certain precision (around 1000 digits) despite being much faster on subsequent calls. I will probably add some user-visible function for controlling this.

Finally, several optimizations are possible still. The most important would be to do the implementation in Cython. (Improving the underlying elementary functions will also speed up the gamma function.) Potential algorithmic improvements include faster handling of near-integers (using Taylor series with a reduced number of terms), avoiding the reflection formula for complex arguments in the left half-plane but not too close to the negative half-axis, and performing the argument transformations with a reduced number of multiplications. But I doubt I will have time (or interest -- the present code was demanding enough to finish) to work of any of these things, at least any time soon.

The speedups naturally also affect performance of many other functions, e.g. hypergeometric functions; I might post some benchmarks of that later.

It is perhaps not so surprising that modular forms mod p do not always behave the way one would expect based on properties of modular forms in characteristic 0. Here is an example of such a discrepancy. (This is of course "well-known", and as with a lot of "well-known" things, it's not always easy to find a reference.)

It is a classical result that the space of modular forms has a basis of Hecke eigenforms. (See for instance Section 4.2.3 in Kilford's Modular forms: a classical and computational introduction.) If you combine this with the q-expansion principle, you get a multiplicity one result, i.e. the eigenspace associated to any system of Hecke eigenvalues occurring in M is one-dimensional. There are similar results for level N, where one needs to differentiate between oldforms and newforms. (See for instance Section 5.8 in Diamond-Shurman's A first course in modular forms.)

Let's now see what happens when we instead work with , the space of modular forms mod p. It turns out that this space need not have a basis of Hecke eigenforms. Here is a sketch of proof: suppose these spaces have a basis of Hecke eigenforms, then each eigensystem has multiplicity one by the q-expansion principle. But a theorem of Serre-Tate says that there are only finitely many Hecke eigensystems mod p (all weights put together). Say this number is H. Choosing a weight k such that the space of modular forms of weight k has dimension greater than H gives us more than H eigensystems in weight k, contradiction.

The argument tells us where to look for such situations. I'll restrict my attention to cusp forms, but this is not important. If p is 5, then there are only 4 Hecke eigensystems in all weights. So taking a space of dimension at least 5 would do the trick, and indeed the space of cusp forms of weight 60 works. But it is possible to find smaller such examples (of course, the space must be at least two-dimensional for this to work). For p=5, the earliest example comes up in weight 36. Sage makes it easy to explore this as follows. First, we need a basis for the space of cusp forms mod 5 of weight 36, and we'll use the corresponding Victor Miller basis:

sage: bas = victor_miller_basis(k=36, prec=20, cusp_only=True)
sage: bas_mod5 = [ f.change_ring(GF(5)) for f in bas ]

The second line reduces the original Victor Miller basis (which lives over the integers) modulo 5. Now we figure out the matrix of the Hecke operator on this basis. This should just be a matter of doing hecke_operator_on_basis(bas_mod5, 2, 36), but at the moment that doesn't work over finite fields (for silly reasons, and this should be an easy bug to fix).

No matter, we can just do it by hand:

sage: f1, f2, f3 = bas_mod5; print f1; print f2; print f3
q + 3*q^4 + 2*q^6 + 2*q^9 + 2*q^11 + O(q^12)
q^2 + 4*q^4 + 2*q^5 + q^6 + 4*q^7 + 4*q^8 + 3*q^9 + 2*q^10 + O(q^12)
q^3 + 3*q^4 + 4*q^5 + 2*q^6 + q^7 + 3*q^8 + q^9 + 4*q^10 + O(q^12)
sage: hecke_operator_on_qexp(f1, 2, 36)
q^2 + 2*q^3 + O(q^6)

And herein lies the beauty of the Victor Miller basis: we immediately see that the last result is really . Similarly:

sage: hecke_operator_on_qexp(f2, 2, 36)
q + 4*q^2 + q^3 + 2*q^4 + 2*q^5 + O(q^6)
sage: hecke_operator_on_qexp(f3, 2, 36)
3*q^2 + 2*q^3 + 3*q^4 + 4*q^5 + O(q^6)

Now we can write down the matrix of and see what happens:

sage: t2 = matrix(GF(5), [[0, 1, 0], [1, 4, 3], [2, 1, 2]]); t2
[0 1 0]
[1 4 3]
[2 1 2]
sage: t2.jordan_form()
[4|0 0]
[-+---]
[0|1 1]
[0|0 1]

And there you have it: since the matrix is not diagonalisable, this space of cusp forms cannot even have a basis of eigenvectors for , let alone a basis of eigenvectors for all the Hecke operators.

Another "small" example is p=67, k=32, where the space of cusp forms is two-dimensional and the matrix of has Jordan normal form

sage: t2 = matrix(GF(67), [[0, 1], [5, 28]])
sage: t2.jordan_form()
[14  1]
[ 0 14]

There doesn’t seem to have been much done in the way of studying fractals of iterated functions involving gamma. At least nothing that I could find.

So I started playing around with it. It’s a bit difficult to work with as sage (4.3.2) keeps crashing when I try to plot it (try to run complex_plot(gamma(gamma(x)), [-10,10], [-10,10]) or complex_plot(gamma(gamma(gamma(x))), [-5,5], [-5,5]) ), but I still got some pictures…

Another new feature in mpmath: support for multidimensional series in nsum. Some very simple examples (finite and infinite summations can be combined in any desired order):

>>> nsum(lambda i,j: 2**(-i-j), [0,inf], [0,inf])
4.0
>>> nsum(lambda i,j,k: 2**(-i-j-k), [0,inf], [0,inf], [0,inf])
8.0
>>> nsum(lambda i,j,k,l: i/2**j*(i+l)/2**k, [1,2], [0,inf], [1,inf], [2,3])
50.0
>>> nsum((lambda i,j,k,l: 1/(2**(i**2+j**2+k**2+l**2))), *([[-inf,inf]]*4))
20.5423960756379
>>> jtheta(3,0,0.5)**4
20.5423960756379

One could of course also make nested calls to nsum, but having a direct syntax is much more convenient. Nested evaluation is also usually inefficient for convergence acceleration, so this is not what nsum does internally. Instead, it combines all the infinite summations to a single summation over growing hypercubes. The distinction is very important for conditionally convergent series.

For example, nsum can now directly evaluate the Madelung constant in three dimensions (ignore=True is to ignore the singular term at the origin):

>>> nsum(lambda i,j,k: (-1)**(i+j+k)/(i**2+j**2+k**2)**0.5,
...     [-inf,inf], [-inf,inf], [-inf,inf], ignore=True)
-1.74756459463318

While this evaluation takes several seconds, it is somewhat remarkable that a precise value can be obtained at all. A better way to compute the Madelung constant is using the following rapidly convergent 2D series:

>>> f = lambda i,j: -12*pi*sech(0.5*pi*sqrt((2*i+1)**2+(2*j+1)**2))**2
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318
>>> mp.dps = 100
>>> nsum(f, [0,inf], [0,inf])
-1.74756459463318219063621203554439740348516143662474175815282535076504062353276
117989075836269460789

Another nice application is to evaluate Eisenstein series directly from the definition:

>>> tau = 1j
>>> q = qfrom(tau=tau)
>>> nsum(lambda m,n: (m+n*tau)**(-4), [-inf,inf], [-inf,inf], ignore=True)
(3.1512120021539 + 0.0j)
>>> 0.5*sum(jtheta(n,0,q)**8 for n in [2,3,4])*(2*zeta(4))
(3.1512120021539 + 0.0j)

I mess around with fractals a lot, and sometimes stumble on some interesting things. Recently I came across some odd jagged/discontinuous fractals:

This tends to happen in fractals as soon as I involved non-Integer powers.

Why is this happening? To understand this, we need to look at some basic complex analysis.

Recall that a complex number can be interpreted as a vector like . But we can also think if a vector in terms of direction and magnitude. For our purposes, we will think of vectors as a magnitude and an angle from the positive portion of the real number line (sometimes called the `argument,’ we will just call it the angle).

Complex numbers have the interesting property that when we multiply them their magnitudes multiply like normal numbers but their angles add. For example, has an angle of so when we multiply two of them we end up with a magnitude of and an angle of and thus .

When we raise a value to a real power, the angle is multiplied by that value. We can visualise this as:

It’s fairly clear that for the inverse, for any value, there are two valid answers.

We can also look at it with an alternative visualisation I just cooked up:

Notice that the angle of the ray of discontinuity is arbitrary, but its existence is inevitable.

(For more on complex analysis, I recommend Visual Complex Analysis by Tristan Needham (Google Books, Amazon). It’s awesome!)

(Sage users: Here is a generalised version of the function I used to make the above image. It takes a value n and returns a visualisation of the nth root:

def root_visualize(n):
 l=[]
 for j in range(n*8):
 l.append(complex_plot(lambda x: (x*exp(-sqrt(-1)*j*pi/4))^(1/n)*exp(sqrt(-1)*j*pi/n/4),[-5,5],[-5,5]))
 print j
 return animate(l)

Be warned: it is slow)

Now, notice that is the same as . Thus, it has the same discontinuity. That seems like a reasonable answer to “why is this so jagged?”

But it leads to some interesting questions. The intuitive reaction is (for me at least): part of the fractal is missing. You’re cutting apart the nice, smooth mutlifunction and thus are only seeing part of the fractal. Well, at iterations there possible ways to have cut this hypothetical super-fractal (as I will refer to it).

Let’s look at the Julia super-fractal of for four iterations…

If you stare at all the possibilities of the forming set you will notice some patterns and what look like reciprocals and other patterns. One may intuitively wish to add them to extract a single fractal, but this isn’t possible since they’d cancel… Perhaps adding the absolute values?

(For more on the formation of fractals, see my previous post.)

Of course, one really should look at a larger number of iterations than I am. There is the small problem of this being , but it should be easy to parallelize…

So that’s the extent that I’ve explored this to so far. I’d be thrilled to get some feedback. Do you have any thoughts about this? Know of similar things that have been studied? (I’ve never formally studied Chaos Theory and suspect I have lots of holes in my knowledge. Reminder to self: take/audit course on chaos theory next year…)

The Floyd-Warshall-Roy algorithm is an algorithm for finding shortest paths in a weighted, directed graph. It allows for negative edge weights and detects a negative weight cycle if one exists. Assuming that there are no negative weight cycles, a single execution of the FWR algorithm will find the shortest paths between all pairs of vertices. This [...]

Here are some notes to myself regarding compiling Sage on my HP Mini 2140 with an Atom N270 CPU. I'm running Linux Ubuntu 9.10. This successfully compiles Sage 4.3 and 4.3.1 and might also work for later version.
Download


Links