## Friday, March 12, 2010

### Computing large zeta zeros with mpmath

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 clockmp.dps = 20for 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.12100+1  (0.5 + 21.022039638771554993j) 0.05101    (0.5 + 49.773832477672302182j) 0.05101+1  (0.5 + 52.970321477714460644j) 0.04102    (0.5 + 236.5242296658162058j)  0.32102+1  (0.5 + 237.769820480925204j)   0.30103    (0.5 + 1419.4224809459956865j) 0.93103+1  (0.5 + 1420.416526323751136j)  0.84104    (0.5 + 9877.7826540055011428j) 3.75104+1  (0.5 + 9878.6547723856922882j) 3.71105    (0.5 + 74920.827498994186794j) 2.53105+1  (0.5 + 74921.929793958414308j) 2.58106    (0.5 + 600269.67701244495552j) 2.43106+1  (0.5 + 600270.30109071169866j) 2.89107    (0.5 + 4992381.014003178666j)  1.63107+1  (0.5 + 4992381.2627112065366j) 2.30108    (0.5 + 42653549.760951553903j) 2.36108+1  (0.5 + 42653550.046758478876j) 2.93109    (0.5 + 371870203.83702805273j) 6.98109+1  (0.5 + 371870204.36631304458j) 6.721010   (0.5 + 3293531632.3971367042j) 11.481010+1 (0.5 + 3293531632.6869557853j) 15.921011   (0.5 + 29538618431.613072811j) 53.671011+1 (0.5 + 29538618432.07777426j)  43.001012   (0.5 + 267653395648.62594824j) 162.011012+1 (0.5 + 267653395648.84752313j) 174.671013   (0.5 + 2445999556030.2468814j) 1262.551013+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 1015th 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 1013:

`>>> 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.13472514173469379045725198356247027078425711569924317568556746014996342980925676494901039317156101277920297154879743676614269146988225458250536323944713778041338123720597054962195586586020055556672583601077370020541098266150754278051744259130625448197865107230493872562973832157742039521572567480933214003499046803434626731442092037738548714137831735639699536542811307968053149168852906782082298049264338666734623320078758761792005604868054356801444424651065597568665903228686510544859444320624072727032094274522213048748720924123851418351460542790152447833835425453344004487936806761697300819000731393854983736215013045167269683892003917628512321285422052396913342583227533516406016976352756375896953767492033612720925999173042707568308795118445348918008630082648312516911271068291052375961797743181517071354531677549515382893784903647470972701994848553220925357435790922612524773659551801697523346121397731600535412592674745572587780147260983080897860071253208750939599796666067537838121489190886j)`

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:

Arkapravo said...

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

Anonymous said...

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.