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.
Programming Praxis 
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.
Oops, Green, not Greene. Silly typo =(
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
fractionsmodule 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 bA 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)))))A solution in Haskell.