# Fun with zeta functions

February 4, 2009

In mpmath-trunk there is now an implementation of the Riemann-Siegel Z function as well as the related Riemann-Siegel theta function. There is also a function for computing Gram points.

A picture is worth a thousand words. The Z function is a kind of real-valued version of the Riemann zeta function on the critical strip (compare to this previous plot):

>>> from mpmath import *>>> plot(siegelz, [0,80])

The implementation is entirely generic, so it works in the complex numbers:

>>> cplot(siegelz, points=50000)

>>> cplot(siegelz, [-25,25], [-25,25], points=50000)

Gram points are useful for locating zeros of the zeta/Z functions. Huge Gram points themselves are easily located:

>>> mp.dps = 50>>> print grampoint(10**10)3293531632.7283354545611526800803306343201329980271

Unfortunately, evaluation of the Riemann zeta or Riemann-Siegel Z functions is not feasible for such large inputs with the present zeta function implementation in mpmath. Lowering expectations a bit, one can still compute some fairly large zeros:

>>> g1 = grampoint(1000)>>> g2 = grampoint(1001)>>> r = findroot(siegelz, [g1, g2], solver='illinois')>>> print r1421.8505671870486539107068075509847506037846486061>>> print g1 < r < g2True>>> nprint(siegelz(r))-1.30567e-51

(I chose the Illinois solver here because it combines bracketing with the fast convergence of the secant method.)

A plot of the initial couple of Gram points:

>>> mp.dps = 15>>> plot(lambda x: grampoint(floor(x)), [0, 30])

The following visualization of Gram points is perhaps more illustrative. It is usually the case that the Gram points and sign changes (i.e. zeros) of Z alternate with each other:

>>> gs = [grampoint(n) for n in range(25)]>>> def marker(x):...     for g in gs:...         if abs(g-x) < 0.2:...             return 1...     return 0...>>> plot([siegelz, marker], [0,50])

Put another way, it is usually the case that (-1)n Z(gn) is positive (this is Gram’s law):

>>> (-1)**10 * siegelz(grampoint(10)) > 0True>>> (-1)**11 * siegelz(grampoint(11)) > 0True>>> (-1)**12 * siegelz(grampoint(12)) > 0True

The first exception occurs at n = 126 (more exceptions are listed in OEIS A114856). The Gram point is very close to a root:

>>> (-1)**126 * siegelz(grampoint(126)) > 0False>>> mp.dps = 50>>> print grampoint(126)282.4547208234621746108397940690599354048302480008>>> print findroot(siegelz, grampoint(126))282.46511476505209623302720118650102420550683628035>>> print siegelz(grampoint(126))-0.02762949885719993669875120344077657187547768743854

Now, it takes only a few lines of code to enumerate exceptions to Gram’s law, up to say n = 2000:

>>> mp.dps = 10>>> for k in range(2000):...     g = grampoint(k)...     s = siegelz(g)...     if (-1)**k * s < 0:...         print k, g, s...126 282.4547208 -0.02762949482134 295.583907 -0.01690039157195 391.4482021 0.0232894207211 415.60146 0.3828895679232 446.8057559 -0.1410432574254 478.9568293 -0.0600779492288 527.6973416 -0.6654588176367 637.320354 0.1579920088377 650.8910448 0.8376010676379 653.5978317 0.217352745397 677.8523216 0.1342801302400 681.8765522 -0.06717029681461 762.6678783 0.1525747623507 822.4194896 0.7419389942518 836.5739092 -0.3982959071529 850.6795334 0.1756176097567 899.0502587 0.8296209263578 912.9534756 -0.3537611108595 934.3572317 0.121967239618 963.1605748 -0.04025206217626 973.1388241 -0.1578260824637 986.8259207 0.2852904149654 1007.905352 -0.5547383392668 1025.199826 -0.2608153866692 1054.715528 -0.07836895313694 1057.167832 -0.1696327251703 1068.189532 1.004550445715 1082.850818 0.1608682601728 1098.690513 -0.6875362205766 1144.741881 -0.1678183794777 1158.005614 9.646277741e-3793 1177.246581 0.7011668936795 1179.647457 0.3577242914807 1194.033228 0.7072513907819 1208.386043 0.2285739236848 1242.939633 -0.6161899199857 1253.626041 0.2131228775869 1267.84791 1.271208219887 1289.124631 0.1830773844964 1379.419269 -1.392200147992 1411.979146 -0.06783933868995 1415.459427 0.26234563541016 1439.777518 -0.93248917151028 1453.639612 -0.31061398741034 1460.561547 -0.65018202151043 1470.933184 0.09874838311046 1474.387414 -0.56801951111071 1503.115534 0.15303108571086 1520.304266 -0.1060373061094 1529.457103 -0.18396543921111 1548.873932 0.021029602581113 1551.155349 0.80035104271135 1576.211048 0.22764199421156 1600.060766 -0.029255812531165 1610.262397 8.575848564e-31178 1624.977566 -0.26558850131207 1657.717899 0.86456039181209 1659.971552 0.366254451231 1684.725787 0.61564240591250 1706.052177 -0.82859665581263 1720.616514 0.82880439981276 1735.158918 -0.027341480291290 1750.795762 -0.19532693341294 1755.258868 -0.26756065031307 1769.750095 0.51636718971319 1783.107961 0.092168807161329 1794.225992 0.30121171791331 1796.448134 0.13199414851342 1808.661254 -0.12585064951344 1810.880254 -0.030613045671402 1875.026023 -0.57585758121430 1905.854681 -1.1847957621443 1920.138276 1.2948050861456 1934.403327 -0.034364044741485 1966.159589 0.31926886461487 1968.346369 0.8645091781495 1977.089273 0.14639222511498 1980.366127 -0.014645964411513 1996.736314 0.25709447381517 2001.097757 0.82470754471532 2017.438527 -0.16917553921543 2029.407189 0.18074419731545 2031.581995 1.2073826881600 2091.233527 -0.039953658821613 2105.289905 0.46952288031620 2112.852035 -0.30587659711643 2137.666434 0.86927607261646 2140.899441 -0.26798368041656 2151.670095 -0.42517029231669 2165.658156 0.92471640911672 2168.883971 -0.016017269171684 2181.779042 -0.14203855761702 2201.097281 -1.2749399571704 2203.241961 -0.25409911721722 2222.528116 -0.38045431791744 2246.06145 -0.15318219771747 2249.267282 1.3662921871760 2263.150267 -0.73497348151773 2277.0188 0.48178311991787 2291.938132 0.069197120651796 2301.520437 -0.16879932991804 2310.032371 -0.6873684151816 2322.790332 -0.13773543181819 2325.977969 0.16457315071843 2351.452601 0.010042270961850 2358.873909 -0.100899971856 2365.231898 -0.19521299431863 2372.645911 1.1203447031876 2386.404453 -0.13973863021892 2403.31974 -0.62200709341902 2413.881627 -0.3678360151921 2433.927875 0.010430431051933 2446.574386 1.4138782851935 2448.681071 1.8816398131953 2467.627613 0.050886371841969 2484.448555 0.29653725341982 2498.101553 -0.3858252138

(Note: the Z values are not quite accurate to the indicated precision, due to the magnitude of the Gram points. Increasing the precision, one finds siegelz(grampoint(1982)) = -0.385825235416...)

For those interested in this topic, there is a very nice book called “Riemann’s Zeta Function” by H. M. Edwards. There is just as much information to be found on the web, e.g. on the amazing “Numbers, constants and computation” site by Xavier Gourdon and Pascal Sebah.