examples/mandelbrot.4th

   1 ( --- Complex arithmetic --- )
   2
   3 ( Location of floating point. )
   4 : precision 10000 ;
   5
   6 : sign dup abs / ;
   7
   8 : >scaled
   9     precision 10 / * over
  10     ?dup 0<> if
  11         sign *
  12     then
  13     swap precision * +
  14 ;
  15
  16 ( Redefine multiplication.  Yay forth! )
  17 : * precision */ ;
  18
  19 : c* ( x1 y1 x2 y2 -- x3 y3 )
  20     swap -rot               ( x1 x2 y1 y2 )
  21     2dup * negate           ( x1 x2 y1 y2 -y1y2 )
  22     4 pick 4 pick * +       ( x1 x2 y1 y2 (x1x2-y1y2))
  23     4 roll 2 roll *         ( x2 y1 (x1x2-y1y2) x1y2 )
  24     3 roll 3 roll * +       ( (x1x2-y1y2) (x1y2+x2y1) )
  25 ;
  26
  27 : c+ ( x1 y1 x2 y2 -- x3 y3 )
  28     rot +
  29     -rot +
  30     swap
  31 ;
  32
  33 : csq 2dup c* ;
  34
  35 : cmagsq ( x1 y1 -- mag )
  36     csq abs
  37 ;
  38
  39 ( --- Mandelbrot set calculations  --- )
  40
  41 : iterate ( cr ci zr zi -- cr ci z'r z'i )
  42     csq c+
  43 ;
  44
  45 : inSet? ( cr ci -- res )
  46     DEBUGON
  47
  48     0 0     ( z_0 = 0 )
  49
  50     true    ( flag indicating set membership )
  51     100 0 do
  52         drop
  53
  54         iterate
  55         2dup cmagsq
  56         4 >scaled > if
  57             false ( not in set )
  58             leave
  59         then
  60
  61         true ( maybe in set )
  62     loop
  63
  64     ( Clear z and c, leaving set membership flag )
  65     -rot 2drop -rot 2drop
  66 ;
  67
  68 : xsteps 100 ;
  69 : ysteps 50 ;
  70
  71 ( Draw the Mandelbrot Set!)
  72 : mandel ( x1 y1 x2 y2 -- )
  73
  74     1 pick 4 pick -
  75     2 pick 5 pick do
  76         i 0 inSet? if
  77             42 emit
  78         else
  79             '.' emit
  80         then
  81     dup +loop
  82 ;
  83
  84 ( Clean up - hide non-standard multiplication def. )
  85 hide *