`zetazero`function which computes the

*n*th 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)

10^{0}(0.5 + 14.13472514173469379j) 0.12

10^{0}+1 (0.5 + 21.022039638771554993j) 0.05

10^{1}(0.5 + 49.773832477672302182j) 0.05

10^{1}+1 (0.5 + 52.970321477714460644j) 0.04

10^{2}(0.5 + 236.5242296658162058j) 0.32

10^{2}+1 (0.5 + 237.769820480925204j) 0.30

10^{3}(0.5 + 1419.4224809459956865j) 0.93

10^{3}+1 (0.5 + 1420.416526323751136j) 0.84

10^{4}(0.5 + 9877.7826540055011428j) 3.75

10^{4}+1 (0.5 + 9878.6547723856922882j) 3.71

10^{5}(0.5 + 74920.827498994186794j) 2.53

10^{5}+1 (0.5 + 74921.929793958414308j) 2.58

10^{6}(0.5 + 600269.67701244495552j) 2.43

10^{6}+1 (0.5 + 600270.30109071169866j) 2.89

10^{7}(0.5 + 4992381.014003178666j) 1.63

10^{7}+1 (0.5 + 4992381.2627112065366j) 2.30

10^{8}(0.5 + 42653549.760951553903j) 2.36

10^{8}+1 (0.5 + 42653550.046758478876j) 2.93

10^{9}(0.5 + 371870203.83702805273j) 6.98

10^{9}+1 (0.5 + 371870204.36631304458j) 6.72

10^{10}(0.5 + 3293531632.3971367042j) 11.48

10^{10}+1 (0.5 + 3293531632.6869557853j) 15.92

10^{11}(0.5 + 29538618431.613072811j) 53.67

10^{11}+1 (0.5 + 29538618432.07777426j) 43.00

10^{12}(0.5 + 267653395648.62594824j) 162.01

10^{12}+1 (0.5 + 267653395648.84752313j) 174.67

10^{13}(0.5 + 2445999556030.2468814j) 1262.55

10^{13}+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

^{15}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

^{13}:

>>> 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.

## 2 comments:

UBER article ! ... keep up the good job !

The new "Digital Library of Mathematical functions" is compiling an index of all software for the compuation of special functions. You should send them an email to get them to include mpmath.

Post a Comment