Babbage

A brief explanation of the method of differences, used in Charles Babbage's Difference Engine #2

Original Mathematica code by Ed Theren, rewritten by Scott Centoni.

The Difference Engines were designs for large machines to automatically print out tables of mathematical functions by extrapolation from a small number of points, using the method of differences.

initialization

CleanSlate[]

imports

Needs["Graphics`Graphics`"]

prologs

Automatically display matrices in MatrixForm.

$PrePrint = If[MatrixQ[#], MatrixForm[#], #] &;

Improve the appearance of text in graphs.

$TextStyle = {FontFamily -> "Palatino", 
          FontSize -> 12} ; $FormatType = 
        TraditionalForm ;

misc

Print non-matrix 2D tables in a format that keeps the items centered on each line.

CheckerInterleave[x_] := 
      Most@Flatten[{#, ""} & /@ x, 1]

CheckerPad[n_][x_] := 
      Module[{
padding = 
            Table["", {k, 
                n - Length[x]}]
}, 
Join[
          padding, CheckerInterleave[x], 
          padding ]
]

CheckerForm[x_] := 
      Module[{
width = 
            Max[Length /@ x]
 (*,strlen = 
              Max@Map[Length@ToString, 
                  x, {2}]*) 
}, 
TableForm[
CheckerPad[width] /@ 
            x, 
TableSpacing -> {4, 
              1}
]
]

Pascal's triangle

Table[Binomial[n, k], {n, 0, 8}, {k, 0, n}]

CheckerForm[Table[Binomial[n, k], {n, 0, 8}, {k, 0, n}]]

method of differences

Perhaps the best way of thinking of the Difference Engine is as a machine for extrapolating tables of numbers. For instance, if we have a sequence of numbers 0,1,8,27,64, we recognize this as a list of the first few cubes. But if we did not make use of this information, how could we predict the next number? One way to approach this is to look at the successive differences from one item to the next:

Differences[x_] := Rest[x] - Most[x]

seq = {0, 1, 8, 27, 64, 125, 216};

Differences@%

Hmm, not quite obvious what the pattern is, so let's repeat the process

Differences@%

Aha...well, we can predict the next item in this sequence, but let's do it again anyway

Differences@%

Differences@%

We can see the pattern more clearly by showing all the differences together.

dat = NestList[Differences,  seq,  5]

CheckerForm[dat]

Of course, all the higher-order differences are exactly zero, so that taking the first item in each list of differences forms a new sequence, 0, 1, 6, 6, 0, ...

Now, with a bit of thought, it should be evident that this process can be reversed: if we have the sequence of coefficients, we can generate the sequence of values. That is precisely what the Difference Engine does: One inputs the difference coefficients, and cranks the machine to repeatedly add them up to obtain the values.  Matrix multiplication is a useful way to express this process (though it was not developed until after the Difference Engine was designed).

crm[x_] := crm[x] = 
        Table[If[0 ≤ j - i < 2, 1, 0], {i, x}, {j, 
            x}]

cri[x_] := cri[x] = 
        Table[If[i ≤ j, (((-1))^(i + j)), 0], {i, x}, {j, 
            x}]

This is the matrix representing cranking through a Difference Engine with 6 wheels:

crm[6]

This is the matrix that represents cranking that machine backward:

cri[6]

If we do first one, then the other, then we get the identity matrix

crm[6] . cri[6]

Crank[x_] := crm[Length[x]]  x

ReverseCrank[x_] := cri[Length[x]] . x

This basically cranks the machine once:

Crank[{0, 1, 6, 6}]

That is, 0+1=1, 1+6=7, 6+6=12.

Let's crank the machine again:

Crank@%

So 1+7=8, 7+12=19, 12+6=18.

Again:

Crank@%

Or, better yet, several steps shown in sequence:

NestList[Crank,  {0, 1, 6, 6, 0},  7]

And we see that the first column has reproduced our table of function values, plus one more--and we could keep on cranking the machine as long as we wanted. In most cases, we're only concerned with the values taken on sequentially by the register shown in the first column.

CrankTimes[regs_, n_] := First /@ NestList[Crank, regs, n]

CrankTimes[{0, 1, 6, 6, 0}, 20]

CrankThrough[regs_] := CrankTimes[regs, Length[regs] - 1]

CrankThrough[{0, 1, 6, 6, 0}]

using the machine to crank out its own coefficients

There is a trick to use the machine to generate its own coefficients from a table of regularly-spaced values. If you reverse the sign of every other item in the sequence before and after cranking through the list, it just churns out the coefficients.

Alternate[x_] := x (((-1))^Range)[0, Length[x] - 1]

DeCrankThrough[regs_] := 
      Alternate@CrankThrough@Alternate@regs

DeCrankThrough[{0, 1, 8, 27, 64}]

If we input n numbers as function values, we can output n numbers as coefficients without making any assumptions about higher-order terms.  If we output more than n numbers, we assume that the higher-order differences are all zero.

Let's look at what happens in the general case:

Array[a,  5]

DeCrankThrough[Array[a,  5]] // TableForm

This is very similar to the definition of a derivative, but note that the derivatives of different orders are being approximated at different points. We do recover the original points by cranking through the machine in the normal way:

CrankThrough@%

scaling the step size

The technique of using the Difference Engine to generate its own coefficients from a sequence of equally-spaced values is of limited practical benefit.  It's more useful if we can adjust the coefficients afterward to produce a more finely-spaced sequence.

To begin with, we consider the coefficients that correspond to the simple integer powers, 1, x, (x^2), (x^3),  ...

polycoeff[0] = {1};

polycoeff[n_] := DeCrankThrough[Table[(k^n), {k, 0, n}]]

CheckerForm[Table[polycoeff[k], {k, 0, 8}]]

We'll try the cube function again, but this time sample it with a step of two instead of one.

func = (x^3);

Table[func, {x, 0, 14}]

dat = Table[func, {x, 0, 14, 2}]

coe = DeCrankThrough[dat]

We see that instead of the coefficient list {0,1,6,6,0,...} we got before, we get larger values. But we can't just multiply the coefficients by a constant to change the step size.

It turns out that we can multiply the coefficient list by a matrix that will scale the step. To scale a 7th-order coefficient list by one half, we multiply it on the left by

halve=(({{1, 0, 0, 0, 0, 0, 0, 0}, {0, (1/2), -(1/8), (1/16), -(5/128), (7/256), -(21/1024), (33/2048)}, {0, 0, (1/4), -(1/8), (5/64), -(7/128), (21/512), -(33/1024)}, {0, 0, 0, (1/8), -(3/32), (9/128), -(7/128), (45/1024)}, {0, 0, 0, 0, (1/16), -(1/16), (7/128), -(3/64)}, {0, 0, 0, 0, 0, (1/32), -(5/128), (5/128)}, {0, 0, 0, 0, 0, 0, (1/64), -(3/128)}, {0, 0, 0, 0, 0, 0, 0, (1/128)}}));

halve  coe

This recovers the coefficient list we saw before, and the same technique works with any function.

CrankThrough[coe]

CrankThrough[halve . coe]

We can repeat this to get a step size of one half:

halve . halve . coe

CrankThrough[halve . halve . coe]

N@%

The general scaling matrix (of the appropriate order) is a function of the scaling factor r.

CoeffMat[n_] := 
      Transpose@Table[PadRight[polycoeff[k], n + 1], {k, 0, n}]

scalmat[r_, n_] := 
      CoeffMat[n] . (IdentityMatrix[n + 1] (r^Range)[0, n]) . 
        Inverse[CoeffMat[n]]

scal[r_][mat_] := scalmat[r, Length[mat] - 1] . mat

Table[MatrixForm@scalmat[r, k], {k, 0, 5}]

polynomial approximation of non-polynomial functions

Now, the real point of the Difference Engine is not to tabulate polynomial functions for their own sake, but because they can be good approximations of non-polynomial functions, like logarithmic and trigonometric functions.  There are several ways to come up with the coefficients to supply the Difference Engine, and they do not all have the same accuracy.

RMS[x_] := (Sqrt[Mean])[(x^2)]

MaxAbs[x_] := Max@Abs@x

f[x_] := Sin[x]; xmin = 0; xmax = π/8; ordermax = 7; 
    deltax = ((xmax - xmin)/ordermax); dx = ((π/180)/60);

Plot[
        f[x], {x, xmin, xmax}, 
Frame -> 
          True, 
AspectRatio -> 
          Automatic
];

coarsedat = 
      Table[{x, f[x]}, {x, xmin, xmax, ((xmax - xmin)/ordermax)}]

finedat = N[Table[{x, f[x]}, {x, xmin, xmax, dx}], 31];

finedat[[347]]

taylorleft[x_] = Normal[Series[f[x], {x, xmin, ordermax}]]

taylorleftdat = 
        N[Table[{x, taylorleft[x]}, {x, xmin, xmax, dx}], 31];

taylormid[x_] = 
      Normal[Series[f[x], {x, ((xmin + xmax)/2), ordermax}]]

taylormiddat = 
        N[Table[{x, taylormid[x]}, {x, xmin, xmax, dx}], 31];

coe = 
      DeCrankThrough[
        N[coarsedat[[All, 2]], 
          31]]

coe2 = scal[(dx/(((xmax - xmin)/ordermax)))][coe]

sd = CrankTimes[coe2, ((xmax - xmin)/dx)];

scaldat = 
        Transpose[{Table[x, {x, xmin, xmax, dx}], sd}];

lse[x_] = 
      Fit[Table[{x, f[x]}, {x, xmin, xmax, ((xmax - xmin)/ordermax)}], 
        Table[(x^k), {k, 0, ordermax}], x]

Note that Fit works only with machine precision. The more general function FindFit allows one to specify precision, though.

lse[x_] = 
      Sum[a[k] (x^k), {k, 0, ordermax}] /. 
        FindFit[Table[{x, f[x]}, {x, xmin, xmax, ((xmax - xmin)/ordermax)}], 
          Sum[a[k] (x^k), {k, 0, ordermax}], Table[a[k], {k, 0, ordermax}], x, 
          WorkingPrecision -> 31]

lsedat = N[Table[{x, lse[x]}, {x, xmin, xmax, dx}], 31];

\(\(namelist = {"taylorleft", "taylormid", "lse", "scal"};\)\)

\(TableForm[
      err = \({RMS[\((# - finedat)\)[[All, 
                  2]]], 
              MaxAbs[\((# - finedat)\)[[All, 
                  2]]]} &\) /@ 
          ToExpression[\(# <> "dat" &\) /@ 
              namelist], 
TableHeadings -> {namelist, {"RMS error", "max abs error"}}
]\)

-Log[10, err]

Our tables should be good for about this many digits:

maxdig = 
      Ceiling[-Log[10, 
            err[[3, 1]]]]

Most@#/Rest@# &@err

We see that the rms error and maximum error are quite poor for the Taylor series starting from one side, somewhat better for the Taylor series centered about the midpoint of the interval, and quite good for the methods of scaled coefficients or least squared error minimization.

ysub[a_, b_] := 
      Transpose@{a[[All, 2]], 
          a[[All, 2]] - 
            b[[All, 2]]}

\(\(\(ListPlot[
            ysub[finedat, 
              ToExpression[# <> "dat"]], 
PlotLabel -> #, 
Frame -> 
              True, 
PlotRange -> 
              All
] &\) /@ namelist;\)\)

We can easily export the table of values as a text file:

Spaces[n_] := StringJoin[Table[" ", {n}]]

\(Export["BabbageSin225.txt", 
\({#1, #2, #3, #4, #5} &\) @@@ \(({{"  d", "  m", " dd", 
                Spaces[maxdig + 4] <> "sin(x)", 
                Spaces[maxdig + 3] <> "error"}}~
            Join~\((\({ToString@
                      NumberForm[IntegerPart[#1\ 180\/π], 2, 
                        NumberPadding -> {" ", "0"}], 
                    ToString@
                      NumberForm[60  FractionalPart[#1\ 180\/π], 2, 
                        NumberPadding -> {" ", "0"}], 
ToString@
                      NumberForm[N[#1, maxdig], {maxdig, maxdig}, 
                        NumberPadding -> {" ", "0"}], 
ToString@
                      NumberForm[#2, \ {maxdig, maxdig}, 
                        NumberPadding -> {" ", "0"}], 
ToString@NumberForm[FortranForm[#2 - f[#1]], 3]} &\) @@@ 
                scaldat[[All, 
                  All]])\))\), 
"TSV"
]\)

Appendix

Taylor series

One of the useful results of calculus is that a function may be approximated to any desired accuracy near a point as a power series whose coefficients are determined by the value of the function and its derivatives at that point.  This is related to the way the Difference Engine works.

Series[g[x], {x, x0, 7}]

Series[Sin[x], {x, 0, 7}]

Plot[
        Evaluate@{Sin[x], 
            Normal@Series[Sin[x], {x, 0, 3}]}, 
{x, 
          0, π/2}, 
Frame -> 
          True
];

Trigonometric functions

Plot[{Sin[θ], Sin[θ]^2}, {θ, 0, 
          2  π}, 
Frame -> 
          True, 
FrameTicks -> {Range[0, 
              2  π, π/4], 
Range[-1, 1, 
              1/2], 
{#, (180/π) #} & /@ 
              Range[0, 
                2  π, π/
                  4], 
Automatic}, 
PlotStyle -> {GrayLevel[0], 
            Hue[0]}
];

The trigonometric functions are so related to each other that we need only square roots and the elementary arithmetic operations to convert one to another, so we will choose sin.

Since the trig functions are functions on a circular domain, we only have to compute them for 0 ≤ θ ≤ 2  π or \(TraditionalForm\`\(-π\) ≤ θ ≤ π\), etc. Since sin(π - θ) = -sin(θ), we only need to consider 0 ≤ θ ≤ π. Then, sin(π/2 - θ) = 
        sin(π/2 + θ), so we can further reduce the domain to 0 ≤ θ ≤ π/2. In fact, we can repeat this halving as many times as we wish, since \(TraditionalForm\`\(sin\^2\)(θ) = \(1 + sin(2  θ)\)\/2\), so we could fold to get 0 ≤ θ ≤ π/4, etc.

Deriving the coefficient matrix by direct brute force:

CoeffMat[omax_]:=Module[{dat, coeff},
dat = 
              Table[
Normal[Series[a[x], {x, 0, omax}]] /. 
                  x -> k dx, 
{k, 0, omax}];
coeff = DeCrankThrough[dat];
PadRight[#, omax + 1] &/@CoefficientList[coeff/.{a[0] -> 1,Derivative[n_][a][0]->(dxa^n) (dx^(-n)) n!},dxa]
]

Table[MatrixForm@CoeffMat[k],{k,8}]

Fixed-point real numbers and ten's complement notation

The Difference Engine does not use floating-point arithmetic, but rather fixed-point.  What is stored in each register is a 31-digit decimal integer. The operator must decide how to scale the coefficents so that overflow and underflow are avoided and as little precision as possible is lost.  In effect, the operator decides where to put the decimal point, but the decimal point does not not exist in the machinery.

Humans typically represent negative numbers by a negative sign followed by the absolute value, such as -736.23, but it is much simpler to build a machine to work by storing negative numbers as ten's complement (or in the case of modern binary machines, two's complement). Each digit is subtracted from one less than the base, then a one is added to the smallest digit.  The result is a number padded on the right by zeroes and on the left by zeros (for a positive number) or nines (for a negative number).  Then the sum of two numbers can be computed using the same machinery regardless of the signs of the numbers.

rep[x_] := IntegerDigits[x 1*^25, 10, 31]

neg[digs_] := 
      IntegerDigits[FromDigits[9 - digs] + 1, 10, 31]

add[a_, b_] := 
      IntegerDigits[FromDigits[a] + FromDigits[b], 10, 31]

rep[3]

neg@rep[3]

rep[5]

add[neg@rep@3, rep@5]

rep[782 + (49/100)]

neg@rep[782 + (49/100)]

rep[982 + (49/100)]

add[neg@rep[741 + (23/100)], rep[982 + (49/100)]]

Original version

endDegree = 360/16

starting degree assumed to be 0

make computing polynomial, assumed to be 8th degree x^7

Table[{i, N[i (π/56)], N@Sin[i (π/56)]}, {i, 0, 7}]

arrayValues = Table[N[Sin[i (π/56)], 31], {i, 0, 7}]

arrayPoly = Table[x^i, {i, 0, 8}]

original polynomial

arrayFitPoly = Fit[arrayValues,arrayPoly,x]

arrayStepValues = Table[ N[Sin[i/60*Pi/180],31], {i, 0, 7}]

(*  arrayStepValues = {0, 1, 8, 27, 64, 125, 216, 343 }  *)

arrayStepValues  (7/arrayStepValues)[[-1]]

For [i=1, i<9, i++,
   {  Print [i," ",N[i*cycleStep*PI/180]," ",N[arrayStepValues[[i]] ,40] ] ;
      Write[stmp, i ,"  " ,arrayStepValues[[i]] ]   
    }]

1 0.00555556 cycleStep PI 0

2 0.0111111 cycleStep PI 0.0002908882045634245963742974157400

3 0.0166667 cycleStep PI 0.0005817763845130676106143952143067

4 0.0222222 cycleStep PI 0.0008726645152351495433045892990738

5 0.0277778 cycleStep PI 0.001163552572115895060465990383306

6 0.0333333 cycleStep PI 0.001454440530541535076274490817111

7 0.0388889 cycleStep PI 0.001745328365898308835778202720850

8 0.0444444 cycleStep PI 0.002036216053572465997614191193873

(* initialize and do difference array *)
arrayDifferences = Table[N[0,31], {i, 1, 8 }]

(* initialize *)
For [ i=1, i<=8, i++,
   { 
    arrayDifferences[[i]] = arrayStepValues[[1]] ;
    For [ ii=1, ii<=8-i, ii++, 
        {arrayStepValues[[ii]]=arrayStepValues[[ii+1]]-arrayStepValues[[ii]]
        }]
    }]

?arrayDifferences

Global`arrayDifferences

({{({{\(arrayDifferences = {0, 
                      0.000290888204563424596374297415740000929832946686`30.999999999999996, \(-2.4613781582134199617173266632039246136`23.325392249211113*^-11\), \(-2.461377949941570409662642742417797322`22.84827096999262*^-11\), 
                      4.165436814809890963434227252073`15.65078449842222*^-18, 2.082717966826953606233897164052`14.951814420589713*^-18, \(-5.2869354849706366055755`7.976176747633331*^-25\), \(-1.76231113243164673568794`7.131078560626081*^-25\)}\)}})}})

let's crank the Babbage machine a few times

degrees = 0;

Print ["endDegree = ", endDegree] ;

endDegree = 22.5

For [ n=0, n<1500, n++,
   {
     minutes = Mod[{n * minutesStep}, 60] ;
     degrees = Floor[{n * minutesStep / 60} ] ;
    (* Print [      degrees, " ",minutes, " ", arrayDifferences[[1]] ] ; *)
     Write [stmp, degrees, " ",minutes, " ", arrayDifferences[[1]] ] ;
     (* do reference value *)
     reference = N[ Sin[n*minutesStep/60*Pi/180] , 31] ;
     (*Print ["   reference = ", reference ]; *)
     If [ reference != 0,
        { 
          error = reference - arrayDifferences[[1]] ;
          percentError = N[error/reference*100, 31] ; 
          (* Print [      "       ", percentError, " %" ] ; *)
           Write [stmp, "       ", percentError, " %"  ] ;
        }] ;
     For [ i=1, i<8, i++,
       {arrayDifferences[[i]]=arrayDifferences[[i]]+arrayDifferences[[i+1]] }
        ] ;
    }] ;

Close [stmp]

Print ["done"]

done

Original text

http://www.ed-thelen.org/bab/MA-S-23code.html

(* Print ["Good Morning Mr. Battage, hope you slept well."] *)

stmp = OpenWrite["MA-s-23.txt"]
Write[stmp, Good Morning Mr. OrganGrinder]
Write[stmp, Please note - this example is _not_ pipelined]

endDegree = 22.5  (* starting degree assumed to be 0  *)
(* make computing polynomial, assumed to be 8th degree x^7 *)
x=. (* clear any previous value *)
arrayValues = Table[N[Sin[i*endDegree/7*Pi/180],31], {i,0,7}]
?arrayValues
Print ["printing ?arrayValues"]
Write[stmp, "writing ?arrayValues"]
For [i=1, i<9, i++,
   {  Print [i,"  ", N[i*endDegree/7*PI/180], "  ", N[arrayValues[[i]] ,31] ] ;
      Write[stmp, i ,"  " ,arrayValues[[i]] ]   
    }]
arrayPoly = Table[x^i, {i, 0, 8}]
?arrayPoly
arrayFitPoly = Fit[arrayValues,arrayPoly,x]  (* orig polynomial *)
?arrayFitPoly

minutesStep = 1 (* minutes *)
arrayStepValues = Table[ N[Sin[i*minutesStep/60*Pi/180],31], {i, 0, 7}]
(*  arrayStepValues = {0, 1, 8, 27, 64, 125, 216, 343 }  *)
Print ["printing arrayStepValues"]
Write[stmp, "writing ?arrayStepValues"]
For [i=1, i<9, i++,
   {  Print [i," ",N[i*cycleStep*PI/180]," ",N[arrayStepValues[[i]] ,40] ] ;
      Write[stmp, i ,"  " ,arrayStepValues[[i]] ]   
    }]

(* initialize and do difference array *)
arrayDifferences = Table[N[0,31], {i, 1, 8 }]  (* initialize *)
For [ i=1, i<=8, i++,
   { 
    arrayDifferences[[i]] = arrayStepValues[[1]] ;
    For [ ii=1, ii<=8-i, ii++, 
        {arrayStepValues[[ii]]=arrayStepValues[[ii+1]]-arrayStepValues[[ii]]
        }]
    }]
?arrayDifferences

(* lets crank the Babbage machine a few times *)
stringText = "lets crank the Babbage machine a few times"
Print [stringText] ;
Write [stmp, stringText] ;
degrees = 0;
Print ["endDegree = ", endDegree] ;
For [ n=0, n<1500, n++,
   {
     minutes = Mod[{n * minutesStep}, 60] ;
     degrees = Floor[{n * minutesStep / 60} ] ;
    (* Print [      degrees, " ",minutes, " ", arrayDifferences[[1]] ] ; *)
     Write [stmp, degrees, " ",minutes, " ", arrayDifferences[[1]] ] ;
     (* do reference value *)
     reference = N[ Sin[n*minutesStep/60*Pi/180] , 31] ;
     (*Print ["   reference = ", reference ]; *)
     If [ reference != 0,
        { 
          error = reference - arrayDifferences[[1]] ;
          percentError = N[error/reference*100, 31] ; 
          (* Print [      "       ", percentError, " %" ] ; *)
           Write [stmp, "       ", percentError, " %"  ] ;
        }] ;
     For [ i=1, i<8, i++,
       {arrayDifferences[[i]]=arrayDifferences[[i]]+arrayDifferences[[i+1]] }
        ] ;
    }] ;

Close [stmp]
Print ["done"]