Sums Of Powers
February 11, 2011
Our first computation of the Bernoulli numbers is a simple transliteration of the given algorithm:
(define (a i j)
(if (zero? i) (/ (+ j 1))
(* (+ j 1) (- (a (- i 1) j) (a (- i 1) (+ j 1))))))
(define (b n) (a n 0))
That works, but its exponential running time renders it unsuitable for use. Here’s a dynamic programming version that builds the array in such a way that all intermediate calculations are available as they are needed, one row at a time:
(define (akiyama-tanigawa limit)
(let loop ((i 0) (bs '()) (xs (map / (range 1 (+ limit 2)))))
(if (null? (cdr xs)) (reverse bs)
(loop (+ i 1) (cons (car xs) bs)
(map * (range 1 (- limit i -1))
(map - (but-last xs) (cdr xs)))))))
Here are some examples:
> (b 19)
43867/798
> (akiyama-tanigawa 19)
(1 1/2 1/6 0 -1/30 0 1/42 0 -1/30 0 5/66 0 -691/2730 0 7/6 0 -3617/510 0 43867/798)
We also have two choices for the function that computes the sums of powers. The brute-force version of the function is O(n):
(define (sum-powers m n)
(do ((i 1 (+ i 1)) (s 0 (+ s (expt i m)))) ((< n i) s))))
The version that uses Bernoulli’s formula is O(m), and requires an auxiliary function choose
that calculates the binomial coefficient:
(define (choose n r)
(if (zero? k) 1
(* n (/ k) (choose (- n 1) (- k 1)))))
(define (bernoulli-formula m n)
(let ((m+1 (+ m 1)))
(let loop ((k 0) (bs (akiyama-tanigawa m+1)) (s 0))
(if (< m k) (/ s m+1)
(loop (+ k 1) (cdr bs)
(+ s (* (choose m+1 k) (car bs)
(expt n (- m+1 k)))))))))
Here are some examples:
> (sum-powers 10 1000)
91409924241424243424241924242500
> (bernoulli-formula 10 1000)
91409924241424243424241924242500
Bernoulli’s formula is strikingly better than the brute-force calculation when n is large. For instance, on my computer (sum-powers 100 1000000)
takes about 12 seconds, but (bernoulli-formula 100 1000000)
makes the same calculation in 150 milliseconds, about two orders of magnitude better. That’s how Jakob Bernoulli was able to calculate S10(1000) in half of a quarter of an hour.
We used but-last
from the previous exercise and range
from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/D2FTIMnB.
I used this as an opportunity to play with NumPy’s (http://numpy.scipy.org/)
n-dimensional arrays. Also uses the fractions module for rational numbers.
Small typo that’s bothering me. It’s equivalent, but line 17 should read
return a_t(n)[n, 0]
instead to be consistent with the rest.
I believe the formula given for s in the exercise is wrong: it should be sum of k = 1 to m + 1 (or 0 to m), not 1 to n. At least that’s the version I found on wikipedia that worked for me after the 1 to n one didn’t.
My Haskell solution (see http://bonsaicode.wordpress.com/2011/02/11/programming-praxis-sums-of-powers/ for a version with comments):
[…] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli […]