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.
#!/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][0] 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)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):
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]][…] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli […]