## 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.

Pages: 1 2

### 4 Responses to “Sums Of Powers”

1. Graham said

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.

```#!/usr/bin/env python

import numpy as np
from fractions import Fraction
from operator import mul

def a_t(n):
A = np.zeros(shape=(n + 1, n + 1), dtype=Fraction)
for j in xrange(n + 1):
A[0, j] = Fraction(1, j + 1)
for i in xrange(1, n + 1):
for j in xrange(n - i + 1):
A[i, j] = (j + 1) * (A[i - 1, j] - A[i - 1, j + 1])
return A

def b(n):
return a_t(n)[n]

def c(n, k):
num = reduce(mul, xrange(n - k + 1, n + 1), 1)
denom = reduce(mul, xrange(1, k + 1), 1)
return num / denom

def s_naive(m, n):
return sum(pow(k, m) for k in xrange(1, n + 1))

def s(m, n):
s = sum(c(m + 1, k) * b(k) * pow(n, m + 1 - k) for k in xrange(m + 1))
return s // (m + 1)

if __name__ == "__main__":
print b(18)
print a_t(18)[:, 0]
print s_naive(10, 3)
print s(10, 3)
print s_naive(10, 1000)
print s(10, 1000)
```
2. Graham said

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.

3. 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.

```import Data.Ratio

a :: (Integral a, Integral b) => a -> a -> Ratio b
a i j = iterate (\xs -> zipWith (*) [1..] \$ zipWith (-) xs (tail xs))
(map (1 %) [1..]) !! fromIntegral i !! fromIntegral j

bernoullis :: (Integral a, Integral b) => a -> [Ratio b]
bernoullis upto = map (flip a 0) [0..upto]

choose :: Integral a => a -> a -> Ratio a
choose n k = product [1..n] % (product [1..n-k] * product [1..k])

s :: Integral a => a -> a -> Ratio a
s m n = 1 % (m+1) * sum [choose (m+1) k * a k 0 * (n%1)^(m+1-k) | k <- [0..m]]
```
4. […] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli […]