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.

About these ads

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

    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
    

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 627 other followers

%d bloggers like this: