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.