The Mandelbrot Set, before Mandelbrot

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:

the squashed bug that started it all: the Mandelbrot set

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:

text plot of the Mandelbrot set with points inside set marked with asterisks
figure 2 (original caption “The set of C’s such that f(z) = z² + C has a stable periodic orbit“) from 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.

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:

mandelbrot set, rendered in asterisks from BRKMTF.F
party like it’s ’78

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.

Update: see the comment from one of the authors of the original paper below.

7 comments

  1. 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.

  2. 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

  3. 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;

  4. To nobody: if you are going to quote my recent disclosure you should give the reference: Grid spacing, iterations used in the 1978 first published rendering of the Mandelbrot set?. It took 40+ years for this to appear in print. To scruss: FORTRAN was supplied by the hardware manufacturer Rand-UNIVAC. 1978 was way too soon to see FORTRAN-77. I did not buffer the output—-one computation one pixel. No complex support used either, real and imaginary parts sufficed. You need to dust off your WRITE(6,10)
    10 FORMAT(2H+*) and 20 FORMAT(2H+ ). By the way single precision was 36 bits and double was 72 bits.

Leave a comment

Your email address will not be published. Required fields are marked *