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 10¹². And the Riemann R approximation improves, percentage-wise, as x increases.

You can run the program at http://programmingpraxis.codepad.org/Fvwfmjc5.

Like this:

Be the first to like this post.

Pages: 1 2

Posted by programmingpraxis

Filed in Exercises

1 Comment »

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

Programming Praxis