%%HP: T(3)A(R)F(.);
DIR                             
@ FRACTALS, by Dan Ciarniello
@
@ These programs will generate the Mandelbrot set.  To speed up calculations, a
@ couple of properties of the Mandelbrot set are used.
@
@ Since the Mandelbrot set is symmetric about the x-axis, if the x-axis is
@ present, the program will generate the image on only one side and then map
@ the image to the other side of the axis.
@
@ The Mandelbrot set is a closed connected set.  This means that if all points
@ on the border are in the Mandelbrot set, then all points interior to the
@ rectangle are also in the Mandelbrot set.  Ditto if all points are not in the
@ Mandelbrot set.  Thus the program calculates all points around the border of
@ a rectangle.  If all points on the border are the same (in or out), then fill
@ the rectangle the appropriate colour and go to the next rectangle.
@ Otherwise, cut the rectangle in half and check the border of the new
@ rectangle.  This is called the Mariani technique.  Note that loop counters
@ run through pixel coordinates.  These pixel coordinates are converted to user
@ coordinates using PX\->C which uses PPAR to determine conversion factors.

  CST { MANDEL SAVE               @ Custom menu.
LOAD { DITH                       @ DITH turns dithering on
    \<< 3 SF                      @ NODITH turns dithering off
    \>> } { NODITH                @                     
    \<< 3 CF                      @                     
    \>> } }                       @                     
  LOAD                            @ LOAD restores a previously saved
    \<< 2 MENU                    @ image and its associated PPAR.
"Enter Filename" ""               @
INPUT OBJ\-> OBJ\->               @
DROP PICT STO                     @
'PPAR' STO 0 MENU                 @
GRAPH                             @
    \>>                           @
  SAVE                            @ SAVE an image along with its PPAR.
    \<< PPAR PICT RCL             @ PPAR is necessary so that the correct
