Counting Primes

Azimuth 2026-10-24

Since the probability that a number $n$ is prime is roughly $1/\log n,$ the number of primes $\le x$ is close to

$\displaystyle{ \mathrm{li}(x) = \int_0^x \frac{d s}{\log s} }$

Indeed, if $\pi(x)$ is the number of primes $\le x,$ the recently proved quasi-Riemann hypothesis implies

$|\mathrm{li}(x) - \pi(x)| \le C x^p$

for some constants $C > 0$ and $p < 1.$ If the Riemann hypothesis is true, any $p > 1/2$ will do.

But here’s what I want to talk about today. For all $x$ where $\pi(x)$ can actually be counted, $\pi(x) < \mathrm{li}(x).$ But in 1914, Littlewood shocked the world by sketching a proof that $\pi(x) - \mathrm{li}(x)$ changes sign infinitely often!

Littlewood’s proof didn’t say when $\pi(x)$ first exceeds $\mathrm{li}(x).$ He gave his student Stanley Skewes the task of figuring that out.

In 1933, Skewes gave an enormous upper bound for where this first happens, assuming the Riemann hypothesis:

$\displaystyle{ e^{e^{e^{79}}} \approx 10^{10^{10^{34}}} }$

This is called Skewes’ number. Hardy called it “the largest number which has ever served any definite purpose in mathematics”. In 1955, still working on his advisor’s assignment, Skewes gave an even larger upper bound not assuming the Riemann hypothesis:

$\displaystyle{ e^{e^{e^{e^{7.705}}}} \approx 10^{10^{10^{964}}} }$

This disgustingly large quantity is called the second Skewes number.

I first read about Skewes’ number as a kid in The Guinness Book of World Records. It always intrigued me, but I never tried to understand the math involved, because I figured it would be too hard. This article is my attempt to finally make some progress.

In 1966, Lehman found a trick that lowered Skewes’ bound to about $1.65 \times 10^{1165}.$ In 1987, te Riele brought it down to a mere $6.66 \times 10^{370}.$

In 2000, Bays and Hudson went further. They proved that there are at least $10^{153}$ integers $x$ near $1.39822 \times 10^{316}$ with $\pi(x) > \mathrm{li}(x).$ Their paper includes many fascinating graphs—but the printing is of very poor quality, unpleasant to look at. I think the world deserves better pictures of their amazing result, so yesterday I decided to redo their work with the help of AI.

To draw these graphs we need to use Riemann’s explicit formula for $\pi(x)$ as a sum over zeros of the Riemann zeta function. My calculations use twice as many zeros as Bays and Hudson had access to, but the computation now takes only a couple of seconds on a laptop. Here is my version of their Figure 2b:

$Narrow window: blue and red curves of (pi(x) − li(x)) log x / sqrt(2) for x from 1.3977 to 1.3989 times 10^316$

Bays and Hudson’s proof does not depend on the Riemann hypothesis, nor do the later proofs mentioned below. Their graphs are another matter—and the same is true for the graphs here. It’s hard to get a really precise approximation to $\pi(x)$ for huge numbers $x$ without using the Riemann hypothesis, since our best formula for it is a sum over zeros of the Riemann zeta function, and even a single zero off the critical line $\mathrm{Re}(z) = \frac{1}{2}$ could have a big effect. However, a smoothed version of $\pi(x)$ is much less sensitive to this problem.

The red curves in my graphs are smoothed. Since they go above zero for $10^{311}$ integers near $1.39822 \times 10^{316},$ and they go up sufficiently high, Bays and Hudson were able to prove that $\pi(x) > \mathrm{li}(x)$ for some $x$ near this number. On the other hand, the blue curves are not deliberately smoothed. Its calculation assumes the Riemann hypothesis, and could change substantially if that was wrong! There could also be other problems, as I’ll explain. But it’s probably close to correct if the Riemann hypothesis holds.

Intriguingly, $\pi(x) - \mathrm{li}(x)$ also becomes positive for some slightly smaller numbers, in a region around $1.39717 \times 10^{316}.$ Patrick Demichel spotted this region, Saouter and Demichel proved it exists in 2010, and Saouter, Trudgian and Demichel pinned it down more sharply in 2015. Here is a wider view that also includes their region. See the sharp spike in the blue curve?

$Wider view: the same curves for x from 1.3955 to 1.3995 times 10^316$

What is plotted

Nobody can compute $\pi(x)$ exactly for numbers so large as $10^{316}.$ What we can compute is an approximation coming from Riemann’s explicit formula, which expresses $\pi(x) - \mathrm{li}(x)$ as a sum over the nontrivial zeros $\rho = \tfrac12 + i\gamma$ of the zeta function. Assume the Riemann hypothesis, so $\gamma$ is real. Write $u = \log x$ and measure $\pi(x) - \mathrm{li}(x)$ in units of $\sqrt{x} / \log x$ :

$\begin{array}{ccl} D(u)\! &:= & \bigl(\pi(x) - \mathrm{li}(x)\bigr)\,\frac{\log x}{\sqrt{x}} \\ \\ &\approx & \!\!-\Bigl(1 + \frac{2}{u} + \frac{8}{u^2}\Bigr) - 2\sum_{0 <\gamma\le T} \mathrm{Re}\!\left[\frac{e^{i\gamma u}}{\rho}\Bigl(1 + \frac{1}{\rho u}\Bigr)\right] \end{array}$

The symbol $\approx$ hides several approximations of quite different kinds; the section “What is rigorous and what is not” explains exactly what it means. The term

$-(1 + 2/u + 8/u^2)$

comes from the squares of primes: $\pi(x)$ counts primes, but the explicit formula naturally counts prime powers, and correcting for squares gives this term. This term is why $\pi(x) - \mathrm{li}(x)$ is usually negative, by about $\sqrt{x} / \log x.$ Higher powers of primes give extra corrections that are negligible when $x$ is huge.

More interesting is the sum over zeros: each zero contributes an oscillation in $D(u)$ with frequency $\gamma$ and amplitude roughly $1/\gamma.$ These oscillations are the key: $\pi(x) - \mathrm{li}(x)$ becomes positive only when enough of these oscillations get in phase and add up to a positive number bigger than 1! Or more precisely, bigger than

$1 + 2/u + 8/u^2$

On the graphs, $D > 0$ means $\pi(x) > \mathrm{li}(x),$ and the vertical axis is in units of $\sqrt{x} / \log x,$ which is about $10^{155}$ near $x = 10^{316}.$ Each graph shows two curves:

• Blue: the sum over the first 2,001,052 zeros, with no smoothing. This approximates $D(u)$ point by point. Leaving out the higher zeros causes an error of roughly $\pm 0.002$ in these units, assuming the omitted terms behave like random noise.

• Red: the choice Bays and Hudson made. It uses the first $10^6$ zeros, each damped by a factor $\exp(-\gamma^2/2\alpha)$ with $\alpha = 10^{10}.$ This is a Gaussian average of $D(u)$ over a width of about $10^{-5}$ in $\log x,$ which is the quantity Lehman’s trick from 1966 controls rigorously.

Data and computation

The graphs here were made with a short Python program, skewes_window.py, using the NumPy, mpmath, Matplotlib and FINUFFT libraries. If you’re good at Python and/or number theory, please check this program!

The Riemann zeta zeros come from a table computed by Andrew Odlyzko: the first 2,001,052 zeros, accurate to within $4 \times 10^{-9}.$ Back in 2000, Bays and Hudson needed a personal favor from Odlyzko to get a million zeros. Today the same table is a free download, and it is even packaged on PyPI as passagemath-database-odlyzko-zeta.

Each point on any curve in the graphs here is a sum of one or two million sine waves. Summing directly at 8000 points would take about $10^{10}$ operations. Instead the program uses a nonuniform fast Fourier transform from the FINUFFT library, which finishes in about a second on one core. Two precautions matter: the program never forms $x$ itself, only $u = \log x,$ and it computes each phase $\gamma u$ modulo $2\pi$ carefully, because $\gamma u$ reaches about $10^9.$

Three graphs

Narrow window: $1.3977$ – $1.3989 \times 10^{316}.$ This graph shows our estimate of $(\pi(x) - \mathrm{li}(x)) \cdot \log x / \sqrt{x}$ for $x$ between $1.3977 \times 10^{316}$ and $1.3989 \times 10^{316}.$

$Narrow window: blue and red curves of (pi(x) − li(x)) log x / sqrt(x) for x from 1.3977 to 1.3989 times 10^316$

This is the same range as Figure 2b in Bays and Hudson’s paper. The region they proved positive is the green hump near $1.39822 \times 10^{316}.$ $D$ usually hovers around $-1,$ but throughout this window the low zeros happen to be nearly in phase, lifting it close to zero. Even so it is a narrow escape: the blue curve rises only about 0.017 above zero. The red curve peaks at about $+0.010,$ which is what lets Lehman’s theorem prove that $\pi(x) > \mathrm{li}(x)$ somewhere in this region. Two smaller bumps of the blue curve, near 1.3981 and 1.3983, barely cross zero; they are about as big as the truncation error, so they could be artifacts.

Narrower window: $1.3980$ – $1.3984 \times 10^{316}.$

This is just a piece of the previous graph, to make the green hump more visible.

Wider view: $1.3955$ – $1.3995 \times 10^{316}.$ This zooms out to show what lies just to the left.

$Wider view: the same curves for x from 1.3955 to 1.3995 times 10^316$

The thin spike at about $1.39717 \times 10^{316}$ is a second positive region, found after Bays and Hudson. Saouter and Demichel proved $\pi(x) > \mathrm{li}(x)$ there in 2010, and Saouter, Trudgian and Demichel later narrowed it to an interval centered at $\log x = 727.951335426.$ In this graph the blue curve rises only about 0.0018 above zero at that spike. That is no larger than the error from using only two million zeros, so this graph alone could not prove the spike is positive. Their proof used about $5 \times 10^8$ zeros. The red curve, smoothed too heavily to see so narrow a spike, stays below zero there.

What is rigorous and what is not

The red curve could be given rigorous error bars; the blue curve cannot, and $\pi(x) - \mathrm{li}(x)$ itself, point by point, is out of reach near $10^{316}.$

Where the formula comes from. The starting point is Riemann’s explicit formula for his prime-power counting function $J(x) = \sum_{p^k \le x} 1/k.$ For $x > 1$ it says

$J(x) = \mathrm{li}(x) - \sum_\rho \mathrm{li}(x^\rho) - \log 2 + \int_x^\infty \frac{d t}{t(t^2-1)\log t}$

where the sum over zeros converges only conditionally and must be taken in order of increasing $|\gamma|.$ Proofs are in the books of Edwards and Ingham. Möbius inversion gives

$\displaystyle{ \pi(x) = \sum_n \mu(n) J(x^{1/n})/n, }$

$\pi(x) - \mathrm{li}(x) = -\tfrac12 \mathrm{li}(\sqrt{x}) - \sum_\rho \mathrm{li}(x^\rho) + (\text{smaller terms}).$

Expanding

$\mathrm{li}(\sqrt{x}) = (2\sqrt{x}/u)(1 + 2/u + 8/u^2 + \cdots)$

and

$\mathrm{li}(x^\rho) = \mathrm{Ei}(\rho u) = (x^\rho/(\rho u))(1 + 1/(\rho u) + \cdots),$

multiplying by $u/\sqrt{x},$ and pairing each zero with its complex conjugate gives the formula for $D(u).$ Apart from the $1/(\rho u)$ correction and the Gaussian damping, this is the function $g(x)$ in Table 1 of Bays and Hudson.

What the $\approx$ means. Three separate approximations are involved, and they are very different in character.

1) Dropped smaller terms. These come from $x^{1/3}$ and higher roots, the constants, and the zero sum belonging to $J(\sqrt{x}).$ Assuming the Riemann hypothesis, they are smaller than the terms kept by factors like $x^{-1/4} \log^2 x,$ which is about $10^{-76}$ near $10^{316}.$ Rigorous bounds for them are routine.

