The Sierpinski curve: a lesson in debugging and conversion. David H. Ahl.

In several recent tutorials, I promised to publish a routine for plotting the Sierpinski curve, one of the most fascinating curves in mathematics (see box). I had seen several ways of approaching the problem, but by far the most elegant was a demonstration program on a NEC 8801 computer which used the N. Wirth algorithm.

I managed to get a listing of the relevant portion of the program, but upon examining it more closely, I realized that it employed several statements with which I was unfamiliar. Several were minor variations of statements in MSX Basic, but LINE - STEP (H,H) and its variants really had me stumped (see lines 7370-7400 in Listing 1). To make matters worse, I had no way of easily obtaining the Basic manual for the system. Moreover, the five books I had on MBasic had no description of the graphics commands.

Sometime in the 1700's, mathematicians described the word "curve" more precisely than had been done previously by defining a curve as the loci of points that satisfy equations that are continuous functions.

If a curve describes a continuous function, it certainly ought to

be possible to draw a tangent to any point on the curve.   However,

by the mid-1800's, mathematicians began to find strange, new curves that had no unique tangent at any point. One of the most famous such curves was one described in 1890 by Giuseppe Peano. This Italian logician showed how a single point, moving continuously over a square, could pass at least once through every point on the surface of the square and its boundary. Peano's curve is a legitimate diagram of a continuous function, yet nowhere on it can a unique tangent be drawn because at no instant can one specify the direction in which the point is moving.

It wasn't long before other logicians proposed curves with similar properties. Two of the most interesting are the curves devised by David Hilbert and Waclaw Sierpinski. Figure 1 shows the method of constructing the basic closed Sierpinski curve. Figure 2 shows the construction of a second order curve.

An interesting problem for computer enthusiasts is to find the area bounded by the Sierpinski curve at its upper limit. From Figure 1, we can see that of the 16 smaller squares, a first order Sierpinski curve occupies four squares completely plus 12 one-eighths, or a total of 11/32 of the total.

Looking at Figure 2, we see that one-quarter of the second order curve occupies the same area as the first order curve plus 7/8 of the corner square, or 11/32 plus 7/8 /16 which equals 51/128.

If we generalize this progression, we find that the denominator of each new term is four times the previous numerator plus 7. Thus the first few terms of the series are: 11/32 51/128 211/512 851/2048 3411/8192 13651/32768

I leave it up to readers to write a short program to determine the area bounded by the upper limit Sierpinski curve. This area should be expressed as a fraction.

There is much more to be said about Sierpinski and other "monster" curves, but Martin Gardner has said it much more elegantly than I can. See especially "Mathematical Games" in Scientific American, Sept. 1976.

This is a situation frequently faced by programmers. You would like to take advantage of the work done by someone else and not, so to speak, reinvent the wheel. However, the routine you want is in another language or a different dialect of the language in which you are working. What to do?

Since I had a good idea of the shape of the Sierpinski curve, I could deduce what program should be doing. Basically, the end point of each new line should be fairly close to the previous point and should differ by a one or two units in the x or y direction, or both. Consequently, it appeared that the two variables in the STEP portion of the LINE statement were x and y increments. Since MSX Basic does not have an incremental plot function, I substituted a subroutine (lines 800-820) and two variables, A and B, for the LINE STEP statement. Cumbersome, but I thought it would work.

Unfortunately, the program in Listing 2 did nothing of the kind, and produced the curve in Figure 1. Obviously something was wrong. But what?

My first thought was that the plotting subroutine perhaps did not need the dummy variables C and D to represent X and Y, so I eliminated them from the subroutine. I was rewarded with the curve in Figure 2--hardly an improvement. 800 ' Plotting Subroutine 810 LINE(X,Y)-(X+A,Y+B) 820 X=X+A:Y=Y+B:RETURN

Well, of course! The initial value of HO (see line 100) should be set to the vertical height of the screen, or 192. At least that's what I thought before I saw the output in Figure 3. Actually, it was slightly better than the previous one, but it occurred to me that I really should start at the beginning--with just one, large curve. 100 H0=192;SP=0:H=H0/4:X=2*H:Y=3*H:I=0

Hence, i changed the main loop to go from 1 to 1 and wound up with the curve in Figure 4. Ah, that's better; at least the program can draw one complete curve successfully. 20 FOR DI=1 TO 1 30 GOSUB 100 40 NEXT DI

Next I enlarged the loop to draw the second curve and got the curve shown in Figure 5. 20 FOR DI=1 TO 2 30 GOSUB 100 40 NEXT DI

Where were the last three segments of the second order curve? I had no idea, so I pored over the original program once again and compared it to my translated version. Lo and behold, I found a line missing in subroutine D, in particular, the one in which it calls itself. I inserted this as line 560, ran the program again, and was rewarded with two complete curves (Figure 6). 500 ' Subroutine D 510 IF TP<=0 THEN RETURN 520 PS=TP-1:GOSUB 600 530 GOSUB 500:A=H:B=H:GOSUB 800 540 GOSUB 200:A=0:B=2*H:GOSUB 800 550 GOSUB 400:A=-H:B=H:GOSUB 800 560 GOSUB 500 570 GOSUB 700 580 RETURN

