## Square root without the √ key

```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Δ```

## Lorenz attractor on a fx-730P

Who said chaos can’t be modelled on an ancient machine with RAM limited to just several kilobytes? Execution time may be hours or weeks, but with a bit of careful programming, fractals and strange attractors will be calculated to the same resolution as on more capable graphic workstations.

The code at the end of this post is an attempt to solve the system of differential equations of the Lorenz attractor on a Casio fx-730P programmable calculator. The machine is low on computational power but has a fair amount of memory for programs and data (8 kb!), so I decided on a quite sophisticated numerical algorithm that offers good accuracy at a relatively low computational cost. The method is the Cash-Karp variant of the fifth-order Runge-Kutta method with adaptive stepsize control. Cash-Karp is an all-round method widely used in scientific computing. The downside is that the code is quite long; with all control and output routines it takes almost 2 kilobytes on the little Casio!

Before going into details of the code, let’s take a look at the results of the program:

The Lorenz attractor: data points were calculated on Casio fx-730P and imported (manually!) to Wolfram Mathematica for 3d plotting

The image shows the first 5000 data points of the Lorenz attractor calculated on the fx-730P. Points were joined with straight line segments and 3D-plotted in Wolfram Mathematica. The exact Mathematica command was this:

```lines = Table[{Hue[.45 + i/20000],
Line[{attr[[i]], attr[[i+1]]}]}, {i, 1, Length[attr] - 1}];
Graphics3D[{Thickness[.002], lines}, Background -> Black,
AspectRatio -> .65, Axes -> True]```

The “attr” variable is a list of 5000 points in 3D-space. This list was manually created from the output of the Casio program!

Now for some details about the Casio program that was used to calculate this. The program lets the user specify the initial value for the calculation. The initial point doesn’t really matter because the same attractor is reached from whatever point you begin, but to produce the image shown above, I started at (0,1,1).  The initial stepsize has to be provided – it again doesn’t matter too much because the algorithm will adjust the stepsize according to the prescribed accuracy. I took 0.01. The absolute error tolerance (accuracy) of the method is hard-coded in a variable called L (line number 4170). I found that the accuracy of 10^-6 is achieved at a stepsize so large that there would be surprisingly long gaps between calculated points in a 3D-plot. It shows how accurate the Cash-Karp method is, but to generate enough points for a smooth 3D-plot, an upper limit is set on the stepsize (line number 4180, variable V).

When executed, the program displays one calculated point at a time and waits until EXE is pressed. It gives you an opportunity to write down the coordinates on paper – the fx-730P doesn’t have enough memory to store all the results and the coordinates are overwritten after each step! One point takes around 12 seconds to calculate, unless the prescribed accuracy is not reached – in which case the calculation is repeated with a shorter stepsize.

The program asks in the beginning for the number of points to be calculated. The next time the program is executed, it lets you continue the calculation from the last point found previously, or start over with a new initial condition. This is how an output may look like:

The FX-730P can only display coordinates of one point at a time. To make a 3D-diagram of the attractor, coordinates have to be manually copied into graph plotting software on a more powerful machine…

The three numbers on the display are, well, x- y- and z- coordinates of a newly calculated point along the attractor. Square brackets at the beginning of the line provide some diagnostics about the stepsize:

• “[ ]” is shown when the stepsize was unchanged or increased from the previous step;
• “[M]” indicates that the stepsize was forced to the hard-coded maximum;
• “[*]” and a beep means that the stepsize was reduced by the adaptive algorithm.

Below is the complete code of the Lorenz attractor solver. When stored in calculator’s memory, the code takes up 1898 bytes (excluding REM-comments at the end). The program uses several one- and two-dimensional arrays, which when allocated, occupy another 504 bytes of RAM.

```1 REM Lorenz attractor
2 REM solved with Runge-Kutta 5th order
3 REM For Casio FX-730P
10 INPUT "[N]ew [C]ontinue",M\$
20 IF M\$="N" THEN GOSUB 4001
30 INPUT "How many steps",P
40 FOR I=0 TO 2
50 X(I)=Y(I)
60 NEXT I
70 FOR K=1 TO P
80 U\$="[ ]": IF S=V THEN U\$="[M]"
90 GOSUB 9001
100 FOR I=0 TO 2
110 G(I)=F(I)
120 NEXT I
130 GOSUB 2001
140 M=0
150 FOR I=0 TO 2
160 IF ABS(E(I))>M THEN M= ABS(E(I))
170 NEXT I
180 M=M/L
190 IM M<=1 THEN 250
200 U\$="[*]": BEEP
210 IF M>6561 THEN S=.1*S
220 IF M<=6561 THEN S=.9*S*M^-.25
230 IF T+S=T THEN PRINT "stepsize underflow"
240 GOTO 130
250 T=T+S
260 FOR I=0 TO 2
270 Y(I)=X(I)
280 NEXT I
290 IF M<=1.89E-4 THEN S=5*S
300 IF M>1.89E-4 THEN S=.9*S*M^-.2
310 IF ABS(S)> ABS(V) THEN S=V
320 PRINT U\$; RND(X(0),-4); RND(X(1),-4); RND(X(2),-4)
330 NEXT K
999 END

2001 REM RK4 algorithm routine
2002 REM Cash-Karp Runge-Kutta step
2010 FOR I=0 TO 2
2020 K(0,I)=G(I)
2030 X(I)=Y(I)+.2*S*K(0,I)
2040 NEXT I
2050 GOSUB 9001
2060 FOR I=0 TO 2
2070 K(1,I)=F(I)
2080 X(I)=Y(I)+S*(.075*K(0,I)+.225*K(1,I))
2090 NEXT I
2100 GOSUB 9001
2110 FOR I=0 TO 2
2120 K(2,I)=F(I)
2130 X(I)=Y(I)+S*(.3*K(0,I)-.9*K(1,I)+1.2*K(2,I))
2140 NEXT I
2150 GOSUB 9001
2160 FOR I=0 TO 2
2170 K(3,I)=F(I)
2180 X(I)=Y(I)+S*(C(0)*K(0,I)+C(1)*K(1,I))
2190 X(I)=X(I)+S*(C(2)*K(2,I)+C(3)*K(3,I))
2200 NEXT I
2210 GOSUB 9001
2220 FOR I=0 TO 2
2230 K(4,I)=F(I)
2240 X(I)=Y(I)+S*(C(4)*K(0,I)+C(5)*K(1,I))
2250 X(I)=X(I)+S*(C(6)*K(2,I)+C(7)*K(3,I)+C(8)*K(4,I))
2260 NEXT I
2270 GOSUB 9001
2280 FOR I=0 TO 2
2290 X(I)=Y(I)+S*(C(9)*K(0,I)+C(10)*K(2,I)+C(11)*K(3,I))
2300 X(I)=X(I)+S*C(12)*F(I)
2310 E(I)=S*(C(13)*K(0,I)+C(14)*K(2,I)+C(15)*K(3,I))
2320 E(I)=E(I)+S*(C(16)*K(4,I)+C(17)*F(I))
2330 NEXT I
2340 RETURN

4001 REM Initialize memory
4010 CLEAR
4020 DIM X(3): DIM F(3): DIM G(3)
4030 DIM Y(3): DIM E(3)
4040 DIM K(5,3): DIM C(18)
4050 C(0)=-11/54:C(1)=2.5:C(2)=-70/27:C(3)=35/27
4060 C(4)=1631/55296:C(5)=175/512:C(6)=575/13824
4070 C(7)=44275/110592:C(8)=253/4096:C(9)=37/378
4080 C(10)=250/621:C(11)=125/594:C(12)=512/1771
4090 C(13)=C(9)-2825/27648:C(14)=C(10)-18575/48384
4100 C(15)=C(11)-13525/55296:C(16)=-277/14336
4110 C(17)=C(12)-.25
4120 INPUT "Initial stepsize",S
4130 INPUT "X=",Y(0)
4140 INPUT "Y=",Y(1)
4150 INPUT "Z=",Y(2)
4160 T=0
4170 L=1E-6
4180 V=.01
4190 RETURN

9001 REM Derivatives of the Lorenz attractor ODEs
9010 F(0)=-10*(X(0)-X(1))
9020 F(1)=28*X(0)-X(1)-X(0)*X(2)
9030 F(2)=X(0)*X(1)-8*X(2)/3
9040 RETURN

10001 REM Variables
10002 REM S stepsize (changed after each step)
10003 REM P number of steps requested
10004 REM L error tolerance
10005 REM T independent variable (incremented after each step)
10006 REM M max absolute error of a step
10007 REM V largest step allowed
10008 REM Y(3) dependent variables x, y, z
10009 REM X(3) temporary storage for new values of Y()
10010 REM F(3) rhs of the attractor ODEs
10011 REM C(18) Cash-Karp coefficients
10012 REM E(3) Error estimate (calculated after each step)
10013 REM K(5,3) Runge-Kutta increments
10014 REM G(3) Temp storage for F() at the beginning of a step
10015 REM U\$ Set to "[*]" and beeps if stepsize reduced
10016 REM U\$ Set to "[M]" if working at max allowed stepsize```