( --- Complex arithmetic --- )
( Location of floating point. )
-: precision 10000 ;
+10000 value precision
-: >scaled precision 10 / * swap precision * + ;
+: sign dup abs / ;
+
+: >scaled
+ precision 10 / * over
+ ?dup 0<> if
+ sign *
+ then
+ swap precision * +
+;
( Redefine multiplication. Yay forth! )
: * precision */ ;
: c* ( x1 y1 x2 y2 -- x3 y3 )
- swap -rot ( x1 x2 y1 y2 )
- 2dup * negate ( x1 x2 y1 y2 -y1y2 )
- 4 pick 4 pick * + ( x1 x2 y1 y2 (x1x2-y1y2))
- 4 roll 2 roll * ( x2 y1 (x1x2-y1y2) x1y2 )
- 3 roll 3 roll * + ( (x1x2-y1y2) (x1y2+x2y1) )
+ swap -rot ( x1 x2 y1 y2 )
+ 2dup * negate ( x1 x2 y1 y2 -y1y2 )
+ 4 pick 4 pick * + ( x1 x2 y1 y2 (x1x2-y1y2))
+ 4 roll 2 roll * ( x2 y1 (x1x2-y1y2) x1y2 )
+ 3 roll 3 roll * + ( (x1x2-y1y2) (x1y2+x2y1) )
;
: c+ ( x1 y1 x2 y2 -- x3 y3 )
- rot +
- -rot +
- swap
+ rot +
+ -rot +
+ swap
;
: csq 2dup c* ;
+
+: conj ( x y -- x -y )
+ negate
+;
+
: cmagsq ( x1 y1 -- mag )
- csq abs
+ 2dup conj c* +
;
( --- Mandelbrot set calculations --- )
: iterate ( cr ci zr zi -- cr ci z'r z'i )
- csq c+
+ 2over 2swap csq c+
;
+100 value maxiter
+
: inSet? ( cr ci -- res )
+ 0 0 ( z_0 = 0 )
- 100 0 do
+ true ( flag indicating set membership )
+ maxiter 0 do
+ drop
- 2swap 2dup 5 roll 5 roll
iterate
2dup cmagsq
- 100 > if
+ 4 0 >scaled > if
+ false ( not in set )
leave
then
+ true ( maybe in set )
loop
+
+ ( Clear z and c, leaving set membership flag )
+ -rot 2drop -rot 2drop
;
+100 value xsteps
+30 value ysteps
+
( Draw the Mandelbrot Set!)
-: mandel ( x1 y1 x2 y2 -- )
+: mandeldraw ( x1 y1 x2 y2 -- )
+ 0 pick 3 pick - ysteps /
+ 1 pick 4 pick do
+
+ 2 pick 5 pick - xsteps /
+ 3 pick 6 pick do
+
+ i j inSet? if
+ 42 emit
+ else
+ space
+ then
+
+ dup +loop
+ drop
+
+ cr
+
+ dup +loop
+ drop
;
( Clean up - hide non-standard multiplication def. )
hide *
+
+( Default picture )
+: mandel
+ -2 0 >scaled -1 0 >scaled 0 5 >scaled 1 0 >scaled
+ mandeldraw
+;
+
+." Enter 'mandel' to draw the Mandelbrot Set."