You might wonder why I didn't just start here and published this debugged program in Creative Computing. After all, you certainly don't care about the incorrect version. My reason for publishing the prior steps is to illustrate the process of conversion from one dialect of Basic to another and to show that, all too often, problems that we blame on the program are something else entirely. You would be amazed at the number of letters and phone calls I get about a program in Basic Computer Games or Creative Computing that "doesn't work" because the person copying it has left out a line or copied incorrectly.

Back to Mr. Sierpinski. Things now appeared to be in order, so I expanded the number of curves to three and got a correct plot. 20 FOR DI=1 TO 3 30 GOSUB 100 40 NEXT DI

What happens when the upper and lower bounds of the main loop are changed? I tried the program with limits of 2 and 4 and got the delightful output shown in Figure 8. 20 FOR DI=2 TO 4 30 GOSUB 100 40 NEXT DI

Could a fifth level curve be plotted? Yes, but the resolution of my Spectra Video 328 and Sakata monitor started to go over the brink (see Figure 9), so I settled for a fourth level curve.

Next, I introduced some color with two additional statements (25 and 35). This colors the background after the first curve is drawn and then uses different colors for the next two curves. The effect is quite striking! 20 FOR DI=2 TO 4 25 CO=4*(DI-1) 30 GOSUB 100 35 IF DI=2 THEN PAINT (191,191),CO 40 NEXT DI

Okay, everything now worked and it was time to tidy up the program with a few remarks, labels, and blank lines. The final program is shown in Listing 3. It is a great deal of fun to watch these curves being traced out on the screen; try it!

Putting the screen plot routine into a subroutine instead of 16 separate LINE statements as in the original had an unexpected benefit, namely that it can be easily changed for different computers or for a plotter. For example, suppose you have a Houston Instruments DMP-29 plotter. To plot the Sierpinski curve on your plotter requires just the following five new lines. 30 PRINT #1,"Z ;: A" 70 PRINT #1,"U" 90 PRINT #1,"P0":END 125 PRINT #1,X,Y "D" 810 PRINT #1,X+A,Y+B

Line 30 initializes the plotter; line 125 moves the pen to the starting point and puts it down; line 810 does the actual plot; line 70 raises the pen; and line 90 moves the pen home. If you have a different plotter, you should have no trouble getting this routine to work. The curve at the opening of the article was done with this program.

This program, incidentally, takes advantage of the feature in MSX Basic and MBasic that permits a subroutine to call itself. This is known as recursion. It is often said that languages such as Basic and Fortran do not permit recursion. This is simply not true. While not all versions of Basic permit a subroutine to call itself, there are other ways of achieving recursion, but that is a subject for another day. (If you wish to delve into it, pick up a copy of any one of my Ideabooks, and you will find an entire chapter devoted to solving problems using recursion in Basic).

Be sure to look at the following article in which William Fujimoto has devised several interesting variations of the Sierpinski curve using the C language on an Altos 8600 with a Houston Instruments DMP-2 plotter. (No, I don't pretend you can easily translate from C to Basic, but using this Basic program you ought to be able to implement some of the variations suggested by Mr. Fujimoto).

For example, a tilted curve can be implemented in the Basic program just by adding and changing a few lines. Dummy plot variables (XP,YP,XQ, and YQ) are set up in lines 115 and 116 and then updated in the plot subroutine. This produces the "tilted" plots in Figures 10 and 11. 115 YP=SQR(Y)*7 116 XP=X*(-Y/(2*H0)+1)+Y/4 800 ' Plotting Subroutine 805 X=X+A:Y=Y+B 810 YQ=SQR(Y)*7 812 XQ=X*(-Y/(2*H0)+1)+Y/4 815 LINE (XP,YP)-(XQ,YQ) 820 XP=XQ:YP=YQ:RETURN

If you would rather have straight sides as in Figure 12, it is simply a matter of changing two lines. 116 XP=X*(-YP/H0+1)+YP/2 812 XQ=X*(-YP/H0+1)+YP/2

But enough. I could go on forever with more variations. But unlike the Sierpinski curve which has no end, this article does. Here it is. Listing 1. Statements from original NEC program. 7300 '-- display routine -- 7310 *SIERPIN

7320 HO=512:SP=O:H=HO 4:X=2*H:Y=3*H:I=O
7330 I=I+1:X=X-H:H=H 2:Y=Y+H

7340 IF I<DI THEN 7330 7350 POINT(X,Y):COLOR ,,,COL 7360 PSVAL=I:GOSUB *PUSH 7370 GOSUB *SUBAA:LINE -STEP(H,-H) 7380 GOSUB *SUBBB:LINE -STEP(-H,-H) 7390 GOSUB *SUBCC:LINE -STEP(-H,H) 7400 GOSUB *SUBDD:LINE -STEP(H,H) 7410 GOSUB *POP 7420 RETURN 7430 7440 *SUBAA