Plotting the Mandelbrot set with Forth

Posted on • 1419 words • 3 minute read

Mandelbrot set is one of the most famous fractal images. The following Forth program is used to plot the Mandelbrot set.

Representation of Fixed-Point Numbers

Fixed-point numbers in Qint.frac format are used for computation. int is the number of bits for the integer part (including the sign bit, if any), and frac is the number of bits for the fractional part. In the following program, the width of the fixed-point number is 16 bits (the same as a single-precision integer in Forth), i.e., int+frac=16.

Fixed-point numbers in Qint.frac format are essentially fractions with a denominator of \(2^{\mathtt{frac}}\), storing only the numerator. Addition and subtraction of fixed-point numbers in this format are the same as integer operations, and multiplication of a fixed-point number by an integer is also the same as multiplication of two integers. When multiplying two fixed-point numbers, they are first multiplied as integers (the intermediate result is a double-precision integer), and then right-shifted by frac bits (i.e., divided by \(2^{\mathtt{frac}}\)). If the multiplied result does not overflow, the lower 16 bits are retained as the result.1

Fixed-Point Arithmetic Implementation in Forth

Using Q6.10 as an example. Define

10 constant FRAC
1 frac lshift constant ONE  \ 2^10=1024
one 1- constant MASK        \ 1023
16 frac - constant INT      \ 6

To convert a single-precision integer to a Qint.frac fixed-point number, use frac LSHIFT ( d -- q ); to convert a single-precision integer and a fractional part (less than one) to a Qint.frac fixed-point number, use frac LSHIFT OR ( d_frac d_int -- q ).

Shift the double-precision intermediate result right by frac bits, retaining the lower 16 bits:

: QRSHIFT   ( L H -- H_L>>frac )
   swap     ( H L )
   frac rshift  ( H L>>frac )
   swap     ( L>>frac H )
   mask and ( L>>frac H_lower_frac_bits )
   int lshift or \ result is (L>>frac)|(H_lower_frac_bits<<int)
   ;

Multiplication of two fixed-point numbers:

: Q* ( q1 q2 -- q1*q2 ) m* qrshift ;

Print a fixed-point number:

: Q. ( q -- ) dup abs ( q |q| )
   <#
      dup mask and 10000 one */ 0 ( q |q| d_|q|_frac_decimal )
      # # # # 2drop \ print frac part in decimal
      [char] . hold \ print decimal point
      frac rshift 0 ( q |q| d_|q|_int ) #s \ print int part
      rot sign  \ print sign
   #> type space ;

Test:

0 q. 0.0000  ok
32767 q. 31.9990  ok
32768 q. -32.0000  ok
one one 1 rshift + q. 1.5000  ok
one one 1 rshift + . 1536  ok
1536 1536 q* q. 2.2500  ok
1536 -1536 q* q. -2.2500  ok
-1536 1536 q* q. -2.2500  ok
-1536 -1536 q* q. 2.2500  ok

Computation of the Mandelbrot Set

For a point \(c\) on the complex plane, compute the mapping

\[ z_{n+1}=z_n^2+c: \quad z, c \in \mathbb{C}, z_0=0 \]

until \(|z|>2\). Set a maximum number of iterations (e.g., 15), and plot a point at the location \(c\) on the complex plane with a color corresponding to the value of \(n\) when either \(|z|>2\) is reached or the maximum iteration count is attained.

Expressing \(z\) in the form \(z_\mathrm{x} + z_\mathrm{y} \mathrm{i}\), we have:

\[ \begin{align*} {z_\mathrm{x}}_{n+1} = & {z_\mathrm{x}}_n^2 - {z_\mathrm{y}}_n^2 + c_\mathrm{x}, \\ {z_\mathrm{y}}_{n+1} = & 2 {z_\mathrm{x}}_n {z_\mathrm{y}}_n + c_\mathrm{y}. \end{align*} \]

Compute the magnitude of a complex number and compare whether it is greater than 2 (to avoid square root, actually compare \(|z|^2>4\)):

: QSQR ( q -- q*q ) dup q* ;
: Q-MAG-SQR ( qa qb -- qmag^2 ) qsqr swap qsqr + ;
4 one * constant FOUR
: MAG2> ( qa qb -- is_mag_>_2 ) q-mag-sqr four > ;

Using ANSI escape sequences to change the background color of characters:

: COLOR-CODE ( c -- ) \ print color escape codes
  0 <# [char] m hold #s [char] [ hold 27 hold #> type ;
: COLOR ( 0~15 -- )
   dup 8 and if
      7 and 100     \ 8~15: 100...107
   else 
      40            \ 0~7: 40...47
   then  
   + color-code ;

One iteration：

: MANDEL-ITER   ( qzx qzy qcx qcy -- qzx1 qzy1 )
   2swap 2dup   ( qcx qcy qzx qzy qzx qzy )
   q* 2*        ( qcx qcy qzx qzy 2*qzx*qzy )
   3 roll +     ( qcx qzx qzy qzy1=2*qzx*qzy+qcy )
   -rot         ( qcx qzy1 qzx qzy )
   qsqr swap    ( qcx qzy1 qzy^2 qzx )
   qsqr swap    ( qcx qzy1 qzx^2 qzy^2 )
   -            ( qcx qzy1 qzx^2-qzy^2 )
   rot +        ( qzy1 qzx1=qzx^2-qzy^2+qcx )
   swap ;       ( qzx1 qzy1 )

Calculate the value of \(n\) at a certain point \(c\):

16 constant MAX-ITER
: MANDEL-POINT ( qcx qcy -- iter )
   0 0 0            ( qcx qcy qzx0=0 qzy0=0 iter0=0 )
   max-iter 0 do    ( qcx qcy qzx qzy iter )
      drop          ( qcx qcy qzx qzy )
      2dup mag2>    ( qcx qcy qzx qzy is_mag_>_2 )
      if 
         i leave    ( qcx qcy qzx qzy iter ) 
      then          ( qcx qcy qzx qzy )
      2over         ( qcx qcy qzx qzy qcx qcy )
      mandel-iter   ( qcx qcy qzx1 qzy1 ) 
      i ( qcx qcy qzx1 qzy1 iter )
   loop
   ( qcx qcy qzxn qzyn iter ) 
   nip nip nip nip ( iter ) ;

When plotting the Mandelbrot set, the \(x\) range is \(-2 \cdots 1\), and the \(y\) range is \(-1 \cdots 1.5\) with appropriate margins. The step sizes are set appropriately:

one 2 * negate constant XMIN
one constant XMAX
one 24 / constant XSTEP
one 8 / constant YSTEP
one negate ystep 2* - constant YMIN
one one 2 rshift + ystep + constant YMAX

Plot the Mandelbrot set with a double loop:

: MANDEL ( -- )
   page
   ymax ymin do
      xmax xmin do
         i j mandel-point color space
      xstep +loop
      cr
   ystep +loop
   0 color-code ;

Plotting the image below took about 2 minutes and 52 seconds on the My4TH at 22.1184MHz.

Mandelbrot Set

The complete program is as follows:

----[200]-------------------------------------------------------
\ C Mandelbrot Set

\ Qint.frac format fixed-point number
10 constant FRAC   1 frac lshift constant ONE
one 1- constant MASK  16 frac - constant INT
: QRSHIFT  ( L H -- H_L>>frac )
   swap frac rshift swap mask and int lshift or ;
: Q* ( q1 q2 -- q1*q2 ) m* qrshift ;
: Q. ( q -- ) dup abs  <# dup mask and 10000 one */ 0 # # # #
   2drop  [char] . hold  frac rshift 0 #s rot sign #>
   type space ;
\ utils for Mandelbrot set calculation
: QSQR ( q -- q*q ) dup q* ;
: Q-MAG-SQR ( qa qb -- qmag^2 ) qsqr swap qsqr + ;
4 one * constant FOUR
: MAG2> ( qa qb -- is_mag_>_2 ) q-mag-sqr four > ;           -->
----[201]-------------------------------------------------------
\ iteration and point calculation for the Mandelbrot set

: MANDEL-ITER ( qzx qzy qcx qcy -- qzx1 qzy1 )
   \ qzx1 = qzx^2 - qzy^2 + qcx
   \ qzy1 = 2*qzx*qzy + qcy
   2swap 2dup q* 2* 3 roll + ( qcx qzx qzy qzy1 )
   -rot qsqr swap qsqr swap - rot + swap ; ( qzx1 qzy1 )
16 constant MAX-ITER
: MANDEL-POINT ( qcx qcy -- iter )
   0 0 0 max-iter 0 do  drop \ drop i in new iter
      2dup mag2> if i leave then ( qcx qcy qzx qzy )
      2over mandel-iter ( qcx qcy qzx1 qzy1 ) i loop
   ( qcx qcy qzxn qzyn i ) nip nip nip nip ( i ) ;


                                                             -->
----[202]-------------------------------------------------------
\ plot the Mandelbrot set

: COLOR-CODE ( c -- ) \ print color escape codes
  0 <# [char] m hold #s [char] [ hold 27 hold #> type ;
: COLOR ( 0~15 -- )
   dup 8 and  if 7 and 100 else 40 then  + color-code ;

one 2 * negate constant XMIN  one constant XMAX
one 24 / constant XSTEP  one 8 / constant YSTEP
one negate ystep 2* - constant YMIN
one one 2 rshift + ystep + constant YMAX

: MANDEL ( -- )  page
   ymax ymin do  xmax xmin do
      i j mandel-point color space
   xstep +loop  cr  ystep +loop  0 color-code ;
----[EOF]-------------------------------------------------------

  1. Division: dividing a fixed-point number by an integer is the same as integer division; when dividing two fixed-point numbers, the dividend is first left-shifted by frac bits (multiplied by \(2^{\mathtt{frac}}\)), and then divided by the divisor. Note that the intermediate result must have sufficient bit width; in Forth, M*/ can be used for the calculation. ↩︎