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

Pages: 1 2

### 5 Responses to “The First Computer Program”

1. Michael said

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.

2. Michael said

Oops, Green, not Greene. Silly typo =(

3. Graham said

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
```
4. Jan said

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)))))
```
5. Jan Van lent said

```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
```