-
-
Save jimmckeeth/193befe1ee2c7c2998666e10c1d36810 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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