Enigmatic Code

Programming Enigma Puzzles

Tag Archives: babbage

Running the first program: Part 3

[ Part 1 | Part 2 | Part 3 ]

[Note: The programs can be run without the need to install anything @repl.it]

Following on from my previous articles on programming the Analytical Engine (see Part 1 and Part 2), in this post we will actually translate Ada Lovelace’s diagram that is considered to be the first published program into a directly executable program for the Analytical Engine emulator, where each row in the diagram corresponds directly to a statement in the assembly language program.

As suggested in Part 2, I have written a simple assembler to make it a bit easier to write programs for the Analytical Engine emulator. It allows you to specify the arithmetic operations for the Engine as a single statement (rather than the usual 4 cards), and it keeps track of symbolic labels and automatically computes branch offsets for you.

The assembler is implemented by the assemble() function, which I’ve included in the latest version of analytical_engine.py.

The assemble() function takes a multi-line string representing the program, and returns a list of cards suitable for use as a program in the AnalyticalEngine class. It also returns a map from the labels used in the assembler program to the index of the corresponding card in the program.

The factorial program from Part 2 would look like this using the assembler:

from analytical_engine import AnalyticalEngine, Column

n = 40

# initialise the engine
ae = AnalyticalEngine(vars=3, number=Column(digits=50), trace=1)

# assemble the program to compute factorial(n)
(program, _) = ae.assemble("""
  :init
  SET 0 <- {n}
  SET 1 <- 1
  SET 2 <- 1
  :loop
  # operation 1: v[2] = v[0] * v[2]
  MUL 0 2 -> 2
  # operation 2: v[0] = v[0] - 1
  SUB 0 1 -> 0
  # branch if non-zero to operation 1
  BRN loop
  # end
  HALT
""".format(n=n))

# load the program
ae.load_program(program)

# run the program
ae.run()

# the result is in v[2]
print("factorial({n}) = {r}".format(n=n, r=ae.v[2]))

[Program 6 - factorial3.py]

The embedded program for the Analytical Engine is highlighted.

Note that I have used Python’s format() function to interpolate the (Python) variable n into the string before it is passed to the assembler. So the assembler just sees “… SET 0 <- 40 …”.

By setting the trace=1 flag on the emulator the actual program cards loaded into the emulator are printed out before execution (as well as the trace of the execution itself).

The format that the assembler accepts should be fairly self-explanatory, from the example above, and from the implementation of Ada’s program given below, but the main points are:

  1. Comments (as in Perl and Python) are indicated by a # (hash) and continue to the end of the line.
  2. A line that starts with a : (colon) declares a label (which can be used as the target of a branch).
  3. Statements start with an op-code, followed by a number of space separated arguments.
  4. The multiple cards required for an arithmetical operation (ADD, SUB, MUL, DIV) are replaced by a single statement. Each takes two arguments to indicate the source of the operands: either the index of a variable in the store, or the special token DATA, which indicates that the value is loaded from the input data stack. Additional arguments are used to indicate indices in the store where the result is to be stored.
  5. The offset in branch statements can be specified using a symbolic label, by using a symbol that is declared as a label elsewhere in the program. (The label can appear either before or after the branch statement).
  6. The arrows (<- and ->) used in the SET statement to indicate the value being assigned to a variable in the store, and in arithmetic operations to indicate the variables in the store that the result of the operation is assigned to, are entirely optional. (They are removed during the parsing process, but I think it makes the program a little more readable).

So we can now write Ada’s program for Bernoulli Numbers like this:

from analytical_engine import AnalyticalEngine, Column
from enigma import raw_input, printf

# initialise the engine
ae = AnalyticalEngine(vars=14, number=Column(digits=10, dp=40), trace=1)
# assemble the program
(program, labels) = ae.assemble("""
  :init
  SET 0 <- 0
  SET 1 <- 1
  SET 2 <- 2
  SET 3 <- 1
  :start
  MUL 2 3 -> 4 5 6
  SUB 4 1 -> 4
  ADD 5 1 -> 5
  DIV 4 5 -> 11
  DIV 11 2 -> 11
  SUB 13 11 -> 13
  SUB 3 1 -> 10
  BRZ finish
  ADD 2 7 -> 7
  DIV 6 7 -> 11
  MUL DATA 11 -> 12
  ADD 12 13 -> 13
  SUB 10 1 -> 10
  BRZ finish
  :loop
  SUB 6 1 -> 6
  ADD 1 7 -> 7
  DIV 6 7 -> 8
  MUL 8 11 -> 11
  SUB 6 1 -> 6
  ADD 1 7 -> 7
  DIV 6 7 -> 9
  MUL 9 11 -> 11
  MUL DATA 11 -> 12
  ADD 12 13 -> 13
  SUB 10 1 -> 10
  BRN loop
  :finish
  SUB 0 13
  PRINT
  ADD 1 3 -> 3
  SET 7 <- 0
  SET 13 <- 0
  HALT
""")

# load the program
ae.load_program(program)

# indices B[k]
k = 1
# input data, initially empty, but each result is added after computation
data = []
# instruction to start execution at
start = labels['init']
# run the program
while True:
  # load the data and run the program
  ae.load_data(data)
  ae.run(start)
  # get the computed result from the output transcript
  r = (ae.output[-1] if ae.output else None)

  # display the computed result
  printf("B[{k}] = {r}")

  # run the program again?
  try:
    raw_input('\n[paused] >>> ')
  except EOFError:
    printf()
    break

  # add the result to the data and run it again
  data.append(r)
  start = labels['start']
  k += 2

[Program 7 - ada6.py]

Again, the highlighted section is the embedded program for the Analytical Engine, and the statements in the program correspond to the lines in Ada’s diagram, with extra statements for the initialisation, branching and output added. So the first row in Ada’s diagram corresponds to the “multiply” statement after the “:start” instruction (line 12 of the above Python program).

We can use the output of the assembler to provide the indices of the cards where we wish to start and subsequently re-start execution.

This, then, is a directly executable representation of the first published program.

Running the first program: Part 2

