[ 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:

- This is an emulation of the function the Analytical Engine, not a simulation of the mechanical workings.
- 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.
- 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.
- 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.
- The “op codes” are the names of methods on the
`AnalyticalEngine` class. I use upper case methods to indicate the instructions.
- 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).
- 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).
- 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:

- 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.
- I’m not sure what happened when an overflow happened in the Analytical Engine, so I throw an exception.
- 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).
- 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

### Like this:

Like Loading...

## Recent Comments