Enigmatic Code

Programming Enigma Puzzles

Running the first program: Part 2

[ Part 1 | Part 2 | Part 3 ]

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
p = AnalyticalEngine(vars=3, number=int)

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

# run the program
p.run()

# the result is in v[2]
print("factorial({n}) = {f}".format(n=n, f=p.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
p = AnalyticalEngine(vars=14, number=Fraction)

# load the program
p.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
  p.load_data(data)
  p.run(start)
  # get the computed result from the output transcript
  r = (p.output[-1] if p.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
p = AnalyticalEngine(vars=3, number=Column(digits=50))

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

# run the program
p.run()

# the result is in v[2]
print("factorial({n}) = {f}".format(n=n, f=p.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
Advertisements

Leave a Comment

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

%d bloggers like this: