100 REM ITERATED FUNCTION SYSTEMS Ver 1.0 (c) 1993 by J. C. Sprott
110 DEFDBL A-Z: DIM A(12)
120 RANDOMIZE TIMER             'Reseed random number generator
130 SCREEN 12                   'Assume VGA graphics
140 GOSUB 300                   'Set parameters
150 GOSUB 400                   'Iterate equations
160 GOSUB 500                   'Display results
170 GOSUB 600                   'Test results
180 ON T% GOTO 130, 140, 150
190 END

300 REM Set parameters
310 X = .05: Y = .05: XE = X + .00001: YE = Y
320 FOR I% = 1 TO 12: A(I%) = .1 * (INT(25 * RND) - 12): NEXT I%
330     J0 = ABS(A(1) * A(4) - A(2) * A(3))
340     J1 = ABS(A(7) * A(10) - A(8) * A(9))
350 IF J0 + J1 = 0 OR J0 > 1 OR J1 > 1 THEN GOTO 320    'Not contracting
360 P = J0 / (J0 + J1)
370 T% = 3: LSUM = 0: N = 1: N1 = 0: N2 = 0
380 XMIN = 1000000#: XMAX = -XMIN: YMIN = XMIN: YMAX = XMAX
390 RETURN

400 REM Iterate equations
410 IF X <> XE THEN IF RND > P THEN R% = 6 ELSE R% = 0
420 XNEW = A(1 + R%) * X + A(2 + R%) * Y + A(5 + R%)
430 YNEW = A(3 + R%) * X + A(4 + R%) * Y + A(6 + R%)
490 RETURN

500 REM Display results
510 IF N < 100 OR N > 1000 THEN GOTO 560
520     IF X < XMIN THEN XMIN = X
530     IF X > XMAX THEN XMAX = X
540     IF Y < YMIN THEN YMIN = Y
550     IF Y > YMAX THEN YMAX = Y
560 IF N = 1000 THEN GOSUB 800
570 IF X > XL AND X < XH AND Y > YL AND Y < YH AND N > 1000 THEN PSET (X, Y)
590 RETURN

600 REM Test results
610 GOSUB 700                   'Calculate Lyapunov exponent
620 GOSUB 900                   'Calculate correlation dimension
630 IF N > 21000 THEN T% = 2    'Candidate IFS found
640 IF ABS(XNEW) + ABS(YNEW) > 1000000# THEN T% = 2     'Unbounded
650 IF N > 998 AND (F < 1 OR L > -.2) THEN T% = 2       'Uninteresting
660 IF LEN(INKEY$) THEN T% = 0  'User key press
670 X = XNEW: Y = YNEW: N = N + 1
690 RETURN

700 REM Calculate Lyapunov exponent
710 XSAVE = XNEW: YSAVE = YNEW: X = XE: Y = YE
720 GOSUB 400                   'Reiterate equations
730 DLX = XNEW - XSAVE: DLY = YNEW - YSAVE: DL2 = DLX * DLX + DLY * DLY
740 DF = 10000000000# * DL2: RS = 1# / SQR(DF)
750 XE = XSAVE + RS * (XNEW - XSAVE): YE = YSAVE + RS * (YNEW - YSAVE)
760 XNEW = XSAVE: YNEW = YSAVE
770 LSUM = LSUM + LOG(DF): L = .721347 * LSUM / N
790 RETURN

800 REM Resize the screen (and discard the first thousand iterates)
810 DX = .1 * (XMAX - XMIN): DY = .1 * (YMAX - YMIN)
820 XL = XMIN - DX: XH = XMAX + DX: YL = YMIN - DY: YH = YMAX + DY
830 IF XH - XL < .000001 OR YH - YL < .000001 THEN GOTO 890
840     WINDOW (XL, YL)-(XH, YH): CLS
850     LINE (XL, YL)-(XH, YH), , B
890 RETURN

900 REM Calculate fractal dimension
910 IF N < 200 THEN GOTO 990                     'Wait for transient to settle
920 IF N = 200 THEN D2MAX = (XMAX - XMIN) ^ 2 + (YMAX - YMIN) ^ 2
930 IF N = 200 OR RND < .02 THEN XS = X: YS = Y  'New reference point
940 DX = XNEW - XS: DY = YNEW - YS
950 D2 = DX * DX + DY * DY
960 IF D2 < .001 * D2MAX THEN N2 = N2 + 1
970 IF D2 < .00001 * D2MAX THEN N1 = N1 + 1: F = .434294 * LOG(N2 / N1)
990 RETURN