2) Truncated asymptotic expansions. The first omitted prime-square term is about $48/u^3 \approx 10^{-7}.$ The error in expanding each $\mathrm{li}(x^\rho)$ is of order $1/(\gamma^3 u^2),$ negligible even when summed over all zeros. This step too can be made rigorous with explicit remainder bounds.

3) Stopping the sum at the 2,001,052nd zero. This is the step that could cause trouble. The best known rigorous bounds for truncating Riemann’s explicit formula at zeros of height $T$ give an error of order $x \log^2(xT)/T.$ Making that smaller than $\sqrt{x}/\log x$ would need $T$ of roughly $\sqrt{x} \log^2 x,$ around $10^{158}$ here, while the 2,001,052nd zero only has height

What can we do about the problem in item 3? Assuming the Riemann hypothesis, $D(u)$ is an almost periodic function of $u,$ which is the setting Rubinstein and Sarnak use. A Parseval identity for such functions should then give the long-run mean square of the omitted tail as about $\sum_{\gamma > T} 2/|\rho|^2 \approx (0.002)^2.$ That is a statement about averages over very long stretches of $u.$ It does not bound the error at any particular point, or even the average error across our window. But we can say $\pm 0.002$ is the typical error at each point of the blue curve if the omitted terms in the sum over zeta zeros act like random noise. This is just a heuristic, though it’s been checked at $x = 10^{20}.$ So the blue curve is just a good guess, not something proved.

The smoothed curve. Every positive region found so far was proved using a theorem of Lehman (1966), stated as Theorem 1 in Bays and Hudson’s paper. Instead of trying to pin down $\pi(x) - \mathrm{li}(x)$ at a single point, Lehman’s theorem handles a weighted average of it: the average of $D(u)$ against a Gaussian bump of width $1/\sqrt{\alpha}$ centred at some point $\omega.$ The theorem says this average equals $-1$ minus a sum over the zeros, with each zero damped by a factor $e^{-\gamma^2/2\alpha},$ plus an error term with explicit bounds that hold for every $\omega.$ The red curve is essentially this damped sum, plotted as a function of $\omega.$ The key point is that an average can only be positive if the function being averaged is positive somewhere. So if the damped sum is above zero by more than the largest possible error at some $\omega,$ then $\pi(x) > \mathrm{li}(x)$ for some $x$ with $\log x$ near $\omega.$ Bays and Hudson proved their result this way, and later authors used refined versions of the same theorem.

The smoothing is what gives this approach its teeth. The damping makes the contribution of the omitted zeros decay like $e^{-T^2/2\alpha},$ so stopping at height $T$ costs almost nothing. With these parameters the total error is a few thousandths. So one could draw an upper and a lower curve that provably bracket the smoothed function. Wherever the lower curve is above zero, $\pi(x) > \mathrm{li}(x)$ somewhere nearby.

This holds without assuming the Riemann hypothesis. Lehman’s theorem needs the Riemann hypothesis only up to a finite height, which has been verified: Platt and Trudgian have now checked it up to height $3 \times 10^{12}.$ Any hypothetical zeros off the critical line beyond that height are covered by an extra error term, which the Gaussian smoothing makes harmless. For Bays and Hudson’s parameters this term is smaller than $10^{-2001}.$ So the positive regions near $1.39822 \times 10^{316}$ and $1.39717 \times 10^{316}$ are established unconditionally.

The method is one-sided, though. It can prove that crossings exist, but a negative average says nothing about a narrow positive spike hiding inside it.

Earlier crossings. Bays and Hudson observed (p. 1287) that current knowledge cannot rule out sign changes of $\pi(x) - \mathrm{li}(x)$ over gaps as large as $10^{300}$ near $10^{316}.$ For the same reason, nobody can currently rule out a crossing somewhere below $1.397 \times 10^{316}.$ Bays and Hudson thought none is likely. Their Figure 1 marks near-misses around $10^{176},$ $10^{179},$ $10^{190},$ $10^{260}$ and $10^{298}.$ Rigorously, it is known only that $\pi(x) < \mathrm{li}(x)$ for all $x \le 10^{19}.$ So between $10^{19}$ and about $1.397 \times 10^{316}$ the question remains open!

In 2025, Michel and Rivers listed some places in this enormous interval where $\pi(x)$ might get bigger than $\mathrm{li}(x).$ If you want to contribute to this fascinating story, and you enjoy heavy-duty computations, you could look at those places. You may need a big table of Riemann zeta zeros.

References

The main papers

• Bays, Carter, and Richard H. Hudson. A new bound for the smallest $x$ with $\pi(x) > \mathrm{li}(x).$ Mathematics of Computation 69, no. 231 (2000): 1285–96. https://doi.org/10.1090/S0025-5718-99-01104-7.

• Lehman, R. Sherman. On the difference $\pi(x) - \mathrm{li}(x).$ Acta Arithmetica 11 (1966): 397–410.

• te Riele, Herman J. J. On the sign of the difference $\pi(x) - \mathrm{li}(x).$ Mathematics of Computation 48, no. 177 (1987): 323–28. https://www.ams.org/mcom/1987-48-177/S0025-5718-1987-0866118-6/.

• Saouter, Yannick, and Patrick Demichel. A sharp region where $\pi(x) - \mathrm{li}(x)$ is positive. Mathematics of Computation 79, no. 272 (2010): 2395–2405. https://www.ams.org/journals/mcom/2010-79-272/S0025-5718-10-02351-3. Also available at https://hal.archives-ouvertes.fr/hal-00515378.

• Saouter, Yannick, Timothy Trudgian, and Patrick Demichel. A still sharper region where $\pi(x) - \mathrm{li}(x)$ is positive. Mathematics of Computation 84, no. 295 (2015): 2433–46. https://www.ams.org/journals/mcom/2015-84-295/S0025-5718-2015-02930-5. Also available at https://hal.archives-ouvertes.fr/hal-01174411.

• Stoll, Douglas A., and Patrick Demichel. The impact of $\zeta(s)$ complex zeros on $\pi(x)$ for $x < 10^{10^{13}}.$ Mathematics of Computation 80, no. 276 (2011). https://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02477-4. • Platt, Dave, and Tim Trudgian. The Riemann hypothesis is true up to $3 \cdot 10^{12}.$ Bulletin of the London Mathematical Society (2021). Also available at arXiv:2004.09765.

• Revers, Michael. New bounds in R.S. Lehman’s estimates for the difference $\pi(x) - \mathrm{li}(x).$ (2025). Available at arXiv:2501.04488.

The explicit formula and its approximations

• Davenport, Harold. Multiplicative Number Theory. 3rd ed. Revised by Hugh L. Montgomery. New York: Springer, 2000. Ch. 17.

• Edwards, H. M. Riemann’s Zeta Function. New York: Academic Press, 1974.

• Ingham, A. E. The Distribution of Prime Numbers. Cambridge: Cambridge University Press, 1932. Reprint, 1990.

• Rubinstein, Michael, and Peter Sarnak. Chebyshev’s bias. Experimental Mathematics 3, no. 3 (1994): 173–97.

Historical background

• Littlewood, J. E. Sur la distribution des nombres premiers. Comptes Rendus de l’Académie des Sciences, Paris 158 (1914): 1869–72.

• Skewes, S. On the difference $\pi(x) - \mathrm{li}(x).$ Journal of the London Mathematical Society 8, no. 4 (1933): 277–83. https://doi.org/10.1112/jlms/s1-8.4.277.

• Skewes, S. On the difference $\pi(x) - \mathrm{li}(x),$ II. Proceedings of the London Mathematical Society, 3rd ser., 5, no. 1 (1955): 48–70. https://doi.org/10.1112/plms/s3-5.1.48.

• Skewes’s number. Wikipedia.

Data and software

• Barnett, Alex H., Jeremy F. Magland, and Ludvig af Klinteberg. A parallel nonuniform fast Fourier transform library based on an “exponential of semicircle” kernel. SIAM Journal on Scientific Computing 41, no. 5 (2019): C479–C504. Also available at arXiv:1808.06736.

• Odlyzko, Andrew. Tables of zeros of the Riemann zeta function. https://www-users.cse.umn.edu/~odlyzko/zeta_tables. Also available as https://pypi.org/project/passagemath-database-odlyzko-zeta/.