2 \->LIST                         @ user units are restored when the
"Enter Filename" {                @ image is restored with LOAD.
"" \Ga } INPUT OBJ\->             @
STO                               @
    \>>                           @
  MANDEL                          @ The main program
    \<< ERASE { # 0d              @ Erase PICT and view it as program
# 0d } PVIEW RCLF                 @ runs; save flag status; start
TICKS \-> c1 c2 f t               @ timing; save initial window corner
      \<< 4 CF 5 CF               @ values.
        IF 'nITTR'                @ Check whether nITTR exists.  If
VTYPE -1 ==                       @ not, create it with a default value
        THEN 100                  @ of 100.
'nITTR' STO                       @
        END                       @
        IF c2 IM c1               @ Correct the aspect ratio of the
IM - c2 RE c1 RE -                @ image.
64 * 131 / >                      @
        THEN c2 IM                @ If y2-y1 > (64/131)(x2-x1) then set
c1 IM - 131 * 64 /                @ x-value of upper right corner (x2) to
c1 RE + c2 IM                     @ x1+(y2-y1)*131/64
        ELSE c2 RE                @
DUP c1 RE - 64 *                  @ If y2-y1 < (64/131)(x2-x1) then set
131 / c1 IM +                     @ y2 to y1+(x2-x1)*64/131
        END R\->C                 @
'c2' STO                          @
        IF 'PPAR'                 @
VTYPE -1 ==                       @ Check whether PPAR exists.  If not,
        THEN {                    @ create it.                          
(-2,-1) (2.09375,1)               @                                     
X 0 (100,100) TRUTH               @                                     
Y } 'PPAR' STO                    @
        END PPAR 2                @                                     
c2 PUT 'PPAR' STO                 @ Store the window corners in PPAR to
PPAR 1 c1 PUT                     @ set the user coordinates.           
'PPAR' STO                        @
        IF c1 IM c2               @                                     
IM * 0 <                          @ Determine whether the x-axis is     
        THEN c1 RE                @ present.  If so determine the pixel
0 R\->C C\->PX c2 RE 0            @ coordinates of the left and right   
R\->C C\->PX \-> lside            @ end.                                
rside                             @                                     
          \<<                     @                                     
            IF c2                 @ Determine whether the x-axis is in  
IM c1 IM ABS >                    @ the upper or lower half of the      
            THEN 0                @ screen.  If its in the lower half,  
130                               @ then generate the image from top of
              FOR x               @ the screen down to the x-axis.      
x R\->B # 0d 2 \->LIST            @ ITERATE around the first rectangle  
PX\->C ITERATE                    @ starting at the top of the screen.  
              NEXT                @                                     
1 rside 2 GET B\->R               @                                     
              FOR y               @ ITERATE down right side to x-axis.  
# 130d y R\->B 2                  @                                     
\->LIST PX\->C ITERATE            @
              NEXT                @
129 0                             @
              FOR x               @ ITERATE along the x-axis from right
x R\->B lside 2 GET 2             @ to left.                            
\->LIST PX\->C ITERATE            @                                     
-1                                @                                     
              STEP                @                                     
lside 2 GET B\->R 1               @                                     
              FOR y               @ ITERATE up the left side of the     
# 0d y R\->B 2 \->LIST            @ screen from the x-axis to the top   
PX\->C ITERATE -1                 @ of the screen.                      
              STEP                @                                     
{ # 0d # 0d } rside               @
MARIANI 63 lside 2                @ Pass the upper left and lower right
GET B\->R                         @ pixel coordinates of the rectangle  
              FOR y               @ to MARIANI.                         
# 0d y R\->B 2 \->LIST            @                                     
PICT OVER PX\->C CONJ             @ Map the area above the x-axis to    
C\->PX DUP 1 # 130d               @ that below the x-axis line by line.
PUT SUB PICT 3                    @
ROLLD REPL -1                     @
              STEP                @
            ELSE 0                @
130                               @ If the x-axis is in the upper half  
              FOR x               @ of the screen, generate the image   
x R\->B lside 2 GET 2             @ from the x-axis to the bottom of    
\->LIST PX\->C ITERATE            @ the screen.                         
              NEXT                @ Start with the x-axis from left to  
lside 2 GET 1 + B\->R             @ right.                              
63                                @                                     
              FOR y               @                                     
# 130d y R\->B 2                  @ Down left side to bottom of screen.
\->LIST PX\->C ITERATE            @                                     
              NEXT                @                                     
129 0                             @                                     
              FOR x               @                                     
x R\->B # 63d 2 \->LIST           @ Along bottom of screen right to     
PX\->C ITERATE -1                 @ left.                               
              STEP                @                                     
62 lside 2 GET 1 +                @
B\->R                             @                     
              FOR y               @                     
# 0d y R\->B 2 \->LIST            @ Up right side to x-axis.
PX\->C ITERATE -1                 @                     
              STEP                @                     
lside { # 130d                    @ Pass upper left and lower right     
# 63d } MARIANI 0                 @ corners of rectangle to MARIANI.    
lside 2 GET B\->R                 @                                     
              FOR y               @ Map area below the x-axis to that   
# 0d y R\->B 2 \->LIST            @ above the x-axis line by line.      
PICT OVER PX\->C CONJ             @                                     
C\->PX DUP 1 # 130d               @                                     
PUT SUB PICT 3                    @                                     
ROLLD REPL                        @                                     
              NEXT                @                                     
            END                   @                                     
          \>>                     @                                     
        ELSE 1 130                @                                     
          FOR x x                 @ If the x-axis is not present, just  
R\->B # 0d 2 \->LIST              @ ITERATE around the edges of the     
PX\->C ITERATE                    @ screen to start.                    
          NEXT 1 63               @ Along top.                          
          FOR y                   @                                     
# 130d y R\->B 2                  @                                     
\->LIST PX\->C ITERATE            @ Down right side                     
          NEXT 129                @                                     
0                                 @                                     
          FOR x x                 @                                     
R\->B # 63d 2 \->LIST             @                                     
PX\->C ITERATE -1                 @ Along bottom                        
          STEP 62 1               @                                     
          FOR y                   @                                     
# 0d y R\->B 2 \->LIST            @                                     
PX\->C ITERATE -1                 @ Up left side                        
          STEP {                  @                                     
# 0d # 0d } {                     @ Pass upper left and lower right     
# 130d # 63d }                    @ coordinates of the screen to        
MARIANI                           @ MARIANI                             
        END TICKS t               @                                     
- B\->R 29491200 / f              @ Done! Calculate execution time.     
STOF                              @ Restore flags.                      
      \>>                         @                                     
    \>>                           @                                     
  MARIANI                         @                                     
@ MARIANI is a recursive routine.
@ It takes a rectangle (defined by its upper left corner (ulc) and lower right
@ corner (rlc), divides it in two and iterates along the dividing line.  It
@ then checks the border of each new rectangle in turn to determine whether all
@ points on the border have the same status (either on or off).  If so, fill
@ the rectangle.  If not, MARIANI calls itself with the coordinates of the
@ smaller rectangle.  This contues until a rectangle is filled or until the
@ current rectangle is less than 6 pixels in area in which case the state of
@ each pixel is determined individually and no further division occurs.
    \<< \-> ulc lrc               @
      \<< lrc 1 GET               @                                     
ulc 1 GET - lrc 2                 @                                     
GET ulc 2 GET - *                 @                                     
        IF # 6d <                 @
        THEN ulc 1                @ If area less than 6 pixels,         
GET 1 + B\->R lrc 1               @ determine state of each pixel in    
GET 1 - B\->R                     @ rectangle individually.             
          FOR x ulc               @                                     
2 GET 1 + B\->R lrc 2             @                                     
GET 1 - B\->R                     @                                     
            FOR y x               @                                     
R\->B y R\->B 2 \->LIST           @                                     
PX\->C ITERATE                    @                                     
            NEXT                  @                                     
          NEXT                    @                                     
        ELSE lrc 1                @ If area greater than 6, divide      
GET ulc 1 GET - lrc               @ rectangle in two.                   
2 GET ulc 2 GET -                 @                                     
          IF <                    @                                     
          THEN lrc                @ If the rectangle is taller than it  
2 GET ulc 2 GET + 2               @ is wide, determine line which       
/ DUP ulc 2 ROT PUT               @ divides it into upper and lower     
SWAP lrc 2 ROT PUT                @ halves.  Set flag 4 to indicate     
4 SF                              @ this.                               
          ELSE lrc                @                                     
1 GET ulc 1 GET + 2               @ If the rectangle is wider than it   
/ DUP ulc 1 ROT PUT               @ is tall, determine the line which   
SWAP lrc 1 ROT PUT                @ divides it into left and right      
          END \->                 @ halves.                             
lines linet                       @                                     
          \<<                     @ lines: start coordinate of line
            IF 4                  @ linet: end coordinate of line       
FS?C                              @                                     
            THEN                  @                                     
lines 1 GET 1 + B\->R             @                                     
linet 1 GET B\->R 1 -             @ If flag 4 set, draw the horizontal  
              FOR x               @ dividing line.                      
lines 1 x R\->B PUT               @                                     
PX\->C ITERATE                    @                                     
              NEXT                @                                     
            ELSE                  @                                     
lines 2 GET 1 + B\->R             @                                     
linet 2 GET B\->R 1 -             @                                     
              FOR y               @ Otherwise draw the vertical one.    
lines 2 y R\->B PUT               @                                     
PX\->C ITERATE                    @
              NEXT                @
            END ulc               @
linet CHECKBORDER                 @ Check the border of the left (or    
            IF 1                  @ upper) rectangle (defined by        
FS? 2 FS? OR                      @ corners ulc and linet).             
            THEN                  @ If either flag 1 or flag 2 is set,  
PICT ulc 1 GET 1 +                @ then all pixels around the border   
ulc 2 GET 1 + 2                   @ of the rectangle have the same      
\->LIST linet 1 GET               @ status.                             
ulc 1 GET - 1 -                   @                                     
linet 2 GET ulc 2                 @
GET - 1 - BLANK                   @
              IF 2                @
FS?                               @ Fill black if flag 2 is set.
              THEN                @ Fill white if flag 1 is set.        
NEG                               @                                     
              END                 @                                     
REPL 1 CF 2 CF                    @                                     
            ELSE                  @                                     
ulc linet MARIANI                 @ If border pixels are not all one    
            END                   @ state, call MARIANI with corners    
lines lrc                         @ ulc and linet.                      
CHECKBORDER                       @                                     
            IF 1                  @ Check border of right (or lower)    
FS? 2 FS? OR                      @ rectangle (defined by lines and     
            THEN                  @ lrc).                               
PICT lines 1 GET 1                @                                     
+ lines 2 GET 1 + 2               @                                     
\->LIST lrc 1 GET                 @                                     
lines 1 GET - 1 -                 @ Fill according to flag status as    
lrc 2 GET lines 2                 @ above.                              
GET - 1 - BLANK                   @                                     
              IF 2                @                                     
FS?                               @                                     
              THEN                @                                     
NEG                               @                                     
              END                 @                                     
REPL 1 CF 2 CF                    @                                     
            ELSE                  @                                     
lines lrc MARIANI                 @ Or call MARIANI with corners lines  
            END                   @ and lrc.                            
          \>>                     @                                     
        END                       @                                     
      \>>                         @                                     
    \>>                           @                                     
  CHECKBORDER                     @                                     
    \<< \-> ulc lrc               @ Determine status of rectangle       
      \<< ulc PIX? \->            @ border.                             
state                             @ First determine state of ulc pixel.
        \<< ulc 1 GET             @                                     
1 + B\->R lrc 1 GET               @                                     
B\->R                             @                                     
          FOR x ulc               @                                     
1 x R\->B PUT PIX?                @ Now check each pixel around border  
            IF                    @ starting with top line and compare
state \=/                         @ to state of ulc.  If they are the
            THEN 5                @ same continue around border         
SF lrc 1 GET B\->R                @ otherwise break out of the loop,    
'x' STO                           @ there is no need to check further.
            END                   @ Set flag 5 to stop further          
          NEXT                    @ checking.                           
          IF 5 FC?                @                                     
          THEN ulc                @
2 GET 1 + B\->R lrc 2             @ If flag 5 still clear check down    
GET B\->R                         @ right side.                         
            FOR y                 @                                     
lrc 2 y R\->B PUT                 @                                     
PIX?                              @                                     
              IF                  @                                     
state \=/                         @                                     
              THEN                @                                     
5 SF lrc 2 GET B\->R              @                                     
'y' STO                           @
              END                 @                                     
            NEXT                  @                                     
          END                     @                                     
          IF 5 FC?                @                                     
          THEN ulc                @ If flag 5 still clear, check left   
2 GET 1 + B\->R lrc 2             @ side.                               
GET B\->R                         @                             
            FOR y                 @                             
ulc 2 y R\->B PUT                 @                             
PIX?                              @                             
              IF                  @                             
state \=/                         @                             
              THEN                @                             
5 SF lrc 2 GET B\->R              @                             
'y' STO                           @
              END                 @                             
            NEXT                  @                             
          END                     @                                     
          IF 5 FC?                @                                     
          THEN ulc                @                                     
1 GET 1 + B\->R lrc 1             @ If flag 5 still clear, check        
GET 1 - B\->R                     @ bottom.                             
            FOR x                 @
lrc 1 x R\->B PUT                 @
PIX?                              @
              IF                  @
state \=/                         @
              THEN                @
5 SF lrc 1 GET B\->R              @
'x' STO                           @
              END                 @
            NEXT                  @
          END                     @
          IF 5 FC?C               @
          THEN                    @ If flag 5 clear then border all one
state 1 + SF                      @ state.  Set flag 1 for all pixels   
          END                     @ off.  Set flag 2 for all pixels     
        \>>                       @ on.                                 
      \>>                         @
    \>>                           @
  ITERATE                         @ ITERATE determines whether a        
    \<< DUP DUP PIXON             @ coordinate is in the Mandelbrot set
PIXOFF C\->R                      @ or not.  Toggle the current pixel   
0 \-> x y i                       @ as a visual aid in following        
      \<< 0 0 0 0                 @ progress of program.                
        WHILE SQ                  @ Initialize Z(0) as (0,0) and        
SWAP SQ SWAP DUP2 +               @ duplicate.  Real and imaginary      
4 < 'i' INCR nITTR                @ parts are treated separately on the
< AND                             @ stack.  The loop calculates         
        REPEAT - x                @ Z(n+1)=SQ(Z(n))+c where c=x+iy.     
+ 3 ROLLD * DUP + y               @ It repeats until SQ(Z)>4 or nITTR   
+ DUP2                            @ iterations are completed.           
        END 4 DROPN               @                                     
x y R\->C i                       @                                     
        IF nITTR \>=              @                                     
        THEN PIXON                @ If nITTR iterations are completed,  
        ELSE                      @ the coordinate is in the Mandelbrot
          IF i 20 <               @ set.  Turn the pixel on.            
3 FS? AND                         @                                     
          THEN i 2                @ If dithering is on, check to see    
MOD                               @ whether the point diverged in less  
            IF                    @ than 20 iterations.  If so and it   
            THEN                  @ is odd, turn the pixel on.          
PIXON                             @
            ELSE                  @
PIXOFF                            @
            END                   @
          ELSE                    @
PIXOFF                            @
          END                     @
        END                       @
      \>>                         @
    \>>                           @
END                               @
