## 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 *S*_{10}(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 […]