使用 Forth 语言绘制 Mandelbrot 集
Mandelbrot 集是最著名的分形图形之一。下面的 Forth 程序用来绘制 Mandelbrot 集。
定点数的表示
采用 Qint.frac 格式的定点数进行计算。int 为整数部分（包含符号位，若有）的位数，frac 为小数部分的位数。在下面的程序中，定点数的宽度为 16 位（与 Forth 中的单精度整数相同），即 int+frac=16。
Qint.frac 格式的定点数也就是以 \(2^{\mathtt{frac}}\) 为分母的分数，存储时只存储分子。显然，这种格式的定点数的加减法与整数运算相同，一个定点数与一个整数的乘法也与两个整数运算相同。计算两个定点数的乘法时，先将它们按整数相乘（中间结果为双精度整数），然后右移 frac 位（也就是除以 \(2^{\mathtt{frac}}\)）。如果相乘后的结果没有溢出，则只保留低 16 位为结果即可。1
定点数计算的 Forth 实现
以 Q6.10 为例。定义
10 constant FRAC
1 frac lshift constant ONE \ 2^10=1024
one 1- constant MASK \ 1023
16 frac - constant INT \ 6
将单精度整数转换为 Qint.frac 定点数使用 frac LSHIFT ( d -- q )；将单精度整数和小数部分（小于 one）转换为 Qint.frac 定点数使用 frac LSHIFT OR ( d_frac d_int -- q )。
将双精度中间结果右移 frac 位，保留低 16 位：
: 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)
;
两个定点数相乘：
: Q* ( q1 q2 -- q1*q2 ) m* qrshift ;
打印定点数：
: 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 ;
测试：
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
Mandelbrot 集的计算
对复平面上的点 \(c\)，计算映射
\[ z_{n+1}=z_n^2+c: \quad z, c \in \mathbb{C}, z_0=0 \]直到 \(|z|>2\). 设定最大迭代次数（如 15），按达到 \(|z|>2\) 或最大迭代次数时的 \(n\) 值为不同颜色在复平面上 \(c\) 点位置画点。
将 \(z\) 写为 \(z_\mathrm{x} + z_\mathrm{y} \mathrm{i}\) 的形式，有：
\[ \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*} \]计算复数的模及比较其是否大于 2（为避免开方，实际比较的是 \(|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 > ;
使用 ANSI 转义序列改变字符的背景色：
: 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 ;
一次迭代：
: 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 )
计算某个 \(c\) 点的 \(n\) 值：
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 ) ;
绘制 Mandelbrot 集时，\(x\) 的范围为 \(-2 \cdots 1\)，\(y\) 的范围为\(-1 \cdots 1.5\) 加适当余量，并合理设置步长：
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
使用二重循环绘制 Mandelbrot 集：
: MANDEL ( -- )
page
ymax ymin do
xmax xmin do
i j mandel-point color space
xstep +loop
cr
ystep +loop
0 color-code ;
绘制下图在 22.1184MHz 的 My4TH 上约用时 2 分 52 秒。
完整的程序如下：
----[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]-------------------------------------------------------
-
除法的计算：一个定点数除以一个整数，与整数运算相同；两个定点数相除时，先将被除数左移
frac位（乘以 \(2^{\mathtt{frac}}\)），然后除以除数。注意中间结果的位数必须足够宽，Forth 中可以用M*/计算。 ↩︎