How many CPU hours did I burn in the early 1990s rendering bits of the Mandelbrot Set? A lot, mainly because I was doing it in BASIC on an unexpanded 8 MHz Commodore Amiga A500. The image below that Fraqtive rendered in almost no time would have taken me days:
But it turns out that the first rendering of what we now call the Mandelbrot set wasn’t produced by Benoit Mandelbrot, but by Brooks & Matelski a year or two earlier:
What I’ve done – and mostly thanks to tweaks by commenter nobody below – is create period-appropriate code to reproduce that graphic. Since the paper was presented in 1978, there’s a fair chance that the authors had access to a machine running FORTRAN-77 or a near equivalent. FORTRAN’s particularly good for this:
- it has a built-in COMPLEX type that extends regular mathematical functions;
- it has just good enough string handling to output a line of spaces/asterisks. I would not have wanted to write this in FORTRAN-66, as that language had no string manipulation facilities at all.
So here’s the code. It should compile on any Fortran compiler:
PROGRAM BRKMTF
! GENERATE FIGURE FROM BROOKS-MATELSKI PAPER C.1978
! THAT EVENTUALLY BECAME KNOWN AS THE MANDELBROT SET
! - SCRUSS, 2022-05
! REF: BROOKS, ROBERT, AND J. PETER MATELSKI.
! "THE DYNAMICS OF 2-GENERATOR SUBGROUPS OF PSL (2, C)."
! RIEMANN SURFACES AND RELATED TOPICS: PROCEEDINGS OF THE
! 1978 STONY BROOK CONFERENCE,
! ANN. OF MATH. STUD. VOL. 97. 1981: FIG. 2, P. 81
REAL MAP, CR, CI
INTEGER I, J, K, M, ROWS, COLS, MAXIT
COMPLEX C, Z
PARAMETER (ROWS=31, COLS=71, MAXIT=200)
CHARACTER*80 OUT
CHARACTER CH*1
DO J=1,ROWS
CI=MAP(REAL(J), 1.0, REAL(ROWS), -0.8715, 0.8715)
DO I=1,COLS
CR=MAP(REAL(I), 1.0, REAL(COLS), -1.975, 0.475)
C=CMPLX(CR, CI)
Z=CMPLX(0.0, 0.0)
CH='*'
DO 100, K=1,MAXIT
Z = Z**2 + C
IF (ABS(Z) .GT. 2) THEN
CH=' '
GO TO 101
END IF
100 CONTINUE
101 OUT(I:I)=CH
END DO
WRITE(*,*)OUT
END DO
END
REAL FUNCTION MAP(X, XMIN, XMAX, YMIN, YMAX)
REAL X, XMIN, XMAX, YMIN, YMAX
MAP = YMIN + (YMAX - YMIN) * ((X - XMIN) / (XMAX - XMIN))
END
The results are spot-on:
Maybe Brooks & Matelski had access to an Apple II and wrote something in BASIC? I could be entirely period-accurate and write something in PDP-8 BASIC on my SBC6120, but not today.
It really is much easier using a language with complex number support when working with the Mandelbrot set. Here’s the same program in Python3, which bears more of a resemblance to FORTRAN-77 than it might admit:
#!/usr/bin/python3
# brkmtf - Brooks-Matelski proto ASCII Mandelbrot set - scruss, 2022-05
# -*- coding: utf-8 -*-
def valmap(value, istart, istop, ostart, ostop):
return ostart + (ostop - ostart) * ((value - istart) / (istop - istart))
rows = 31
cols = 71
maxit = 200
for y in range(rows):
ci = valmap(float(y + 1), 1.0, float(rows), -0.8715, 0.8715)
for x in range(cols):
cr = valmap(float(x + 1), 1.0, float(cols), -1.975, 0.475)
c = complex(cr, ci)
z = complex(0.0, 0.0)
ch = '*'
for k in range(maxit):
z = z**2 + c
if (abs(z) > 2.0):
ch = ' '
break
print(ch, end='')
print()
I found out about the Brooks-Matelski paper from this article: Who Discovered the Mandelbrot Set? – Scientific American. It’s none too complimentary towards Benoit Mandelbrot.
I did this on a Dragon 32 in 1987. I also did a manual one on graph paper using a set of coloured felt tips using values calculated on said Dragon 32.
I’d forgotten you once had a Dragon, Rob. I appear to have fallen amongst Tandy CoCo (== Dragon-ish) nerds here in Canada. Of course there’s a still-maintained multitasking OS for it: NitrOS-9
This will get you the exact replication, with the iteration max set at ~200. (I used C#)
for (int iy = 15; iy >= -15; iy–) {
double cy = 1.66 * delta * iy;
for (int ix = -35; ix <= 35; ix++) {
double cx = -0.75 + delta * ix;
Forgot the delta: delta = 0.035;
thanks! Seems to work perfectly