Lady Ada Lovelace is generally considered the world’s first computer programmer. Born Augusta Ada Byron, in 1815, she was the daughter of Lord Byron, the English romantic poet. She took an interest in mathematics and her talents lead her to work with Charles Babbage, “the father of computers,” whom she met in 1833.
Babbage started work on his Difference Engine, as a means to automatically calculate and print the values of polynomial functions. While he never saw his Difference Engine completed, two were finished in modern times based on his original designs. One of which is in the Computer History Museum in Silicon Valley, where I saw it in operation.
In 1837 Babbage first described his Analytical Engine, as a general-purpose computer. In 1843, Lovelace translated a French paper that Italian mathematician Luigi Menabrea wrote about Charles Babbage’s Analytical Engine. While translating it she added thousands of words of her own notes to the paper. As an example she wrote of one such sequence—how to calculate Bernoulli numbers. This is regarded by computer historians as the first computer program.
Two Bit History has a great write up on Lovelace’s program: What Did Ada Lovelace’s Program Actually Do? which includes a C translation and a solution for the bug it included.
Naturally I wanted to see what her code would look like in Delphi’s Object Pascal.
program LovelaceBernoulli;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
function B7: extended;
(*
* Calculates what Lady Ada Lovelace labeled "B7",
* which today we would call the 8th Bernoulli number.
* Based on https://gist.github.com/sinclairtarget/ad18ac65d277e453da5f479d6ccfc20e
* More info https://twobithistory.org/2018/08/18/ada-lovelace-note-g.html
* Bernoulli Numbers: https://en.wikipedia.org/wiki/Bernoulli_number
*)
begin
// ------------------------------------------------------------------------
// Data
// ------------------------------------------------------------------------
var v1: Single := 1; // 1
var v2: Single := 2; // 2
var v3: Single := 4; // n
// ------------------------------------------------------------------------
// Working Variables
// ------------------------------------------------------------------------
var v4: Single := 0;
var v5: Single := 0;
var v6: Single := 0; // Factors in the numerator
var v7: Single := 0; // Factors in the denominator
var v8: Single := 0;
var v10: Single := 0; // Terms remaining count, basically
var v11: Single := 0; // Accumulates v6 / v7
var v12: Single := 0; // Stores most recent calculated term
var v13: Single := 0; // Accumulates the whole result
// ------------------------------------------------------------------------
// Result Variables
// ------------------------------------------------------------------------
var v21: Single := 1.0 / 6.0; // B1
var v22: Single := -1.0 / 30.0; // B3
var v23: Single := 1.0 / 42.0; // B5
var v24: Single := 0; // B7, not yet calculated
// ------------------------------------------------------------------------
// Calculation
// ------------------------------------------------------------------------
// ------- A0 -------
(* 01 *) v6 := v2 * v3; // 2n
(* 02 *) v4 := v6 - v1; // 2n - 1
(* 03 *) v5 := v6 + v1; // 2n + 1
// In Lovelace's diagram, the below appears as v5 / v4, which is incorrect.
(* 04 *) v11 := v4 / v5; // (2n - 1) / (2n + 1)
(* 05 *) v11 := v11 / v2; // (1 / 2) * ((2n - 1) / (2n + 1))
(* 06 *) v13 := v13 - v11; // -(1 / 2) * ((2n - 1) / (2n + 1))
(* 07 *) v10 := v3 - v1; // (n - 1), set counter?
// A0 = -(1 / 2) * ((2n - 1) / (2n + 1))
// ------- B1A1 -------
(* 08 *) v7 := v2 + v7; // 2 + 0, basically a MOV instruction
(* 09 *) v11 := v6 / v7; // 2n / 2
(* 10 *) v12 := v21 * v11; // B1 * (2n / 2)
// A1 = (2n / 2)
// B1A1 = B1 * (2n / 2)
// ------- A0 + B1A1 -------
(* 11 *) v13 := v12 + v13; // A0 + B1A1
(* 12 *) v10 := v10 - v1; // (n - 2)
// On the first loop this calculates B3A3 and adds it on to v13.
// On the second loop this calculates B5A5 and adds it on.
while (v10 > 0) do
begin
// ------- B3A3, B5A5 -------
while (v6 > 2 * v3 - (2 * (v3 - v10) - 2)) do
begin // First Loop:
(* 13 *) v6 := v6 - v1; // 2n - 1
(* 14 *) v7 := v1 + v7; // 2 + 1
(* 15 *) v8 := v6 / v7; // (2n - 1) / 3
(* 16 *) v11 := v8 * v11; // (2n / 2) * ((2n - 1) / 3)
// Second Loop:
// 17 v6 := v6 - v1; 2n - 2
// 18 v7 := v1 + v7; 3 + 1
// 19 v8 := v6 / v7; (2n - 2) / 4
// 20 v11 := v8 * v11; (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4)
end;
// A better way to do this might be to use an array for all of the
// "Working Variables" and then index into it based on some calculated
// index. Lovelace might have intended v14-v20 to be used on the
// second iteration of this loop.
//
// Lovelace's program only has the version of the below line using v22
// in the multiplication.
if (v10 = 2) then
(* 21 *) v12 := v22 * v11 // B3 * A3
else
(* 21 *) v12 := v23 * v11; // B5 * A5
// B3A3 = B3 * (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4)
// ------- A0 + B1A1 + B3A3, A0 + B1A1 + B3A3 + B5A5 -------
(* 22 *) v13 := v12 + v13; // A0 + B1A1 + B3A3 (+ B5A5)
(* 23 *) v10 := v10 - v1; // (n - 3), (n - 4)
end;
(* 24 *) v24 := v13 + v24; // Store the final result in v24
(* 25 *) v3 := v1 + v3; // Move on to the next Bernoulli number!
// This outputs a positive number, but really the answer should be
// negative. There is some hand-waving in Lovelace's notes about the
// Analytical Engine sorting out the proper sign.
Result := v24;
end;
begin
try
Writeln('Calculates what Ada Lovelace labeled "B7", which today we would call the 8th Bernoulli number.');
Writeln('--------------');
Writeln('This outputs a positive number, but really the answer should be negative.');
Writeln('There is some hand-waving in Lovelace''s notes about the Analytical Engine sorting out the proper sign.');
Writeln('==============');
Writeln(Format('A0 + B1A1 + B3A3 + B5A5: %.4f',[B7]));
except
on E: Exception do
Writeln(E.ClassName, ': ', E.Message);
end;
Writeln('--------------');
Writeln('Press [enter] to close');
readln;
end.
While this is one way to calculate the 8th Bernoulli number, what would a modern implementation look like? Using Rudy’s Big Numbers library I created a sample application to calculate any Bernoulli number.
function Bernoulli(n: Uint64): BigRational;
begin
var a: TArray<BigRational>;
SetLength(a, n + 1);
for var m := 0 to High(a) do
begin
a[m] := BigRational.Create(1, m + 1);
for var j := m downto 1 do
begin
a[j - 1] := (a[j - 1] - a[j]) * j;
end;
end;
Result := a[0];
end;
You can install the BigNumbers library via GetIt and this sample application is included.