Approximating Pi

July 29, 2011

In the two previous exercises we have seen three different methods for computing the prime-counting function π(n), all based on a formula of Legendre. We have also seen the difficulty of the calculation, as three eminent mathematicians got their sums wrong, and the difficulty continues even today with the use of computers, as Xavier Gourdon had to stop an attempt to compute π(10²³) because of an error. Given the difficulty of calculation, in some cases it might make sense to calculate an approximate value of π, instead of an exact value. In today’s exercise we look at two methods for approximating the prime-counting function.

The first method dates to the German mathematician Carl Friedrich Gauss, who gives an estimate π(x) ∼ Li(x), the logarithmic integral—this is the celebrated prime number theorem, which Gauss figured out when he was fifteen years old! The logarithmic integral can be calculated by expanding the infinite series

$\mathrm{Li}(x) = \int_0^x \frac{d\ t}{\log_e t} = \gamma + \log \log x + \sum_{k=1}^\infty \frac{(\log x)^k}{k!\ k}$

to the desired degree of accuracy; somewhere around a hundred terms ought to be sufficient for arguments up to 10¹². Beware the singularity at x=0.

In 1859, Bernhard Riemann, whose hypothesis about the prime numbers is one of the greatest open questions in mathematics, described a different, and much more accurate, approximation to the prime counting function:

$\mathrm{R}(x) = \sum_{n=1}^\infty \frac{\mu(n)}{n} \mathrm{Li}(x^{1/n})$

where μ(n) is the Möbius function that takes the value 1 when n is 1, 0 when the factorization of n contains a duplicated factor, and (−1)^k, where k is the number of factors in the factorization of n, otherwise. There is no better way to compute the Möbius function than to compute the factors of n. As with the logarithmic integral, about a hundred terms of the infinite series ought to be sufficient to get a good approximation for arguments up to 10¹².

Your task is to write functions to compute the logarithmic integral, the Möbius function, and Riemann’s R function, and to compare the results to the values you calculated for the prime-counting function in the previous exercise. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Posted by programmingpraxis

Filed in Exercises

4 Comments »

4 Responses to “Approximating Pi”

Graham said

August 1, 2011 at 2:51 PM

Here’s my Python solution; it reuses primes, a variation of the Python Cookbook’s version of the Sieve of Eratosthenes.

from math import log
from operator import mul

GAMMA = 0.57721566490153286061

def li(x):
    return GAMMA + log(log(x)) + sum(pow(log(x), n) /
                                     (reduce(mul, xrange(1, n + 1)) * n) for
                                      n in xrange(1, 101))

def factors(n):
    return (p for p in primes(n + 1) if n % p == 0)

def mu(n):
    m = lambda p: -1 if (n / p) % p else 0
    return reduce(mul, (m(p) for p in factors(n)), 1)

def r(x):
    return sum(mu(n) / n * li(pow(x, 1.0 / n)) for n in xrange(1, 101))

David said

August 25, 2012 at 3:36 AM

Clojure version, uses O’Neil sieve to get list of primes. Results for Riemann function matches table above.

(load "lazy-sieve")

(defn factors
   "list of factors of n"
   [n]
   (loop [n n, prime lazy-primes, l (vector-of :int)]
      (let [p (first prime)]
         (cond
            (= n 1)  l
            (= (mod n p) 0)
               (recur (/ n p) prime (conj l p))
            :otherwise
               (recur n (rest prime) l)))))

(defn calc-mu
   "Calculate mu(n) for use in building a static array, n > 1"
   [n]
   (let [fac (factors n)]
      (if (not= (count fac) (count (distinct fac)))
         0
         (reduce * (map (fn [x] -1) fac)))))

(def mu (into (vector-of :int 0 1)
   (map calc-mu (range 2 101))))

(defn powers [x]
   (iterate (partial * x) x))

(defn float-factorials []
   (map first (iterate (fn [[n,i]] [(* n i) (inc i)]) [1.0 2])))

(defn logint
   "Compute gauss' logarithmic integral to N terms"
   [n, x]
   (let [gamma 0.5772156649015329,
         ln-x  (Math/log x),
         sum   (reduce +
                  (map #(/ %1 (* %2 %3))
                     (take n (powers ln-x))
                     (take n (float-factorials))
                     (range 1 (inc n))))]
      (+ gamma (Math/log ln-x) sum)))

(defn riemann
   "Compute riemann function, max 100 terms"
   [x]
   (reduce +
      (map #(* (/ (mu %1) %1) (logint 100 (Math/pow x (/ 1 %1))))
         (range 1 101))))

Robert said
January 30, 2015 at 11:35 AM
You have definied series expansion for li(x).
Li(x) = li(x) – li(2).
li(x) ~ pi(x)

Lucas Brown said

February 22, 2016 at 10:23 AM

I’ll answer this question in a different manner than it was posed in, since the Riemann R-function can be rewritten in a much more numerically-tractable form as the Gram series.

def gramseries(x, n=None):
    # Evaluates the Riemann R-function at x, using n terms of the Gram series
    if n is None: n = 6 * int(log(x, 10)+1) # Appears to be a good value for 64-bit floats and x < 10**9 or so
    lnx = log(x)
    total, lnxp, kfac = 1, 1, 1
    for k in range(1, n+1):
        lnxp *= lnx
        kfac *= k
        total += lnxp / (k * kfac * riemannzeta(k+1))
    return total

This then requires implementing the Riemann zeta function. To compute this, I use the Dirichlet eta function to convert it to an alternating series and then use an elementary convergence-acceleration technique on the alternating series:

def altseriesaccel(a, n):
    # Convergence acceleration for an alternating series.
    # The absolute values of the terms are yielded by the iterable a.
    # The number of terms to use is n.
    d = (3 + sqrt(8))**(n)
    d = (d + 1/d) / 2
    b, c, s = -1, -d, 0
    for k in range(n):
        c = b - c
        s = s + c * next(a)
        b = (k + n) * (k - n) * b / ((k + 0.5) * (k + 1))
    return s / d

# Evaluates the Riemann zeta function at n > 1 by feeding k terms of the Dirichlet eta
# function into a technique for accelerating the convergence of alternating series.
# The default value of k=24 appears to be sufficient for 64-bit floats.
def riemannzeta(n, k=24): return altseriesaccel((1/j**n for j in count(1)), k+1) / (1-2**(1-n))

Programming Praxis