[ Part 1 | Part 2 | Part 3 ]

[Note: The programs can be run without the need to install anything @repl.it]

Yesterday was Ada Lovelace Day, so following on from my previous post about Ada’s program for computing Bernoulli Numbers, I decided to implement an emulation of the Analytical Engine so that I can directly run the “assembly language” program I produced in that post, rather than translate it into a modern programming language in order to run it.

I haven’t finished reading all the materials I have found about the Analytical Engine, but I think I’ve learned enough to have a crack at doing an emulation sufficient for my needs.

So without further ado, here’s my emulation of the Analytical Engine in Python:

# emulation of the Analytical Engine

class HaltException(Exception): pass

class AnalyticalEngine(object):

  # vars = number of variables
  # number = a function to implement the variables
  def __init__(self, **kw):
    # number of variables in the store
    self.vars = 20
    # a method to implement the variables
    self.number = float

    # set options
    for (k, v) in kw.items():
      if hasattr(self, k):
        setattr(self, k, v)
      else:
        print('AnalyticalEngine: ignoring "{k}={v}"'.format(k=k, v=v))
    
    self.reset()

  # reset the machine
  def reset(self):
    # representation of zero
    self.zero = self.number(0)
    # the variables in the store
    self.v = [self.zero] * self.vars
    # the registers in the mill
    self.index = 0
    self.input = [None, None]
    self.result = None
    # current operation
    self.op = None
    # the program and program counter
    self.program = None
    self.pc = None
    # input data
    self.data = None
    self.dc = None
    # output transcript
    self.output = None

  # load the program
  def load_program(self, program):
    self.program = program
    self.pc = 0

  # load the data
  def load_data(self, data):
    self.data = data
    self.dc = 0

  # run the program (starting at instruction start)
  def run(self, start=0):
    print(">>> Running Analytical Engine")
    self.output = []
    self.pc = start
    try:

      while True:
        # get the next instruction
        assert not(self.pc < 0), "Invalid PC" p = self.program[self.pc] # execute it getattr(self, p[0])(*p[1:]) # increment the program counter self.pc += 1 except Exception as e: print(">>> {e}".format(e=e))
      print(">>> Execution halted")


  # implement the opcodes

  # SET <i> <x>
  # set variable <i> in the store to value <x>
  def SET(self, i, x):
    self.v[i] = self.number(x)

  # ADD
  # set the engine to perform addition
  def ADD(self):
    self.op = (lambda a, b: a + b)

  # SUB
  # set the engine to perform subtraction
  def SUB(self):
    self.op = (lambda a, b: a - b)

  # MUL
  # set the engine to perform multiplication
  def MUL(self):
    self.op = (lambda a, b: a * b)

  # DIV
  # set the engine to perform division
  def DIV(self):
    self.op = (lambda a, b: a / b)

  # execute an operation
  def _execute(self):
    self.result = self.op(self.input[0], self.input[1])

  # load value v into the input register
  def _load(self, v):
    self.input[self.index] = v
    # loading the second register triggers the execution
    if self.index == 1:
      self._execute()
    # next time load the other input register
    self.index ^= 1

  # LOAD <i>
  # load the input register from variable <i> in the store
  def LOAD(self, i):
    self._load(self.v[i])

  # LOAD_DATA
  # load the input register from next value in the input data stack
  def LOAD_DATA(self):
    self._load(self.data[self.dc])
    self.dc += 1

  # STORE <i>
  # store the result to variable <i> in the store
  def STORE(self, i):
    self.v[i] = self.result

  # PRINT
  # print the result
  def PRINT(self):
    print(self.result)
    self.output.append(self.result)

  # HALT
  # halt execution of the engine
  def HALT(self):
    raise HaltException("HALT instruction encountered")

  # branch
  def _branch(self, offset):
    self.pc += offset

  # BRZ <offset>
  # branch if zero:
  # if the result is zero move the program instructions by <offset>
  def BRZ(self, offset):
    if self.result == self.zero:
      self._branch(offset)

  # BRN <offset>
  # branch if non-zero:
  # if the result is non-zero move the program instructions by <offset>
  def BRN(self, offset):
    if self.result != self.zero:
      self._branch(offset)

[Program - analytical_engine.py]

Click the link for an expanded version of this implementation that includes code to trace the execution of the Engine. [link]

Before we look at running a program on the Analytical Engine, some notes:

  1. This is an emulation of the function the Analytical Engine, not a simulation of the mechanical workings.
  2. You can specify the number of variables (columns) in the store of the Analytical Engine with the vars parameter to the AnalyticalEngine class. One of Babbage’s designs was for a store with 1000 columns.
  3. You can specify an implementation for the columns using the number parameter to the AnalyticalEngine class. Later on I’ll implement a class that more closely emulates the behaviour of the columns in the design for the Engine, but for now we can use Python’s float built-in (to use floating point numbers) or fractions.Fraction for exact results.
  4. The program is presented as a sequence of instructions (which represent the variable cards and operation cards used to program the Analytical Engine). Each instruction consists of a sequence where the first element is the “op code” and the remaining elements are the corresponding parameters. The program is loaded into the Engine using the load_program() method, and can then be executed using the run() method.
  5. The “op codes” are the names of methods on the AnalyticalEngine class. I use upper case methods to indicate the instructions.
  6. I’ve implemented a stack of “input data” cards as I suggested in my comment on the “Running the first program post” [link]. The data is loaded into the Engine using the load_data() method, and the op code used in the program cards to load a register in the mill from the input data is LOAD_DATA (rather than LOAD input as I suggested in the comment).
  7. I’ve implemented two conditional instructions: “BRZ offset” (branch if zero) and “BRN offset” (branch if non-zero). In each case the result register of the mill is checked, and if the condition holds the program counter is moved by offset instructions. The offset may be positive (to implement conditional execution) or negative (to implement looping).
  8. The Engine will continue to run the program until a (Python) exception occurs, at which point the Engine will stop. There it a HALT instruction provided specifically to raise an exception that will stop the Engine.

This is my own interpretation of the Analytical Engine, for a more rigorous interpretation and Java implementation see the pages at fourmilab.ch.

Below is a simple program to compute factorial(12). The Analytical Engine program is wrapped up in a Python program that creates a suitable instantiation of the AnalyticalEngine class (with 3 variables in the store, and using Python’s built-in int class to implement the columns).

We then load the program for the Analytical Engine, run it and the computed result is in variable 2 from the store.

from analytical_engine import AnalyticalEngine

# initialise the engine
ae = AnalyticalEngine(vars=3, number=int)

# load the program to compute factorial(12)
n = 12
ae.load_program([
  # initialisation
  ['SET', 0, n],
  ['SET', 1, 1],
  ['SET', 2, 1],
  # operation 1: v[2] = v[0] * v[2]
  ['MUL'],
  ['LOAD', 0],
  ['LOAD', 2],
  ['STORE', 2],
  # operation 2: v[0] = v[0] - 1
  ['SUB'],
  ['LOAD', 0],
  ['LOAD', 1],
  ['STORE', 0],
  # branch if non-zero to operation 1
  ['BRN', -9],
  # end
  ['HALT'],
])

# run the program
ae.run()

# the result is in v[2]
print("factorial({n}) = {r}".format(n=n, r=ae.v[2]))

[Program 3 - factorial1.py]

When run you see output like this:

$ python factorial1.py 
>>> Running Analytical Engine
>>> HALT instruction encountered
>>> Execution halted
factorial(12) = 479001600

And indeed this is the correct value for factorial(12).

The difficulty in writing these programs is making sure the branch offsets go to the right place. As statements occupy multiple cards it is easy to branch into the middle of a statement with unexpected results. (What’s really needed is an assembler that can compile a statement like “mul(0, 2, 2)” and produce the four instructions that make up the first operation. This would make a program correspond much more closely with the diagrams published to describe the operation of the Analytical Engine. Such an assembler could also allow for symbolic labels to be used in branches so the offsets are computed automatically. I might add this functionality to a future version of analytical_engine.py).

Here is the full program to execute Ada Lovelace’s algorithm for computing Bernoulli Numbers. As you can see the actual program for the Analytical Engine is a direct translation from the assembly language program I gave in my comment to the last post (except I’ve actually had to work out the correct offsets for the “branch” instructions).

In this instance we use an AnalyticalEngine with 14 variables, and the numbers are implemented using Python’s fractions.Fraction objects, which give exact results for the computation.

The controlling loop collects the results as they are generated and adds them to the input data stack, and then re-runs the program to compute the next Bernoulli Number in the sequence. It pauses between each run of the program to allow the user examine the output and to halt the program if necessary.

(Note that the enclosing Python program uses a couple of utility routines from the enigma.py library, available here [link]).

from analytical_engine import AnalyticalEngine
from fractions import Fraction
from enigma import raw_input, printf

# initialise the engine
ae = AnalyticalEngine(vars=14, number=Fraction)

# load the program
ae.load_program([
  # initialisation
  ['SET', 0, 0],
  ['SET', 1, 1],
  ['SET', 2, 2],
  ['SET', 3, 1],
  # operation 1
  ['MUL'],
  ['LOAD', 2],
  ['LOAD', 3],
  ['STORE', 4],
  ['STORE', 5],
  ['STORE', 6],
  # operation 2
  ['SUB'],
  ['LOAD', 4],
  ['LOAD', 1],
  ['STORE', 4],
  # operation 3
  ['ADD'],
  ['LOAD', 5],
  ['LOAD', 1],
  ['STORE', 5],
  # operation 4
  ['DIV'],
  ['LOAD', 4],
  ['LOAD', 5],
  ['STORE', 11],
  # operation 5
  ['DIV'],
  ['LOAD', 11],
  ['LOAD', 2],
  ['STORE', 11],
  # operation 6
  ['SUB'],
  ['LOAD', 13],
  ['LOAD', 11],
  ['STORE', 13],
  # operation 7
  ['SUB'],
  ['LOAD', 3],
  ['LOAD', 1],
  ['STORE', 10],
  # branch if zero to operation 24
  ['BRZ', +66],
  # operation 8
  ['ADD'],
  ['LOAD', 2],
  ['LOAD', 7],
  ['STORE', 7],
  # operation 9
  ['DIV'],
  ['LOAD', 6],
  ['LOAD', 7],
  ['STORE', 11],
  # operation 10
  ['MUL'],
  ['LOAD_DATA'],
  ['LOAD', 11],
  ['STORE', 12],
  # operation 11
  ['ADD'],
  ['LOAD', 12],
  ['LOAD', 13],
  ['STORE', 13],
  # operation 12
  ['SUB'],
  ['LOAD', 10],
  ['LOAD', 1],
  ['STORE', 10],
  # branch if zero to operation 24
  ['BRZ', +45],
  # operation 13
  ['SUB'],
  ['LOAD', 6],
  ['LOAD', 1],
  ['STORE', 6],
  # operation 14
  ['ADD'],
  ['LOAD', 1],
  ['LOAD', 7],
  ['STORE', 7],
  # operation 15
  ['DIV'],
  ['LOAD', 6],
  ['LOAD', 7],
  ['STORE', 8],
  # operation 16
  ['MUL'],
  ['LOAD', 8],
  ['LOAD', 11],
  ['STORE', 11],
  # operation 17
  ['SUB'],
  ['LOAD', 6],
  ['LOAD', 1],
  ['STORE', 6],
  # operation 18
  ['ADD'],
  ['LOAD', 1],
  ['LOAD', 7],
  ['STORE', 7],
  # operation 19
  ['DIV'],
  ['LOAD', 6],
  ['LOAD', 7],
  ['STORE', 9],
  # operation 20
  ['MUL'],
  ['LOAD', 9],
  ['LOAD', 11],
  ['STORE', 11],
  # operation 21
  ['MUL'],
  ['LOAD_DATA'],
  ['LOAD', 11],
  ['STORE', 12],
  # operation 22
  ['ADD'],
  ['LOAD', 12],
  ['LOAD', 13],
  ['STORE', 13],
  # operation 23
  ['SUB'],
  ['LOAD', 10],
  ['LOAD', 1],
  ['STORE', 10],
  # branch if non-zero to operation 13
  ['BRN', -45],
  # operation 24
  ['SUB'],
  ['LOAD', 0],
  ['LOAD', 13],
  # print
  ['PRINT'],
  # operation 25
  ['ADD'],
  ['LOAD', 1],
  ['LOAD', 3],
  ['STORE', 3],
  # reset working variables
  ['SET', 7, 0],
  ['SET', 13, 0],
  # end
  ['HALT'],
])

