L1 A:X=1 L2 Lbl 0:Y=(X+A/X)/2 L3 X≠Y⇒X=Y:Goto 0:≠⇒YΔ
The above is close to the shortest possible computer program that performs a useful numerical iteration. It prompts the user for a number A and computes its square root, correct to ten significant figures. The program was written for the Casio fx-4500P calculator, and it takes just 39 bytes of calculator memory!
An “iteration” is a way of solving a mathematical problem by producing a sequence of improving approximate answers. Every calculated answer is fed back to the same calculation as input, and repeating the process yields another, better approximation. The process is continued until subsequent answers become identical, within some prescribed error tolerance. (A mathematician would say that the procedure is “convergent”.) The drawback of iterative methods is that to begin the iteration, some initial guess as to the solution has to be provided, but in many practical cases the starting value can be quite far from the correct solution.
The problem that is being solved here is finding the square root of a positive number A. Surely, just pressing the “√” key works at least equally good on most calculators, but it doesn’t leave you with satisfaciton of having learned something new…
The algoritm is based on the observation that, if X is an estimate of √A, then the arithmetic average of X and A/X is another, better estimate. The method converges relatively fast: replacing X with (X+A/X)/2 roughly doubles the number of correct decimal places. The starting point is not so important as the method would converge for any initial value of X. In the program above, X=1 is the first approximation regardless of the value of A. (That’s right! “One” is my universal first estimate of the square root of anything!)
Surprisingly, this algorithm was already known some 2500 years ago to ancient Babylonians – they didn’t have electronic calculators but we keep finding recorded samples of calculations based on this method on old Babylonian clay tablets!
The code above is for a Casio fx-4500P programming calculator. The whole program takes just 39 bytes of memory (including a one-letter filename). Early Casio calculators are known for a highly tokenized, idiosynchratic programming language, well suited for devices that are low on memory. For example, take a look at line 3 of the program:
L3 X≠Y⇒X=Y:Goto 0:≠⇒YΔ
This is a conditional instruction, which can be translated into a simpler to understand metacode as
if (X≠Y) then X:=Y goto 0 else display Y endif
The whole line 3 takes up just 16 bytes of calculator memory (called “program steps” in Casio terminology). In particular, the “Goto” instruction would be saved internally as a one-byte token. The “⇒” symbol represents the if..then instruction, it is another one-byte token. The “Δ” token at the end of the line is an equivalent of “endif”. (On a real fx-4500P keyboard this symbol is represented by an open right-angled triangle, not to be confused with a filled triangle which temporarily stops calculation when the output is displayed.)
Plain “A” in the beginning of the first line instructs the calculator to read the value of A from user input (the program will be paused until the value is entered from the keyboard).
The iteration is continued as long as X≠Y, which means as long as the new value is different from the answer obtained in the previous step. The fx-4500P stores all numerical values with 12-digit precision but comparisons are based on 10 displayed digits, hence the calculation is stopped once all the displayed digits are correct.
The program converges really fast, the solution correct to ten significant figures is obtained in just seven iterations (taking about 2.5 seconds on my fx-4500P). To investigate the convergence of the algorithm, just add the iteration counter, as shown below. After the program has returned, check the value of “I”: it will be the number of iterations actually taken.
L1 A:X=1:I=0 L2 Lbl 0:Y=(X+A/X)/2:I=I+1 L3 X≠Y⇒X=Y:Goto 0:≠⇒YΔ