The First Computer Program

February 8, 2011

Here is our version of the Countess’ program:

(define (bernoulli limit) (let loop ((n 1) (bs '(-1/2)) (ns '()) (d 1) (ds '())) (if (= limit n) bs (let* ((2n (* n 2)) (2n-1 (- 2n 1)) (2n+1 (+ 2n 1)) (ns (cons 2n (map (lambda (x) (* 2n 2n-1 x)) ns))) (b (- (apply + (map * bs (cons 2n-1 (but-last ns)) (map / (cons 2n+1 ds)))))) (bs (append bs (list b))) (d (* d 2n 2n-1)) (ds (append ds (list d)))) (loop (+ n 1) bs ns d ds)))))

Here, n is as in the formula, ns, ds and bs are the accumulating lists of numerators, denominators, and Bernoulli numbers, and 2n, d, and b are the next elements of ns, ds and bs, respectively. Note that the order in which the elements are computed in the let* is important.

We used an auxiliary function but-last that is like cdr except that it operates at the end of the list instead of the beginning; but-last is useful from time to time, but not frequently enough to warrant inclusion in the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/USXM9t9Y, where you can also see the source code for but-last.

Posted by programmingpraxis

Filed in Exercises

5 Comments »

5 Responses to “The First Computer Program”

Michael said
February 8, 2011 at 9:19 AM
This is a great programming assignment, and I’ll enjoy trying it, but I thought I’d mention that Michael Greene is the current holder of the Lucasian Chair, Hawking stepped down in 2009.
Michael said
February 8, 2011 at 9:21 AM
Oops, Green, not Greene. Silly typo =(

Graham said

February 9, 2011 at 3:23 AM

My Python solution.
The first version is a translation of the provided solution into scheme. The
second is a memoized version of the recursive relation that for B_n that can be
found in a bunch of places (I got mine from p. 180 of the “CRC Handbook of
Discrete and Combinatorial Mathematics”). I had to throw in the filter statement
due to differing ideas about the proper enumeration of the Bernoulli numbers.
Note: requires the fractions module for rational numbers.

#!/usr/bin/env python

import functools
from fractions import Fraction


def b_list1(n):
    bs, ns, ds, d = [Fraction(-1, 2)], [], [], 1
    for k in xrange(1, n):
        x, y, z = 2 * k, 2 * k - 1, 2 * k + 1
        ns = [x] + [j * x * y for j in ns]
        b = -sum(p * q * r  for (p, q, r) in
                zip(bs, [y] + ns[:-1], [1 / j for j in ([Fraction(z)] + ds)]))
        bs.append(b)
        d *= Fraction(x * y, 1)
        ds.append(d)
    return bs


def memoize(func):
    func.memo = {}
    def memoizer(*args):
        try:
# Try using the memo dict, or else update it
            return func.memo[args]
        except KeyError:
            func.memo[args] = result = func(*args)
            return result
    return functools.update_wrapper(memoizer, func)


def c(m, k):
    n = d = 1
    for j in xrange(1, m - k + 1):
        n *= k + j
        d *= j
    return n / d


@memoize
def b_rec(n):
    if n == 0:
        return Fraction(1)
    else:
        return -sum(c(n + 1, j) * b_rec(j) for j in xrange(n)) / (n + 1)


def b_list2(n):
    return filter(lambda x: x != 0, (b_rec(k) for k in xrange(2 * n)))


if __name__ == "__main__":
    for b in b_list1(10):
        print b

    print "-" * 17
    for b in b_list2(10):
        print b

Jan said

February 13, 2011 at 11:19 PM

A solution in Common Lisp.

#|
0 = -1/2 (2n-1)/(2n+1) + B_1 (2n/2)
                       + B_3 (2n (2n-1) (2n-2) / (2 3 4))
                       + B_5 (2n (2n-1) ... (2n-4)/(2 3 4 5 6)
                       ...
                       + B_{2n-1}

n = 1, 0 = -1/2 1/3 + B_1
n = 2, 0 = -1/2 (3/5) + B_1 (4/2) + B_3
n = 3, 0 = -1/2 (5/7) + B_1 (6/2) + B_3 ((6 5 4)/(2 3 4)) + B_5

B_1 = -1/2
B_3 = 1/6
B_5 = -1/30
B_7 = 1/42
|#

(use-package :iterate)

(defun dot (xs ys)
  (iter (for x in xs) (for y in ys) (summing (* x y))))

(defun cumprod (xs)
  (iter (for x in xs)
        (for p first x then (* p x))
        (collect p)))

(defun odd-elements (xs)
  (iter (for x in xs)
        (for i from 1)
        (when (oddp i) (collect x))))

(defun factors (n)
  (iter (with 2n = (* 2 n))
        (for i from 2n downto 2)
        (for j from 2 to 2n)
        (collect (/ i j))))

(defun coeffs (n)
  (cons 1 (cumprod (factors n))))

(defun bernoulli-from-previous (n B0 Bs)
  (let* ((2n (* 2 n))
         (A0 (* B0 (/ (- 2n 1) (+ 2n 1))))
         (As (odd-elements (rest (coeffs n)))))
    (- (- A0) (dot As Bs))))

(defun bernoulli (nmax &optional (B0 -1/2))
  (cons B0 (iter (for n from 1 to nmax)
                 (collect (bernoulli-from-previous n B0 Bs) into Bs)
                 (finally (return Bs)))))

Jan Van lent said

February 13, 2011 at 11:28 PM

A solution in Haskell.

import Data.Ratio

-- construct Ratio of arbitrary precision Integers from pair of machine Ints
r :: Int -> Int -> Ratio Integer
r i j = (fromIntegral i) % (fromIntegral j)

-- return the odd elements of a list
-- odds [0,1,2,3,4] == [1,3]
odds [] = []
odds [x] = []
odds (x:y:xs) = y : odds xs

cumprod = scanl (*) 1

dot xs ys = sum (zipWith (*) xs ys)

ratios :: Int -> [ Ratio Integer ]
ratios n = zipWith r [2*n,2*n-1..3] [2..2*n-1]

coeffs = odds . cumprod . ratios

coeff0 n = (r (-1) 2) * (r (2*n-1) (2*n+1))

bernoulli :: [ Ratio Integer ]
bernoulli =  [ -(coeff0 n) - (dot (coeffs n) (tail bernoulli)) | n <- [0..] ]

test = take 9 bernoulli

Programming Praxis