# indices B[k]
k = 1
# input data, initially empty, but each result is added after computation
data = []
# instruction to start execution at, initially 0, but subsequently 4
start = 0
# run the program
while True:
  # load the data and run the program
  ae.load_data(data)
  ae.run(start)
  # get the computed result from the output transcript
  r = (ae.output[-1] if ae.output else None)

  # display the computed result
  printf("B[{k}] = {r}")

  # run the program again?
  try:
    raw_input('\n[paused] >>> ')
  except EOFError:
    printf()
    break

  # add the result to the data and run it again (from instruction 4)
  data.append(r)
  start = 4
  k += 2

[Program 4 - ada5.py]

By running this we can check that the Engine does indeed produce the list of Bernoulli Numbers:

B[1] = 1/6
B[3] = −1/30
B[5] = 1/42
B[7] = −1/30
B[9] = 5/66
B[11] = −691/2730

The computation that Ada outlines in Note G is the calculation of B[7] = −1/30. This involves 10 ADD operations, 11 SUB operations, 8 MUL operations and 7 DIV operations. The Analytical Engine was estimated to do be able to perform an addition or subtraction in about 1 second, but multiplication and division could take up to a minute (the actual time would vary depending on the operands), and also branching would take some time to step through the required number of cards. So altogether we might expect the Engine to take around 15 to 20 minutes to compute B[7]. And as we compute further numbers in the sequence the number of operations and the corresponding time would increase.


So far we have been using Python’s built-in number classes to implement the columns, which can provide arbitrary precision. But the actual columns for the Analytical Engine would have been able to store 50 decimal digits, and a sign (to indicate if the number was positive or negative).

Here’s an implementation of a class suitable for use as the number parameter when constructing an AnalyticalEngine instance, that more closely replicates the columns of the Analytical Engine.

# implementation of the columns in the Analytical Engine

class OverflowException(Exception): pass

# a column with <digits> whole decimal digits and <dp> fractional decimal digits
def Column(digits=50, dp=0):

  shift = (10 ** dp)
  overflow = (10 ** (digits + dp)) - 1
  fmt = '<{s}{m:0' + str(digits) + 'd}' + ('.{d:0' + str(dp) + 'd}' if dp > 0 else '') + '>'

  class Column(object):

    # create a value, and check for overflow
    def __init__(self, n=0, shift=shift):
      if shift:
        n *= shift
      if abs(n) > overflow:
        raise OverflowException("Overflow in column")
      self.n = n

    # output format
    def __repr__(self):
      n = self.n
      (m, d) = divmod(abs(n), shift)
      s = ('-' if n < 0 else '+')
      return fmt.format(s=s, m=m, d=d)

    # arithmetic operations

    def __add__(self, value):
      return Column(self.n + value.n, shift=0)

    def __sub__(self, value):
      return Column(self.n - value.n, shift=0)

    def __mul__(self, value):
      return Column((self.n * value.n) // shift, shift=0)

    def __div__(self, value):
      return Column((self.n * shift) // value.n, shift=0)

    # Python 3 uses __truediv__
    __truediv__ = __div__

    # branch tests

    def __eq__(self, value):
      return self.n == value.n

    def __ne__(self, value):
      return self.n != value.n

  return Column

Notes:

  1. We can specify the number of whole decimal digits for the column with the digits parameter. The columns can also be used to implement fixed point decimal numbers by specifying how many additional fractional decimal places are required with the dp parameter.
  2. I’m not sure what happened when an overflow happened in the Analytical Engine, so I throw an exception.
  3. The columns support the four arithmetical operations that the Analytical Engine supports (addition, subtraction, multiplication and division), and also the two conditionals for the branches (is zero, is non-zero).
  4. The columns display as a fixed with string “<{sign}{digits}.{digits}>”.

The Column class is included in analytical_engine.py.

So we can use the following program to compute factorial(40) using columns containing 50 decimal digits:

from analytical_engine import AnalyticalEngine, Column

# initialise the engine
ae = AnalyticalEngine(vars=3, number=Column(digits=50))

# load the program to compute factorial(40)
n = 40
ae.load_program([
  # initialisation
  ['SET', 0, n],
  ['SET', 1, 1],
  ['SET', 2, 1],
  # operation 1: v[2] = v[0] * v[2]
  ['MUL'],
  ['LOAD', 0],
  ['LOAD', 2],
  ['STORE', 2],
  # operation 2: v[0] = v[0] - 1
  ['SUB'],
  ['LOAD', 0],
  ['LOAD', 1],
  ['STORE', 0],
  # branch if non-zero to operation 1
  ['BRN', -9],
  # end
  ['HALT'],
])

# run the program
ae.run()

# the result is in v[2]
print("factorial({n}) = {r}".format(n=n, r=ae.v[2]))

[Program 5 - factorial2.py]

$ python factorial2.py
>>> Running Analytical Engine
>>> HALT instruction encountered
>>> Execution halted
factorial(40) = <+00815915283247897734345611269596115894272000000000>

And similarly by using number=Column(10, 40) in the implementation of Ada’s program we can compute Bernoulli numbers to 40 decimal places.

In divisions we will potentially lose precision in the final decimal places of the result, and as later numbers in the sequence are calculated from earlier numbers we will see the precision deteriorate as we go on.

Here are the results of running the program:

B[1] = <+0000000000.1666666666666666666666666666666666666666>
B[3] = <-0000000000.0333333333333333333333333333333333333332>
B[5] = <+0000000000.0238095238095238095238095238095238095233>
B[7] = <-0000000000.0333333333333333333333333333333333333302>
B[9] = <+0000000000.0757575757575757575757575757575757575464>
B[11] = <-0000000000.2531135531135531135531135531135531131568>
B[13] = <+0000000001.1666666666666666666666666666666666593360>
B[15] = <-0000000007.0921568627450980392156862745098037431432>
B[17] = <+0000000054.9711779448621553884711779448621498550887>
B[19] = <-0000000529.1242424242424242424242424242422111802933>
B[21] = <+0000006192.1231884057971014492753623188306059856832>
B[23] = <-0000086580.2531135531135531135531135525557263098091>
B[25] = <+0001425517.1666666666666666666666666299288036638424>
...

By the time we reach B[25] the last 15 digits of our 40 digits of decimal precision have been lost (B[25] = 8553103/6 = 1425517 + 1/6), but the computation is still correct to 25 decimal places.

 


 

Note: The Python programs presented should run under Python 2 or Python 3, but as I used the new-style print() function, so if you are running them under Python 2 you will need to put the following line at the start of the programs:

from __future__ import print_function

Running the first program

“Notes on Note G”

[ Part 1 | Part 2 | Part 3 ]

[Note: WordPress sometimes has problems including program code in posts, so I’ve made the code used in these articles available on GitHub – [link]]

[Note: The programs can also be run without the need to install anything @repl.it]

I recently watched the BBC4 documentary Calculating Ada: The Countess of Computing (Hannah Fry’s 10-part radio series Computing Britain is also an interesting listen)[¹]. And although I know that Ada Lovelace is credited with publishing the first computer program, I hadn’t ever seen the form that this program took before. It’s not hard to track down – on Ada’s Wikipedia page there is an image of the diagram ([image] – “Diagram for the computation by the Engine of the Numbers of Bernoulli”).

The diagram itself was part of a note (Note G) that was attached to a translation of a French transcript (by an Italian military engineer named Luigi Frederico Menabrea – who would later become the Italian Prime Minister) of a seminar given by Charles Babbage at the University of Turin in 1840. It was the first and only time Babbage gave a public lecture on the Analytical Engine. Ada was commissioned to translate the paper into English. After nine months of working she published her translation, along with a set of notes (which were more than twice as long as the original paper). This was the first (and, for a hundred years, the only) text on computer programming. It is in one of the notes that the diagram that is considered to be the worlds first published computer program appeared. (The full text of Ada’s paper is available here [link], along with a large version of the diagram [link]).

Here it is:

Diagram for the computation of Bernoullinumbers

“Diagram for the computation by the Engine of the Numbers of Bernoulli”

It’s not what we think of now as a computer program – after all it pre-dates computer programming languages. Instead it essentially the description of a sequence of operations performed by the Analytical Engine. So it is more like a trace of the program execution than the source code of the program itself, and it only includes the actual arithmetic operations that the engine performs. The control flow of the execution is dealt with in the attached text.

Of course the Analytical Engine was never completed, so the actual form of the instructions is not completely clear. However, I’m not trying to construct copy of the punched cards necessary to run Ada’s program. I’m attempting to transliterate her algorithm into a modern programming language in order to understand it, and the Analytical Engine, better.

Like Ada I have found that my own notes have become somewhat longer than I anticipated.

First I’ll explain what the program does. The purpose of the program is to compute Bernoulli Numbers. Ada’s program is concerned with calculating the odd-numbered Bernoulli Numbers[²], and is based on the following equation (given in Note G) for calculating B[2n − 1], given the already computed values of B[1], B[3], …, B[2n − 3]:

A[0] + A[1]B[1] + A[3]B[3] + … + B[2n − 1] = 0

where the values of the A[] coefficients are defined as follows:

A[0] = (−1/2) (2n − 1)/(2n + 1)
A[1] = 2n/2 = n
A[3] = 2n(2n − 1)(2n − 2) / (2.3.4)
A[5] = 2n(2n − 1)(2n − 2)(2n − 3)(2n − 4) / (2.3.4.5.6)

For k > 1, each A[k] is derived from A[k − 2] by multiplying the next two decreasing terms into the numerator, and the next two increasing terms into the denominator. (Note that for each value of n the sequence of A[] coefficients are different. We might more properly think of them as A[nk]).

Ada’s program works by starting at n=1, the equation for deriving B[1] is:

A[0] + B[1] = 0

and:

A[0] = (−1/2) (2n − 1)/(2n + 1) = −1/6

So:

B[1] = 1/6

We then consider n=2, using the previously computed value for B[1] we have:

A[0] + A[1]B[1] + B[3] = 0

and:

A[0] = (−1/2) (2n − 1)/(2n + 1) = −3/10
A[1] = 2

So:

B[3] = −1/30

We then consider n=3, and use the previously computed value for B[1] and B[3] to calculate B[5]. And the procedure can be repeated as many times as we like (or until we run out of space to store the results – the proposed design for the Analytical Engine would have had room for 1,000 variables, each containing 50 decimal digits). Each time we use the results from all the previous runs to calculate the next number in the sequence.

To explain the implementation of the algorithm here is my own implementation of it using modern Python idioms. This program is constructed using the same basic structure as Ada’s program.

from fractions import Fraction

# my implementation of the algorithm
def bernoulli():

  # results
  Bs = list()

  # start at n=1
  n = 1
  # calculate the sequence
  while True:

    # result = A0
    r = -Fraction(2 * n - 1, 2 * n + 1) / 2

    # A1 = n
    A = n
    # for each B[k] already determined calculate the corresponding A[k]
    for (i, B) in enumerate(Bs):
      if i > 0:
        # multiply in the 2 additional terms
        j = 2 * i - 1
        A *= Fraction(2 * n - j, 2 + j)
        j += 1
        A *= Fraction(2 * n - j, 2 + j)
      # add A[k] * B[k] into the result
      r += A * B

    # the computed bernoulli number is -r
    B = -r
    # return the number
    yield B
    # add it to the result list
    Bs.append(B)
    # increase n
    n += 1

# run N iterations of the algorithm
N = 10

# allow N to be specified as a command line argument
import sys
if len(sys.argv) > 1:
  N = int(sys.argv[1])

k = 1 # ada's numbering (normal numbering is k + 1)
K = 2 * N + 1 # max k
for r in bernoulli():
  print("B[{k}] = {r}".format(k=k, r=r))
  k += 2
  if k == K: break

[Program 1 - ada1.py]

Before we examine a more direct implementation of Ada’s diagram, first some notes on my own implementation:

  1. (Line 1) I’m using Python’s Fraction library to generate the Bernoulli sequence. This makes it easy to check the computation against existing lists (which give them in rational form), and provides exact answers. The Analytical Engine would have computed the numbers as fixed point decimal fractions. (So small errors would accumulate in the final decimal places, and as these results are fed back into further computations the errors would compound. It would be best to use more decimal places than required in the results so that suitably accurate truncated results can be generated).
  2. (Line 4) The bernoulli() function will generate successive Bernoulli numbers. I’ve implemented it as Python co-routine, so results are returned with the yield operator as they are determined. It’s up to the calling code to consume as many results as it requires. The diagram given by Ada describes the generation of a single Bernoulli number. And while her program itself sets up some of the data required for the next run, it would also require some additional modification of the cards before it is run again to compute the next number in the sequence.
  3. (Line 18) I’ve simplified the calculation of A1 from (2n/2) to just (n).
  4. (Line 41) We can specify how many numbers in the sequence we want the program to compute (the default is 10). They are printed out according to Ada’s numbering scheme.

The program is certainly non-trivial. It involves two loops (one nested in the other) and conditional execution, as well as the storing of the computed results to use in subsequent calculations.

The diagram published by Ada shows the computation of B[7], the fourth number in the sequence. The table consists of the following columns:

1. Line number.
2. Operation.
3. Operands.
4. Assignment of results.
5-6. Comments.
7-9. Parameters.
10-19. Variables.
20-23. Results.

Only columns 2, 3 and 4 are instructions to the computer, the rest are essentially supporting documentation to explain the program to the reader.

There are a few apparent errors (the first programming bugs) in the originally published diagram:

  1. At line 4 the operation is given as v5 ÷ v4, although the comments and text make it clear that this should be v4 ÷ v5.
  2. At line 23 the formula for A3 has two 3’s in the denominator, instead of a 3 and a 4 (as in the previous line). This is just an error in the comment and would not affect execution of the program.
  3. At line 25 it looks like v6 and v7 are reset to their initial values (of zero). In fact it’s v7 and v13 that need to be reset for the algorithm to work correctly. In the accompanying text Ada refers to resetting v6, v7, and v13 to zero (although it is not necessary to zero v6, as the first thing the program does (operation 1) is to assign 2n to v4, v5, and v6).

These are all “typos” rather than “thinkos”, and may not even have been present in Ada’s original notes.

Here is my Python program with the bernoulli() function replaced with a modified version that reflects a more direct transliteration of Ada’s description:

from collections import defaultdict
from fractions import Fraction

# ada lovelace's algorithm in Python
def bernoulli():

  # initially all varibles in the store are at 0
  v = defaultdict(int)

  # program start

  # initialise the variables
  v[1] = Fraction(1) # constant
  v[2] = Fraction(2) # constant
  v[3] = Fraction(1) # n = 1

  # outer loop (compute B[2n - 1])
  while True:

    # pseudo-block to permit "break"
    while True:

      # 0: set index register
      i = 1

      # 1: v4 = v5 = v6 = 2n
      v[4] = v[5] = v[6] = v[2] * v[3]

      # 2: v4 = 2n - 1
      v[4] = v[4] - v[1]

      # 3: v5 = 2n + 1
      v[5] = v[5] + v[1]

      # 4: v11 = (2n - 1)/(2n + 1) (the diagram seems to say v[5] / v[4]) [FIX]
      v[11] = v[4] / v[5]

      # 5: v11 = (1/2) (2n - 1)/(2n + 1)
      v[11] = v[11] / v[2]

      # 6: v13 = -(1/2) (2n - 1)/(2n + 1) = A0
      v[13] = v[13] - v[11]

      # 7: v10 = n - 1
      v[10] = v[3] - v[1]
      # branch if zero to operation 24
      if v[10] == 0: break

      # 8: v7 = 2
      v[7] = v[2] + v[7]

      # 9: v11 = (2n)/2 = A1 [why not just set v[11] = v[3] instead of 8 & 9 ?]
      v[11] = v[6] / v[7]

      # 10: v12 = A1 * B1
      v[12] = v[20 + i] * v[11]
      i += 1

      # 11: v13 = A0 + A1 * B1
      v[13] = v[12] + v[13]

      # 12: v10 = n - 2
      v[10] = v[10] - v[1]
      # branch if zero to operation 24
      if v[10] == 0: break

      # for each computed result, B = B3 [1], B5 [2], ...
      while True:

        # 13: v6 = 2n - 1 [1], 2n - 3 [2], ...
        v[6] = v[6] - v[1]

        # 14: v7 = 3 [1], 5 [2], ...
        v[7] = v[1] + v[7]

        # 15: v8 = (2n - 1)/3 [1], (2n - 3)/5 [2], ...
        v[8] = v[6] / v[7]

        # 16: v11 = (2n)/2 * (2n - 1)/3 [1], (2n)/2 * (2n - 1)/3 * (2n - 3)/5 [2], ...
        v[11] = v[8] * v[11]

        # 17: v6 = 2n - 2 [1], 2n - 4 [2], ...
        v[6] = v[6] - v[1]

        # 18: v7 = 4 [1], 6 [2], ...
        v[7] = v[1] + v[7]

        # 19: v9 = (2n - 2)/4 [1], (2n - 4)/6
        v[9] = v[6] / v[7]

        # 20: v11 = (2n)/2 * (2n - 1)/3 * (2n - 2)/4 = A3 [1], (2n)/2 * (2n - 1)/3 * (2n - 2)/4 * (2n - 3)/5 * (2n - 4)/6 = A5 [2], ...
        v[11] = v[9] * v[11]

        # 21: v12 = A3 * B3 [1], A5 * B5 [2], ...
        v[12] = v[20 + i] * v[11]
        i += 1

        # 22: v13 = A0 + A1 * B1 + A3 * B3 [1], A0 + A1 * B1 + A3 * B3 + A5 * B5 [2], ...
        v[13] = v[12] + v[13]

        # 23: v10 = n - 3 [1], n - 4 [2], ...
        v[10] = v[10] - v[1]
        # branch if non-zero to operation 13
        if v[10] == 0: break

      # terminate the pseudo-block
      break

    # 24: result (-v13) is copied into the results
    r = v[20 + i] = v[20 + i] - v[13]

    # the result could be printed here
    # we pass it back to the calling routine
    yield r

    # 25: increase n, and reset working variables
    v[3] = v[1] + v[3]
    v[7] = v[13] = 0

    # branch to operation 0

# run N iterations of the algorithm
N = 10

# allow N to be specified as a command line argument
import sys
if len(sys.argv) > 1:
  N = int(sys.argv[1])

k = 1 # ada's numbering (normal numbering is k + 1)
K = 2 * N + 1 # max k
for r in bernoulli():
  print("B[{k}] = {r}".format(k=k, r=r))
  k += 2
  if k == K: break

[Program 2 - ada2.py]

Ada’s diagram shows how the program would run when n=4, and it not a direct indication of the punched cards that would be fed to Analytical Engine. Of course as the Analytical Engine was never completed it’s actual instruction set is not something we can know for certain. But the way it is described as working is as follows:

  1. An operation card indicates the operation that the engine is to perform. (One of +, −, ×, ÷).
  2. A variable card indicates which variable is transferred from the store (memory) to the first Ingress Axis of the mill (CPU). This provides the first operand.
  3. A second variable card indicates which variable is transferred from the store to the second Ingress Axis if the mill. This provides the second operand.
  4. The operation is then performed, the result appears in the Egress Axis of the mill.
  5. A variable card indicates which variable in the store the contents of the Egress Axis are transferred to.

At stage 5 there may be multiple variable cards instructing that the result is stored in several variables, as we see in the first operation of Ada’s table. If we were to write an assembly language version of this operation it would look something like this:

# operation 1
MUL
LOAD v[2]
LOAD v[3]
STORE v[4]
STORE v[5]
STORE v[6]

And each of these opcodes corresponds to a card that acts as instructions to the Analytical Engine.

Ada’s accompanying text describes the need for some common programming idioms that we take for granted in modern programming languages, but are not present in the diagram (which just shows the arithmetic operations performed).

First there is the need for conditional execution. She mentions that when computing B[1] operation 6 would complete the calculation. In my program at line 47 I exit the block when v10 becomes 0. Similarly, as Ada says, operation 12 works in a similar way. Again in my program I conditionally exit the block if v10 becomes 0.

Secondly there is the need for repeated execution of a sequence of instructions (looping). This is dealt with extensively in the notes where the repeated execution of operations 13-23 is described.

Both these could be achieved with a conditional branch instruction. Although the mechanism for branching is not explicitly mentioned in the paper it is clear that the Analytical Engine implements such a function. For example, if the Analytical Engine set a flag when a result was 0 (as microprocessors often do), then there could be an instruction to skip forward or backwards a certain number of program steps and resume operation from there. The instructions were to be provided on a punched cards, so this would require the card reader to skip forward or rewind a certain number of cards.

(In the Fourmilab simulation of the Analytical Engine a “run-up lever” is used. This is set by an operation if the sign of the result of the operation being different from the sign of the first operand (or if an overflow occurs). A conditional branch can then be made based on the state of the lever. So, we could test if a variable was zero by subtracting 1 from it. The only way the sign will change is if the variable initially held the value of 0 (so the result is −1), and the lever will be set. As Ada doesn’t show these operations in her table I’m going to simplify by referring to a “branch if zero” instruction. In a similar way branching on other conditions can be implemented).

This would allow us to place an instruction after operation 7 and operation 12 of “branch if zero to operation 24”. The inner loop would be implemented by placing an instruction after operation 23  of “branch if non-zero to operation 13”. The outer loop would be implemented by placing an instruction after operation 25 of “branch always to operation 1” (or more likely to branch if n was less than a predefined constant value, so that a Bernoulli sequence of the specified length is generated).

The text also hints at the notion of indirection, where one variable (or register) can tell us which of a number of other variables we need to access. This is necessary as each time we iterate through the operations 13-23, we have computed the A[] coefficient after executing operation 20, but then in operation 21 we multiply this the corresponding number already computed in the B[] sequence. The first time round we would need B[3], stored in variable v22 (as specified in the operands for operation 21 in the diagram), but the next time round the loop we would need the B[5] value stored in v23, and so forth. Likewise when the final result (B[7]) is calculated it is stored by operation 24 into v24 and n is increased. The next time we get to operation 24 the final result will be B[9] and it will be stored into v25. Although the diagram itself does not indicate how the instructions of operations 21 and 24 are modified as the program is executed it is clear that it is an issue that Ada has considered.

From the way Ada describes the program it seems that way this would happen is that for each value of n the program is run once. Then it must be set up slightly differently for next run. Ada notes that exactly the same sequence of operation cards can be reused, but the variable cards must be modified slightly on each iteration to achieve the required results.

One way the Analytical Engine could do indirection would be to allow the index of a variable to be specified in another variable, so instead of a “LOAD v[12]” (load the Ingress Axis from variable 12 in the store), we would have something like “LOADI v[12]” (load the Ingress Axis from the variable in the store whose index is stored in v12). So if v[12] = 21, “LOADI v[12]” would be the same as “LOAD v[21]”.

But variables in the store hold 50 decimal digits and indices into the store are only 3 decimal digits, so it is a bit odd. Instead we might consider that the Analytical Engine could be equipped with a few special purpose three decimal digit index registers, these could then be used to provide an additional offset in a load instruction.

For the purposes of Ada’s program it would be necessary to only have a single additional index register. We could imagine a that “LOADI v[20]” would correspond to “LOAD v[20 + i]” where i represents the value in the index register. So if i = 1 the instruction “LOADI v[20]” would be the same as “LOAD v[21]”. And there would be a corresponding “STOREI” instruction. To increase the offset we would need an “INC” instruction, which would increment the number in the index register by 1. This would be sufficient to implement Ada’s program.

But perhaps this approach is too influenced by modern computing. A different approach would be to allow the LOAD and STORE instructions to get a variable from an alternate stack of cards. “LOAD alt” and “STORE alt”. That way the alternate stack of cards could be filled with the addresses of v21 to (say) v30. And in operation 10 the address that the first LOAD card uses isn’t punched on the card itself, instead it is read from the alternate stack (it’s the first card in the stack so it would be v21), and the stack moves to the next card. When we reach operation 21 we would again use a “LOAD alt” instruction to get the index of the variable from the alternate stack. The first time round the loop it will be v22, then v23 and so on. When the loop ends, and we have processed all the previously calculated results the next card in the alternate stack will be the correct address to store the newly computed number. So operation 24 uses “STORE alt” and the number is stored.

Finally at the end of the diagram as well as looping back to the beginning we also rewind the alternate card deck and we are ready to go again without needing to alter any of the cards. After 10 iterations the first 10 Bernoulli numbers will have been stored in the variables referenced by the alternate stack, and the program can halt.

To produce the transliteration of Ada’s program that can run without intervention I have adopted the addition of an index register to allow indirection, although, as far as I know, this is not historically accurate.

Finally Ada notes that the process described in her diagram produces results of the correct magnitude, but not the correct sign. If there is no simple way to change the sign of a value, this can easily be remedied by subtracting the result (in v[13]) from the variable that the result is to be stored in (as it will initially be 0), and then storing it back into the result variable. Consequently I’ve adjusted operation 24 from “v[24] = v[13] + v[24]” to be “v[24] = v[24] − v[13]”, the sign is then correctly set for value to be reused by the algorithm.

Writing an equivalent pseudo-assembly language version of my transliteration would look something like this:

# initialise the variables
SET v[1] 1 # constant 1
SET v[2] 2 # constant 2
SET v[3] 1 # n = 1

# operation 0: initialise the index register
SET i 1

# operation 1
MUL
LOAD v[2]
LOAD v[3]
STORE v[4]
STORE v[5]
STORE v[6]

# operation 2
SUB
LOAD v[4]
LOAD v[1]
STORE v[4]

# operation 3
ADD
LOAD v[5]
LOAD v[1]
STORE v[5]

# operation 4
DIV
LOAD v[4]
LOAD v[5]
STORE v[11]

# operation 5
DIV
LOAD v[11]
LOAD v[2]
STORE v[11]

# operation 6
SUB
LOAD v[13]
LOAD v[11]
STORE v[13]

# operation 7
SUB
LOAD v[3]
LOAD v[1]
STORE v[10]

# branch if zero
BRANCH ZERO (operation 24)

# operation 8
ADD
LOAD v[2]
LOAD v[7]
STORE v[7]

# operation 9
DIV
LOAD v[6]
LOAD v[7]
STORE v[11]

# operation 10
MUL
LOAD v[20 + i] # indirection
LOAD v[11]
STORE v[12]
INC i # index register is incremented

# operation 11
ADD
LOAD v[12]
LOAD v[13]
STORE v[13]

# operation 12
SUB
LOAD v[10]
LOAD v[1]
STORE v[10]

# branch if zero
BRANCH ZERO (operation 24)

# operation 13
SUB
LOAD v[6]
LOAD v[1]
STORE v[6]

# operation 14
ADD
LOAD v[1]
LOAD v[7]
STORE v[7]

# operation 15
DIV
LOAD v[6]
LOAD v[7]
STORE v[8]

# operation 16
MUL
LOAD v[8]
LOAD v[11]
STORE v[11]

# operation 17
SUB
LOAD v[6]
LOAD v[1]
STORE v[6]

# operation 18
ADD
LOAD v[1]
LOAD v[7]
STORE v[7]

# operation 19
DIV
LOAD v[6]
LOAD v[7]
STORE v[9]

# operation 20
MUL
LOAD v[9]
LOAD v[11]
STORE v[11]

# operation 21
MUL
LOAD v[20 + i] # indirection
LOAD v[11]
STORE v[12]
INC i # increment index register

# operation 22
ADD
LOAD v[12]
LOAD v[13]
STORE v[13]

# operation 23
SUB
LOAD v[10]
LOAD v[1]
STORE v[10]

# branch if non-zero
BRANCH NON-ZERO (operation 13)

# operation 24
SUB
LOAD v[20 + i] # indirection
LOAD v[13]
STORE v[20 + i] # indirection

# the result can be output at this point
PRINT

# operation 25
ADD
LOAD v[1]
LOAD v[3]
STORE v[3]

# reset working variables
SET v[7] 0
SET v[13] 0

# repeat
# (or if v[3] < N)
BRANCH ALWAYS (operation 0)

HALT

Each program line would correspond to a punched card that would be fed to the Analytical Engine, so this is the closest we can get to what the original stack of punched cards, and hence the first program, would have looked like.

Note: Whilst writing up this article I have found a lot more information out about the Analytical Engine than when I started, and there is much more for me to go through. I intend to investigate the Analytical Engine further, and I will correct any glaring errors written above as I become more familiar with the subject. Please let me know of any problems you find and I will endeavour to sort them out.

— J.M.R.

Here are some useful resources that I have found about Ada, Babbage and the Analytical Engine:


[¹] Apologies if you can’t access these programmes. BBC iPlayer is only available in the UK, and programmes are only available for a limited time.

[²] Ada’s indices for the Bernoulli Numbers are one less than the indices now generally used.

%d bloggers like this: