Approximating Pi

July 29, 2011

Here is our version of the logarithmic integral:

(define (logint x)
  (let ((gamma 0.57721566490153286061) (log-x (log x)))
    (let loop ((k 1) (fact 1) (num log-x)
               (sum (+ gamma (log log-x) log-x)))
      (if (< 100 k) sum
        (let* ((k (+ k 1))
               (fact (* fact k))
               (num (* num log-x))
               (sum (+ sum (/ num fact k))))
          (loop k fact num sum))))))

We need to be able to factor integers to compute the Möbius function. For general usage, you may want to choose some advanced factoring function, but since we only need to factors the integers to a hundred, a simple factorization by trial division suffices:

(define (factors n)
  (let loop ((n n) (fs (list)))
    (if (even? n) (loop (/ n 2) (cons 2 fs))
      (if (= n 1) (if (null? fs) (list 1) fs)
        (let loop ((n n) (f 3) (fs fs))
          (cond ((< n (* f f)) (reverse (cons n fs)))
                ((zero? (modulo n f))
                  (loop (/ n f) f (cons f fs)))
                (else (loop n (+ f 2) fs))))))))

The Möbius function runs down the list of factors, returning 0 if any factor is the same as its predecessor or ±1 depending on the count:

(define (mobius-mu n)
  (if (= n 1) 1
    (let loop ((fs (factors n)) (prev 0) (m 1))
      (cond ((null? fs) m)
            ((= (car fs) prev) 0)
            (else (loop (cdr fs) (car fs) (- m)))))))

The Riemann function computes a list of Möbius numbers once, when the function is instantiated, then runs through the list accumulating the sum:

(define riemann-r
  (let ((ms (let loop ((n 1) (k 100) (ms (list)))
              (if (zero? k) (reverse ms)
                (let ((m (mobius-mu n)))
                  (if (zero? m) (loop (+ n 1) k ms)
                    (loop (+ n 1) (- k 1) (cons (* m n) ms))))))))
    (lambda (x)
      (let loop ((ms ms) (sum 0))
        (if (null? ms) sum
          (let* ((m (car ms)) (m-abs (abs m)) (m-recip (/ m)))
            (loop (cdr ms) (+ sum (* m-recip (logint (expt x (/ m-abs))))))))))))

Now for the payoff. The table below shows the differences between π(x) and the values of the logarithmic integral and Riemann’s function at various values of x:

           pi     li-pi  r-pi
      ----------- ----- -----
10^3          168    10     0
10^6        78498   130    29
10^9     50847534  1701   -79
10^12 37607912018 38263 -1476

The logarithmic integral approximation is good, but the Riemann R approximation is stunning; the leading seven digits are right for 1012. And the Riemann R approximation improves, percentage-wise, as x increases.

You can run the program at


Pages: 1 2

4 Responses to “Approximating Pi”

  1. Graham said

    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))
  2. David said

    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"
       (loop [n n, prime lazy-primes, l (vector-of :int)]
          (let [p (first prime)]
                (= n 1)  l
                (= (mod n p) 0)
                   (recur (/ n p) prime (conj l p))
                   (recur n (rest prime) l)))))
    (defn calc-mu
       "Calculate mu(n) for use in building a static array, n > 1"
       (let [fac (factors n)]
          (if (not= (count fac) (count (distinct fac)))
             (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"
       (reduce +
          (map #(* (/ (mu %1) %1) (logint 100 (Math/pow x (/ 1 %1))))
             (range 1 101))))
  3. Robert said

    You have definied series expansion for li(x).
    Li(x) = li(x) – li(2).
    li(x) ~ pi(x)

  4. Lucas Brown said

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

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: