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.

About these ads

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][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)
    
  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.

    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]]
    
  4. [...] today’s Programming Praxis exercise, our goal is to implement an algorithm that calculates the bernoulli [...]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 612 other followers

%d bloggers like this: