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
  36 : conj ( x y -- x -y )
  37     negate
  38 ;
  39
  40 : cmagsq ( x1 y1 -- mag )
  41     2dup conj c* +
  42 ;
  43
  44 ( --- Mandelbrot set calculations  --- )
  45
  46 : iterate ( cr ci zr zi -- cr ci z'r z'i )
  47     2over 2swap csq c+
  48 ;
  49
  50 : inSet? ( cr ci -- res )
  51     0 0     ( z_0 = 0 )
  52
  53     true    ( flag indicating set membership )
  54     100 0 do
  55         drop
  56
  57         iterate
  58         2dup cmagsq
  59         4 0 >scaled > if
  60             false ( not in set )
  61             leave
  62         then
  63
  64         true ( maybe in set )
  65     loop
  66
  67     ( Clear z and c, leaving set membership flag )
  68     -rot 2drop -rot 2drop
  69 ;
  70
  71 : xsteps 100 ;
  72 : ysteps 30 ;
  73
  74 ( Draw the Mandelbrot Set!)
  75 : mandel ( x1 y1 x2 y2 -- )
  76
  77     0 pick 3 pick - ysteps /
  78     1 pick 4 pick do
  79
  80         2 pick 5 pick - xsteps /
  81         3 pick 6 pick do
  82
  83             i j inSet? if
  84                 42 emit
  85             else
  86                 space
  87             then
  88
  89         dup +loop
  90         drop
  91
  92         cr
  93
  94     dup +loop
  95     drop
  96 ;
  97
  98 ( Clean up - hide non-standard multiplication def. )
  99 hide *