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:

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 tried to do — and come close, but not actually managed to exactly replicate — 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 BRKMAT
! GENERATE FIGURE FROM BROOKS-MATELSKI PAPER C.1978
!  THAT EVENTUALLY BECAME KNOWN AS THE MANDELBROT SET
! - SCRUSS, 2020-06
! 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=70, MAXIT=36)
         CHARACTER*80 OUT
         CHARACTER CH*1

         DO J=1,ROWS
            CI=MAP(REAL(J), 1.0, REAL(ROWS), -0.89, 0.89)
            DO I=1,COLS
               CR=MAP(REAL(I), 1.0, REAL(COLS), -2.0, 0.45)
               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 close:

but not quite right. 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
# brkmat - Brooks-Matelski proto ASCII Mandelbrot set - scruss, 2020-06
# -*- coding: utf-8 -*-

def valmap(value, istart, istop, ostart, ostop):
    return ostart + (ostop - ostart) * ((value - istart) / (istop - istart))

rows = 31
cols = 70
maxit = 36

for y in range(rows):
    ci = valmap(float(y + 1), 1.0, float(rows), -0.89, 0.89)
    for x in range(cols):
        cr = valmap(float(x + 1), 1.0, float(cols), -2.0, 0.45)
        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.