13
$\begingroup$

This question is not profound, but I can't figure this out myself and thought I'd ask here. Although the paper is of great historical importance, I don't think History of Science and Mathematics Stackexchange is appropriate.

Is there any historical information about the grid size used in the first rendering of the Mandelbrot set shown in The Dynamics of 2-Generator Subgroups of PSL(2, C), Robert Brooks and J. Peter Matelski, 1978? I'm just trying to reproduce this pattern and while I can get close, I can't quite nail it.

Since this is from about 40 years ago, computing time was several orders of magnitude slower (laptops are GigaFlops), so I understand I may need to play with the number of iterations. Also I haven't switched to higher precision yet (I'm just using python's float). But before I get too involved in that, I'd at least like to know I'm using the same grid points as they are.

EDIT: Ideally the answer would be the actual numbers known to be used by the authors when generating this historic image. But it seems more fun to deduce them, so either way is allowed.

Stopping at 1000 iterations, points with markers maintaned $\vert z \vert< 2$, just for example:

$\hskip3.3cm$enter image description here

enter image description here

$\endgroup$
10
  • 1
    $\begingroup$ I guess you could count the dots from $-2$ to $1/4$ to figure out the resolution. $\endgroup$
    – Mark McClure
    Jul 14, 2017 at 17:53
  • 2
    $\begingroup$ @uhoh Sorry, I didn't mean to imply that it's a pointless question. I was just curious why it would be important. If it's just for historicity - I guess that's a good enough reason. $\endgroup$
    – Jair Taylor
    Jul 14, 2017 at 18:40
  • 1
    $\begingroup$ @uhoh I've also spent time rendering the Mandelbrot set (once even on a TI-83 calculator) and I can confirm that it is indeed fun. $\endgroup$
    – Jair Taylor
    Jul 14, 2017 at 18:52
  • 1
    $\begingroup$ It looks like there is a grid point at $(-0.75, 0)$, otherwise there would be a thicker bridge between the cardioid and the circle - this gives the X alignment of the grid. The antenna visibility shows the Y alignment of the grid is at 0. Counting dots (inaccurate as you say, so posting as a comment rather than an answer) across the circle gives around $0.035$ for the X spacing, the Y spacing looks around $0.055$ - so the grid is not square. $\endgroup$
    – Claude
    Jul 14, 2017 at 19:06
  • 1
    $\begingroup$ @Claude thanks! I hadn't even though about using different scales for $x$ and $y$, I wonder if the grid was pre-scaled so that it would come out square after printing on the line printer? (6 lines, 10 characters per inch perhaps?) That would be fun to know! In fact, fun reigns, so let's drop the restriction of factual numbers and open it up to allow reverse-engineered values as well. I'll edit the question accordingly. $\endgroup$
    – uhoh
    Jul 14, 2017 at 19:19

3 Answers 3

Reset to default
15
$\begingroup$

I can confirm that the programming constants used were indeed $\Delta x = .035$ , $\Delta y = 1.66 * \Delta x$ . The grid center is $(-.75 , 0)$ which is a boundary point by the way: it is exactly represented in floating point. The original iteration cap was 250.

The run date was 8 Jan 1979. It was my first computer $M$ image using the now standard program. There has been lots of criticism of it, but it is clearly fine.

The font Lucida-console regular has aspect ratio very close to 1.66 Take Notepad or WinVi with this font and you have an effective previewer and printer in imitation of the lineprinter at Stony Brook in 1978.

Note: the UNIVAC mainframe had 36 bit single and 72 bit double precision. Of course, I ran both varieties. The modern choices are 32, 64, or 80 bit. Using 32 bit is problematic.

$\endgroup$
6
  • 6
    $\begingroup$ Thank you for your confirmation, what a wonderful first post @JPM! $\endgroup$
    – uhoh
    Jan 15, 2020 at 21:37
  • $\begingroup$ To confirm, are you the actual author of the original image? $\endgroup$
    – Neinstein
    Jan 14, 2021 at 9:49
  • 1
    $\begingroup$ @Neinstein, if they answer your question, I should hope it to be negative -- given that Mr. Mandelbrot died a bit over 10 years ago. Although, never say never ... $\endgroup$
    – noughtnaut
    Jan 14, 2021 at 13:29
  • 1
    $\begingroup$ @Noughtnaut Oh. aha. hm.... The sentence "It was my first computer 𝑀 image using the..." , plus the lot of uncited insight, kind of made me suspicious. Perhaps a young contributor? $\endgroup$
    – Neinstein
    Jan 14, 2021 at 16:52
  • 8
    $\begingroup$ @Neinstein OP is asking about a 1978 paper by Robert Brooks and J. Peter Matelski. $\endgroup$
    – att
    Feb 8, 2021 at 23:46
10
$\begingroup$

I wrote a small Haskell program using the Diagrams library:

import Diagrams.Prelude hiding (aspect)
import Diagrams.Backend.SVG.CmdLine (B, defaultMain)

import Data.Complex

main :: IO ()
main = defaultMain (diagram # centerXY # bg white # lw thin)

diagram :: Diagram B
diagram = vcat . map hcat $
  [ [ if mandelbrot (x :+ y) then asterisk else space
    | i <- [-36 .. 33], let x = spacing * i - 0.75
    ]
  | j <- [-16 .. 16], let y = spacing * j * aspect
  ]

aspect :: Double
aspect = 1.66

spacing :: Double
spacing = 0.035

iterations :: Int
iterations = 200

mandelbrot :: Complex Double -> Bool
mandelbrot c
  = null
  . dropWhile (<= 2)
  . map magnitude
  . take iterations
  . iterate (\z -> z^2 + c)
  $ 0

asterisk :: Diagram B
asterisk = withEnvelope space $ mconcat
  [ p2 (-2, -2) ~~ p2 (2,  2)
  , p2 (-2,  2) ~~ p2 (2, -2)
  , p2 ( 0, -3) ~~ p2 (0,  3)
  ]

space :: Diagram B
space = phantom' (rect 6 10)

phantom' :: Diagram B -> Diagram B
phantom' = phantom

Here is the output:

first Mandelbrot set image recreation

I found the important magic values aspect, spacing and iterations by trial and improvement: I used the spacing from dot counting, and first an aspect of $\frac{10}{6}$, but that was a little too high to get the right shape at the top and bottom so I reduced it a little bit by bit (9.99/6, 9.98/6, ...). Finally I tuned the iterations using binary search to get the remaining pixels the same as the image in the question (and I initially made a mistake, thanks for the correction in the comments).

$\endgroup$
4
  • 3
    $\begingroup$ Well done!..... $\endgroup$
    – lhf
    Jul 14, 2017 at 23:05
  • 1
    $\begingroup$ Thank you! I've heard many fond references to Haskell but to see it here helps explain why it's so well liked by those who use it. Noticing the $x=-0.75$ point and unequal spacing right off was the key, sweet! Is 180 roughly the minimum number of iterations for a stable pattern? $\endgroup$
    – uhoh
    Jul 15, 2017 at 4:45
  • 1
    $\begingroup$ @uhoh thanks, well spotted, I fixed it now $\endgroup$
    – Claude
    Jul 21, 2017 at 13:17
  • 1
    $\begingroup$ Beautiful code <3 $\endgroup$
    – Jay
    Jan 13, 2021 at 14:17
3
$\begingroup$

Claude's solution in Haskell is indeed beautiful, but that language wasn't available at the time the paper was written. In communication with one of the authors elsewhere, FORTRAN 66 was used. The following code would be more like the code that ran on the authors' UNIVAC mainframe:

C     *M*-SET PLOT FROM BROOKS AND MATELSKI PAPER
C     OS/8 FORTRAN IV FOR DEC PDP-8 - SCRUSS, 2023-04
C     THANKS TO JPM FOR PARAMETERS
C     HTTPS://MATH.STACKEXCHANGE.COM/A/3510304/578897
C     LPT IS ON LOGICAL UNIT 3
      INTEGER IX,IY,K
      DOUBLE PRECISION CR,CI,ZR,ZI,TI,D
      D=.035D0
      DO 40 IY=-15,15
C     FORTRAN IV CANNOT COUNT BACKWARDS, SO NEGATE Y
      CI=-1.66D0*D*IY
      DO 30 IX=-35,35
      CR=-.75D0+D*IX
      ZR=0D0
      ZI=0D0
      DO 10 K=1,200
      TI=2D0*ZR*ZI
      ZR=ZR*ZR-ZI*ZI+CR
      ZI=TI+CI
      IF (ZR*ZR+ZI*ZI-4D0) 10,10,20
 10   CONTINUE
C     PRINT '*' IF C IS IN SET
      WRITE (3,50)
      GO TO 30
C     PRINT SPACE IF C IS NOT IN SET
 20   WRITE (3,60)
 30   CONTINUE
C     PRINT NEW LINE
      WRITE (3,70)
 40   CONTINUE
C     FORTRAN CARRIAGE CONTROL: '+' SUPPRESSES NEW LINE, ' ' DOES NOT
C     'H' DENOTES HOLLERITH CONSTANTS: OBSOLETE STRING HANDLING
 50   FORMAT (2H+*)
 60   FORMAT (2H+ )
 70   FORMAT (1H )
      END

The code can be compiled with a modern Fortran compiler (use gfortran -std=legacy), but the mechanics of converting the output from ASA carriage control is left as an exercise.

                                                     *
                                                   ****
                                                  ******
                                                   *****
                                              * *********
                                         *** ****************
                                           ******************** **
                                        *************************
                                       ****************************
                                     *******************************
                                     ******************************
                        * *****     ********************************
                       ***********  ********************************
                      ************* *******************************
                  ** ************** ******************************
****************************************************************
                  ** ************** ******************************
                      ************* *******************************
                       ***********  ********************************
                        * *****     ********************************
                                     ******************************
                                     *******************************
                                       ****************************
                                        *************************
                                           ******************** **
                                         *** ****************
                                              * *********
                                                   *****
                                                  ******
                                                   ****
                                                     *

While most dialects of the language supported a COMPLEX numeric type which would have made the code much clearer, the author confirmed that they hadn't used that feature. Note, too, the use of D in double precision constants. If these are omitted, single precision constants are converted to doubles, resulting in inaccuracies.

$\endgroup$
3
  • $\begingroup$ Thanks! That's interesting the PDP-8 could not reproduce the result, sounds like something that could be explored further in another site. I learned on FORTRAN, punch cards and line printers myself so happy to see something familiar :-) How would you calculate the Mandelbrot set in hardware? would have been interesting but it was closed during EESE's "darker ages". $\endgroup$
    – uhoh
    Apr 7 at 20:12
  • 1
    $\begingroup$ @uhoh - that was a mistake on my part. It seems that explicitly stating double precision constants (such as D=.035D0) is required for correct operation. I've corrected the code, and the PDP-8 produces the same results as an x86_64 Linux machine $\endgroup$
    – scruss
    Apr 8 at 23:48
  • $\begingroup$ excellent, good catch - all is well in the universe again :-) $\endgroup$
    – uhoh
    Apr 8 at 23:50

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .