-- An implementation of the algorithm derived by Ada Lovelace to calculate -- the Bernoulli Numbers, as detailed in her "Note G", an addendum to her -- translation of "Sketch of The Analytical Engine" by Luigi Federico -- Menabrea. -- -- You may find a scan of the original table detailing the algorithm here: -- -- There is a more readable transcription of Note G (as well as the original -- translated text and her other Notes) available here: -- -- -- It should be noted that she intended for this to be a memoized, iterative -- process; i.e, the Engine would calculate Bernoulli number 2, then store -- that on a "variable card" that would be fed back in to calculate Bernoulli -- number 4, and so on. So in this algorithm the result variables were -- initialized to the known values of Bernoulli numbers 2, 4, and 6; she -- emphasizes that these could've been easily calculated in previous -- iterations by the Engine. However, the version detailed in the Note and -- implemented here is specific towards calculating Bernoulli number 8, and -- does not have the proper conditional or looping constructs included in -- order to make it fully generic towards calculating any Bernoulli number. -- -- Compile this with -- gnatmake note_g.adb -- and then run -- ./note_g -- -- -- Originally written in 1843 by Ada Augusta, Countess of Lovelace -- -- The original text, diagrams, and algorithms of Note G are free of known -- copyright restrictions. -- -- This version written in 2022 by nytpu -- -- To the extent possible under law, the author(s) have dedicated all -- copyright and related and neighboring rights to this software to the public -- domain worldwide. This software is distributed without any warranty. -- -- You should have received a copy of the CC0 Public Domain Dedication along -- with this software. If not, see . with Ada.Decimal; with Ada.Text_IO; use Ada.Text_IO; procedure Note_G is ----------- -- Types -- ----------- -- The Analytical Engine operated on 40-decimal-digit numbers, however a -- significand of 40 is rather large so we instead use the maximum -- available by the system if 40 digits are not available. On my x86-64 -- GNAT installation, 38 digits are available which is very close to the -- true number anyways. Num_Digits : constant := Positive'Min(Ada.Decimal.Max_Decimal_Digits, 40); -- In Note G Lovelace states: -- -- "In any of the examples already given in the translation and in the -- Notes, some of the data, or of the temporary or permanent results, -- might be fractional, quite as probably as whole numbers. But the -- arrangements are so made, that the nature of the processes would be -- the same as for whole numbers." -- -- Basically, she's saying that the calculation would use (decimal-based) -- fixed-point numbers, with the point placed so the result will be able -- to accomodate both the integral and fractional parts. Unfortunately -- where the point is supposed to be placed is not stated. Since all the -- Bernoulli numbers up to 10 had already been calculated at that time, it -- was known that all the numbers up to "B7" had an integral part of 0 or -- ±1---it happens to be that all Bernoulli numbers up to 18 have an -- absolute value of less than 10. So it would be safe to assume we could -- place the point so there's one integral digit and the remainder are -- fractional, right? -- -- However, it turns out that there are a few intermediary calculations -- that are > 10. Hence, we place the point so there is two integral -- decimal digits and the remainder are fractional, for maximum accuracy. type AE_Digit is delta 10.0 ** (-(Num_Digits - 2)) digits Num_Digits; -- Uninitialized variables would be set to zero. Uninitialized : constant AE_Digit := 0.0; ---------- -- Data -- ---------- -- Helper constants, since there was no built-in "subtract by 1" or -- "divide by 2" operations in the Analytical Engine nor the ability to -- provide constant literals as parameters to operations. V1 : constant AE_Digit := 1.0; V2 : constant AE_Digit := 2.0; -- Stated in the diagram as "n" -- The algorithm will calculate the Bernoulli number 2n---stated by -- Lovelace as calculating number "2n - 1", since she considered 1/6 -- rather than 1/2 to be the first Bernoulli number (alternately she -- could've also been the inventor of zero-indexing). -- -- Since the algorithm is generic enough to calcualate any Bernoulli -- number, I personally believe (with no corresponding evidence) that -- she picked n=4 because the resulting number of steps and required -- variables fit nicely on the page. V3 : AE_Digit := 4.0; ----------------------- -- Working Variables -- ----------------------- V4, V5, V6, V7, V8, V9, V10, V11, V12, V13 : AE_Digit := Uninitialized; -- We roll operations 13-16 and 17-20 into a single loop, so V9 is not -- used in the second iteration as it is used in the origninal table. pragma Unreferenced(V9); -- She does not explicitly mention these variables, but implies their -- existence by saying "V13..." and placing the result variables at V21 -- onwards. It can be guessed that they should be used for the second -- iteration of steps 13-23, however it would be an unecessary -- complication. -- V14, V15, V16, V17, V18, V19, V20 : AE_Digit := Uninitialized; ---------------------- -- Result Variables -- ---------------------- -- Precalculated B1 (Bernoulli number 2) V21 : constant AE_Digit := 1.0 / 6.0; -- Precalculated B3 (Bernoulli number 4) V22 : constant AE_Digit := -(1.0 / 30.0); -- Precalculated B5 (Bernoulli number 6) V23 : constant AE_Digit := 1.0 / 42.0; -- Will hold the final B7 (Bernoulli number 8) V24 : AE_Digit := Uninitialized; begin -------------------- -- Section 1 -- -- Operations 1-7 -- -- Calculate A0 -- -------------------- -- Operation 1 -- 2n V4 := V2 * V3; V5 := V2 * V3; V6 := V2 * V3; -- Operation 2 -- 2n - 1 V4 := V4 - V1; -- Operation 3 -- 2n + 1 V5 := V5 + V1; -- Operation 4 -- The first programming bug? -- Both within the explanatory notes and in the diagram she lists that the -- results stored in V11 should be "(2n - 1) / (2n + 1)". However, either -- due to a typesetting error or due to a bug introduced on her part, the -- operation is incorrectly listed as "V5 / V4", or "(2n + 1) / (2n - 1)"! -- V11 := V5 / V4; -- The corrected Operation 4: -- (2n - 1) / (2n + 1) V11 := V4 / V5; -- Operation 5 -- (1/2) * ((2n - 1) / (2n + 1)) V11 := V11 / V2; -- Operation 6 -- -(1/2) * ((2n - 1) / (2n + 1)) V13 := V13 - V11; -- A0 := -(1/2) * ((2n - 1) / (2n + 1)) -- Operation 7 -- n - 1 -- I believe this is to decrement the "counter" of the remaining logical -- steps to execute. This counter was likely intended to be used in the -- looping iterative version of the algorithm. V10 := V3 - V1; -- If V10 = 0 (i.e. V3 = 1, or B1 had not been calculated yet), we would -- store V13 to V21, then add one to V3 and restart from the beginning. ------------------------- -- Section 2 -- -- Operations 8-12 -- -- Calculate A0 + B1A1 -- ------------------------- -- Operation 8 -- 2 + 0 -- (basically a hack to emulate a load/move instruction) V7 := V2 + V7; -- Operation 9 -- 2n / 2 V11 := V6 / V7; -- A1 := (2n / 2); -- Operation 10 -- B1 * A1 V12 := V21 * V11; -- B1A1 := B1 * (2n / 2); -- Operation 11 -- A0 + B1A1 V13 := V12 + V13; -- Operation 12 -- n - 2 -- Yet another decrement of the "step counter" V10 := V10 - V1; -- If V10 = 0 (i.e. V3 = 2, or B3 had not been calculated yet), we would -- store V13 to V22, then add one to V3 and restart from the beginning. ---------------------- -- Section 3 -- -- Operations 13-25 -- -- Calculate B7 -- ---------------------- -- Operations 13-23, repeated twice -- On the first iteration it will calculate B3A3 and add it to V13 -- On the second iteration it will calculate B5A5 and add it to V13 -- -- If V10 = 1 (i.e. V3 = 3, or B5 had not been calculated yet), then -- this loop would only execute once and we would instead store the result -- in V23. However, Operations 13-16 would need to be modified to execute -- the correct number of times based off of the current value of V10. while V10 > 0.0 loop -- Operations 13-16, then 17-20 -- 13: 2n - 1 -- 14: 2 + 1 -- 15: (2n - 1) / 3 -- 16: (2n / 2) * ((2n - 1) / 3) -- 17: 2n - 2 -- 18: 3 + 1 -- 19: (2n - 2) / 4 -- 20: (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4) for I in 1 .. 2 loop V6 := V6 - V1; V7 := V1 + V7; V8 := V6 / V7; V11 := V8 * V11; end loop; -- A3 := (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4); -- Operation 21 if V10 = 2.0 then -- B3 * A3 V12 := V22 * V11; -- B3A3 := B3 * (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4); else -- B5 * A5 V12 := V23 * V11; -- B5A5 := B5 * (2n / 2) * ((2n - 1) / 3) * ((2n - 2) / 4); end if; -- Operation 22 -- V10 = 2: A0 + B1A1 + B3A3 -- V10 = 1: A0 + B1A1 + B3A3 + B5A5 V13 := V12 + V13; -- Operation 23 -- V10 = 2: n - 3 -- V10 = 1: n - 4 -- Decrement the step counter again V10 := V10 - V1; end loop; -- Operation 24 -- Without this negation, then B7 will be incorrectly output as a negative -- number. Lovelace says: -- "In the above table and diagram we are not considering the signs of any -- of the B's, merely their numerical magnitude. The engine would bring -- out the sign for each of them correctly of course, but we cannot enter -- on every additional detail of this kind as we might wish to do." So she -- does some hand-waving to dismiss sign issues, however as far as I can -- tell, negating it here will result in the correct answer. -- V24 := V13 + V24; -- The corrected Operation 24: -- 0 - (A0 + B1A1 + B3A3 + B5A5) V24 := V24 - V13; -- Operation 25 -- 4 + 1 -- Add one to the counter, so the Engine could continue calculating the -- next Bernoulli number! V3 := V1 + V3; -- Output the final result. Put_Line("B7 = A0 + B1A1 + B3A3 + B5A5 = " & V24'Image); end Note_G;