Category: Delphi Projects

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.

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;