The First Computer Program

February 8, 2011

By Mirko Tobias Schaefer. Licensed under CC-BY-2.0. Beginning in the 1820s, and continuing until his death in 1871, Charles Babbage, an English mathematician (he was for ten years the Lucasian Professor of Mathematics at Cambridge, a chair once held by Sir Isaac Newton and currently held by Stephen Hawking), engineer, scientist and inventor was fascinated with the idea of a machine that could perform arithmetic calculations. His first machine, the Difference Engine, first described in 1822, used finite differences to compute the values of polynomials. His second machine, the Analytical Engine, first described in 1837, was much more ambitious, consisting of a control unit, an arithmetic processor, and a storage space for intermediate calculations, and could be programmed to perform specific calculations by use of punched cards. Both machines were mechanical, not electrical, powered by steam. Neither machine was built prior to his death, due to a lack of money and the expertise to manufacture the needed parts with the required precision.

Augusta Ada King (née Byron), Countess of Lovelace, daughter of the English poet Lord Byron and an amateur mathematician, was fascinated by Babbage’s plans for the Analytical Engine. In 1842, at Babbage’s request, she translated a manuscript about the Analytical Engine by the Italian mathematician Luigi Menabrea (who later became the Prime Minister of Italy), adding her own notes, which were more extensive than the original manuscript. Note G describes the use of the Analytical Engine to compute the Bernoulli numbers; it is complete and rigorous, and is now recognized as the first computer program, making the Countess the first computer programmer.

The Countess’ translated manuscript and notes, available from FourmiLab, make fascinating reading, a terrific way to pass a cold winter evening. She was clearly in full command of her subject, and anticipated much of computer architecture and modern computer programming languages a full century before they were discovered. For instance, if you ignore the stilted language, this quote could have been written last week:

In studying the action of the Analytical Engine, we find that the peculiar and independent nature of the considerations which in all mathematical analysis belong to operations, as distinguished from the objects operated upon and from the results of the operations performed upon those objects, is very strikingly defined and separated.

It is hard to overstate the Countess' accomplishment. She had no high-level programming language, no libraries, no integrated development environment, not even an assembler or a debugger; she was working on "bare metal" — literally! The only documentation she had was that which she wrote herself. And she was working in untrodden territory, as there were no textbooks about programming, no tutors, no "Learn X in 21 Days" to help her. She didn't even have a machine to run her program. It is amazing that she wrote the program at all; incredibly, modern computer scientists have examined the program and declared it bug-free, and today’s programmers must simply stand in awe. The Countess was not only the first computer programmer, she was also an extremely talented one.

The program that the Countess chose to demonstrate the Analytical Engine was a program to compute the Bernoulli numbers and is based on Equation 8 of Note G:

$0 = -\frac{1}{2} \cdot \frac{2n-1}{2n+1}$ $+ \mathrm{B}_1\left(\frac{2n}{2}\right)$ $+ \mathrm{B}_3\left(\frac{2n\cdot(2n-1)\cdot(2n-2)}{2\cdot3\cdot4}\right)$ $+ \mathrm{B}_5\left(\frac{2n\cdot(2n-1)\cdot\ldots\cdot(2n-4)}{2\cdot3\cdot4\cdot5\cdot6}\right)$ $+ \ldots$ $+ \mathrm{B}_{2n-1}$

Thus, B₁ = -1/2 when n = 1, B₃ = 1/6 when n = 2, B₅ = -1/30 when n = 3, B₇ = 1/42 when n = 4, B₉ = -1/30 when n = 5, B₁₁ = 5/66 when n = 6, B₁₃ = -691/2730 when n = 7, B₁₅ = 7/6 when n = 8, B₁₇ = -3617/510 when n = 9, B₁₉ = 43867/798 when n = 10, and so on. The successive Bernoulli numbers are computed in order, each feeding the computation of the next, the various numerators and denominators being computed according to the pattern in the formula.

Your task is to write a program to calculate the Bernoulli numbers, using the method of the Countess’ Note G. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

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