Skip to content

Instantly share code, notes, and snippets.

@jimmckeeth
Created March 29, 2023 21:24
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save jimmckeeth/193befe1ee2c7c2998666e10c1d36810 to your computer and use it in GitHub Desktop.
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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment