编写Demo程序不仅要求编程人员具备熟练的数学技巧、卓越的计算机绘图水平,而且还非常考验他们手写x86代码的基本功。
本节介绍的程序都是MS-DOS环境下的.COM程序。
在参考资料[a12]里,我们找到了一个非常简单的随机图案生成程序。虽然它的功能只是不停地在屏幕上打印斜杠和反斜杠,但是这两种字符最终可构成一种几何图案,如图83.1所示。
图83.1 一种随机几何图案
在16位的x86平台上,这种算法非常多。
Trixter在他的网站上[1]公开了具备42字节大小的程序。笔者把它摘录出来,并标注上了自己的注释:
00000000: B001 mov al,1 ; set 40x25 video mode
00000002: CD10 int 010
00000004: 30FF xor bh,bh ; set video page for int 10h call
00000006: B9D007 mov cx,007D0 ; 2000 characters to output
00000009: 31C0 xor ax,ax
0000000B: 9C pushf ; push flags
; get random value from timer chip
0000000C: FA cli ; disable interrupts
0000000D: E643 out 043,al ; write 0 to port 43h
; read 16-bit value from port 40h
0000000F: E440 in al,040
00000011: 88C4 mov ah,al
00000013: E440 in al,040
00000015: 9D popf ; enable interrupts by restoring IF flag
00000016: 86C4 xchg ah,al
; here we have 16-bit pseudorandom value
00000018: D1E8 shr ax,1
0000001A: D1E8 shr ax,1
; CF currently have second bit from the value
0000001C: B05C mov al,05C ;'\'
; if CF=1, skip the next instruction
0000001E: 7202 jc 000000022
; if CF=0, reload AL register with another character
00000020: B02F mov al,02F ;'/'
; output character
00000022: B40E mov ah,00E
00000024: CD10 int 010
00000026: E2E1 loop 000000009 ; loop 2000 times
00000028: CD20 int 020 ; exit to DOS
实际上,上述程序里的伪随机数是Intel 8253计时器芯片(硬件)回传过来的时间信息。它选用了零号计时器,而时钟决定这个计时器每秒递增18.2次。
向0x43端口发送零字节,相当于发送了“选定#0号(通用)计数器”“计数器的输出持续可读”和“采用二进制计数(返回值是二进制数字,而非BCD码)”这三条指令。[2]
当程序执行POPF指令时,CPU恢复IF标识的同时会恢复终端功能。
在使用IN指令读取数据时,返回值必须写到AL寄存器里,所以后面出现了数据交换指令xchg。
这个程序并没有使用计时器查询确切时间,而是用它来生成伪随机数。因此,我们没有必要屏蔽系统中断。此外,我们只需要返回值的低8位数据,所以读这低8位数据即可。
笔者对Trixter的程序稍作精简,把它改进为27字节的程序:
00000000: B9D007 mov cx,007D0 ; limit output to 2000 characters
00000003: 31C0 xor ax,ax ; command to timer chip
00000005: E643 out 043,al
00000007: E440 in al,040 ; read 8-bit of timer
00000009: D1E8 shr ax,1 ; get second bit to CF flag
0000000B: D1E8 shr ax,1
0000000D: B05C mov al,05C ; prepare '\'
0000000F: 7202 jc 000000013
00000011: B02F mov al,02F ; prepare '/'
; output character to screen
00000013: B40E mov ah,00E
00000015: CD10 int 010
00000017: E2EA loop 000000003
; exit to DOS
00000019: CD20 int 020
MS-DOS系统完全没有内存保护技术的概念。换句话说,应用程序可以任意访问内存地址。不仅如此,在使用LODSB指令从DS:SI读取单字节数据的时候,即使程序没有预先给寄存器赋值也没有问题——它会从任意地址读取一个字节!
Trixter甚至在其网页里[3]推荐不初始化相关寄存器就直接使用LODSB指令。
原文同时建议使用SCASB指令替代LODSB指令,因为前者可以在读取数据的同时,根据数据直接设置标识位。
此外,使用DOS系统的syscall INT 29h还能对程序进行进一步精简。这个中断可以在屏幕上输出AL寄存器里的字符。
Peter Ferrie和Andrey“hermlt”Baranovich分别写出了11字节和10字节的程序。[4]
指令清单83.1 Andrey“herm1t”Baranovich: 11 bytes
00000000: B05C mov al,05C ;'\'
; read AL byte from random place of memory
00000002: AE scasb
; PF = parity(AL - random_memory_byte) = parity(5Ch - random_memory_byte)
00000003: 7A02 jp 000000007
00000005: B02F mov al,02F ;'/'
00000007: CD29 int 029 ; output AL to screen
00000009: EBF5 jmp 000000000 ; loop endlessly
SCASB指令计算“AL-[随机地址]”,并设置相应标识位。JP指令比较少见,触发JP转移的条件是奇偶标识位PF为1(偶数)。在这个程序里,输出的字符不再由随机字节的某个比特位决定,而是由这个字节的各个比特位共同决定。因此,这个程序的散列程度有望更好一些。
如果使用x86未公开的SALC指令(即SETALC)、单指令完成“SET AL CF”,那么整个程序还可以更短。这个指令最初出现在NEC v20 CPU上。用自然语言解释的话,它的功能就是“有CF标识位填充AL寄存器”若CF为1,则AL的值将会是0xFF;否则AL的值就是零。受到SALC适用性的影响,任8086/8088平台上这个程序应该跑不起来。
指令清单83.2 Peter Ferrie: 10 bytes
; AL is random at this point
00000000: AE scasb
; CF is set according subtracting random memory byte from AL.
; so it is somewhat random at this point
00000001: D6 setalc
; AL is set to 0xFF if CF=1 or to 0 if otherwise
00000002: 242D and al,02D ;'-'
; AL here is 0x2D or 0
00000004: 042F add al,02F ;'/'
; AL here is 0x5C or 0x2F
00000006: CD29 int 029 ; output AL to screen
00000008: EBF6 jmps 000000000 ; loop endlessly
因此,完全有可能彻底抛弃条件转移指令。反斜杠(“\”)和斜杠(“/”)的ASCII值分别是0x5C和0x2F。余下的问题就是:可否根据CF的(伪随机)状态把AL寄存器的值设置为0x5C和0x2F。
实际上解决方法十分简单:无论AL的值是0还是0xFF,我们把它和0x2D进行“与”运算,即可得到0和0x2D。再把这个值与0x2F相加,就得到了0x2F和0x5C。然后把AL寄存器里的值打印到屏幕上即可。
应当注意的是:在DOSBox、Windows NT主机、甚至是MS-DOS主机上运行同一个程序,看到的图案都可能是不同的。影响程序结果的因素有:模拟器对Intel 8253计时器的不同模拟方式、寄存器的不同初始值,以及其他因素。
多少年来,编程人员不懈地钻研曼德博集合[5]的各种算法。本文将要介绍的,是Sir_Lagsalot在2009年发表的曼德博集合的Demo程序[6]。这个程序由30个16位x86指令构成,文件大小仅为64字节。
它绘制的图案如图83.2所示。
图83.2 曼德博集合Demo程序绘制的图案
本节将介绍这个程序的工作原理。
复数
复数是二元有序实数对,由实部(Re())和虚部(Im())两部分构成。在复数概念的二维空间里,任意复数的实部和虚部都可表示为二维坐标,从而把该复数表示为平面的一个点。
本节用到的复数运算的基本法则有:
加法:(a+bi)+(c+di)=(a+c)+(b+d)i。
即:
Re(sum) = Re(a) + Re(b)
Im(sum) = Im(a) + Im(b)
乘法:(a + bi)(c + di) = (ac − bd) + (bc + ad)i。
即:
Re(product) = Re(a) · Re(c)−Re(b) · Re(d)
Im(product) = Im(b) · Im(c) + Im(a) · Im(d)
平方:(a+bi)2 = (a+bi)(a+bi)=(a2 −b2)+(2ab)i。
即:
Re(square) = Re(a)2−Im(a)2
Im(square) = 2 · Re(a) · Im(a)
曼德博集合的绘制方法
曼德博集合可以由复二次多项式来定义:对于由复数Z() 构成的递归序列来说,不同的参数c可能使序列的绝对值逐渐发散到无限大,也可能收敛在有限的区域内。曼德博集合就是使序列不延伸至无限大的所有复数c的集合。
简单地说,普通程序的做法大致如下:[7]
把屏幕划分为像限/取值区域。
将每个坐标点视为一个c值,并验证它是否属于曼德布洛特集合。
如果这个点所代表的c值不会使递归序列的值逃逸到无限大,那么就在屏幕上用某种颜色标注这个点。
笔者编写了两个程序,分别以复数(宏观参数)和代数二项式(表达形式)两种角度实现上述算法。
指令清单83.3 For complex numbers
def check_if_is_in_set(P):
P_start=P
iterations=0
while True:
if (P>bounds):
break
P=P^2+P_start
if iterations > max_iterations:
break
iterations++
return iterations
# black-white
for each point on screen P:
if check_if_is_in_set (P) < max_iterations:
draw point
# colored
for each point on screen P:
iterations = if check_if_is_in_set (P)
map iterations to color
draw color point
当参数为复数时,就要使用代数二项式“代入”上述参数,并应用前文介绍过的复数计算法则。
指令清单83.4 For integer numbers
def check_if_is_in_set(X, Y):
X_start=X
Y_start=Y
iterations=0
while True:
if (X^2 + Y^2 > bounds):
break
new_X=X^2 - Y^2 + X_start
new_Y=2*X*Y + Y_start
if iterations > max_iterations:
break
iterations++
return iterations
# black-white
for X = min_X to max_X:
for Y = min_Y to max_Y:
if check_if_is_in_set (X,Y) < max_iterations:
draw point at X, Y
# colored
for X = min_X to max_X:
for Y = min_Y to max_Y:
iterations = if check_if_is_in_set (X,Y)
map iterations to color
draw color point at X,Y
维基百科的网站介绍了一种以C#语言实现的曼德博集合算法[9]。不过,那个源程序只能用符号显示迭代次数的信息。笔者把它改得更为直观了一些,让它能够显示出序列的迭代次数[10]:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Mnoj
{
class Program
{
static void Main(string[] args)
{
double realCoord, imagCoord;
double realTemp, imagTemp, realTemp2, arg;
int iterations;
for (imagCoord = 1.2; imagCoord >= -1.2; imagCoord -= 0.05)
{
for (realCoord = -0.6; realCoord <= 1.77; realCoord += 0.03)
{
iterations = 0;
realTemp = realCoord;
imagTemp = imagCoord;
arg = (realCoord * realCoord) + (imagCoord * imagCoord);
while ((arg < 2*2) && (iterations < 40))
{
realTemp2 = (realTemp * realTemp) - (imagTemp * imagTemp) - realCoord;
imagTemp = (2 * realTemp * imagTemp) - imagCoord;
realTemp = realTemp2;
arg = (realTemp * realTemp) + (imagTemp * imagTemp);
iterations += 1;
}
Console.Write("{0,2:D} ", iterations);
}
Console.Write("\n");
}
Console.ReadKey();
}
}
}
运行上述程序,可得到文件:http://beginners.re/examples/mandelbrot/result.txt
。
这个程序限定迭代次数的上限为40。输出结果中的“40”表示在40次迭代之内序列仍然收敛;而40以内的数字(比方说是n),则表示在该点的c值在迭代n次之后令序列逃逸。
除此之外,http://demonstrations.wolfram.com/MandelbrotSetDoodle/
公开了一个精彩的demo程序。它能够用彩色线条显示指定c值产生的递归序列zn。如果彩线落在灰色圆圈(模为2)之外,那么该点代表的c值就不属于曼德博集合。
首先,笔者在黄色区域之内选取一个c值,观察它产生的递归序列,如图83.3所示。
图83.3 在黄色区域之内选取c值
上述线条是收敛的。这表明笔者选取的c值属于曼德博集合。
然后,笔者在黄色区域之外选取c值。线条立刻混乱、甚至超出边界。如图83.4所示。
图83.4 在黄色区域之外选取c值
这种图案是发散的。这表示刚才选取的c值不属于曼德博集合。
这个网站还开发了另一款曼德博集合的研究程序,有兴趣的读者可访问:http://demonstrations. wolfram.com/IteratesForTheMandelbrotSet/
。
本节将要介绍的这个曼德博集合的demo程序,仅由30条指令、64个字节构成。虽然它的算法与前文的算法一致,但是它用到了很多编程技巧。
整个程序的源代码是公开的。为了便于读者阅读,笔者添加了一些注释。
指令清单83.5 程序源代码+注释
1 ; X is column on screen
2 ; Y is row on screen
3
4
5 ; X=0, Y=0 X=319, Y=0
6 ; +------------------------------->
7 ; |
8 ; |
9 ; |
10 ; |
11 ; |
12 ; |
13 ; v
14 ; X=0, Y=199 X=319, Y=199
15
16
17 ; switch to VGA 320*200*256 graphics mode
18 mov al,13h
19 int 10h
20 ; initial BX is 0
21 ; initial DI is 0xFFFE
22 ; DS:BX (or DS:0) is pointing to Program Segment Prefix at this moment
23 ; ... first 4 bytes of which are CD 20 FF 9F
24 les ax,[bx]
25 ; ES:AX=9FFF:20CD
26
27 FillLoop:
28 ; set DX to 0. CWD works as: DX:AX = sign_extend(AX).
29 ; AX here 0x20CD (at startup) or less then 320 (when getting back after loop),
30 ; so DX will always be 0.
31 cwd
32 mov ax,di
33 ; AX is current pointer within VGA buffer
34 ; divide current pointer by 320
35 mov cx,320
36 div cx
37 ; DX (start_X) - remainder (column: 0..319); AX - result (row: 0..199)
38 sub ax,100
39 ; AX=AX-100, so AX (start_Y) now is in range -100..99
40 ; DX is in range 0..319 or 0x0000..0x013F
41 dec dh
42 ; DX now is in range 0xFF00..0x003F (-256..63)
43
44 xor bx,bx
45 xor si,si
46 ; BX (temp_X)=0; SI (temp_Y)=0
47
48 ; get maximal number of iterations
49 ; CX is still 320 here, so this is also maximal number of iteration
50 MandelLoop:
51 mov bp,si ; BP = temp_Y
52 imul si,bx ; SI = temp_X*temp_Y
53 add si,si ; SI = SI*2= (temp_X*temp_Y)*2
54 imul bx,bx ; BX = BX^2= temp_X^2
55 jo MandelBreak ; overflow?
56 imul bp,bp ; BP = BP^2= temp_Y^2
57 jo MandelBreak ; overflow?
58 add bx,bp ; BX = BX+BP = temp_X^2 + temp_Y^2
59 jo MandelBreak ; overflow?
60 sub bx,bp ; BX = BX-BP = temp_X^2+temp_Y^2 - temp_Y^2 = temp_X^2
61 sub bx,bp ; BX = BX-BP = temp_X^2 - temp_Y^2
62
63 ; correct scale:
64 sar bx,6 ; BX=BX/64
65 add bx,dx ; BX=BX+start_X
66 ; now temp_X = temp_X^2 - temp_Y^2 + start_X
67 sar si,6 ; SI=SI/64
68 add si,ax ; SI=SI+start_Y
69 ; now temp_Y = (temp_X*temp_Y)*2 + start_Y
70
71 loop MandelLoop
72
73 MandelBreak:
74 ; CX=iterations
75 xchg ax,cx
76 ; AX=iterations. store AL to VGA buffer at ES:[DI]
77 stosb
78 ; stosb also increments DI, so DI now points to the next point in VGA buffer
79 ; jump always, so this is eternal loop here
80 jmp FillLoop
这个程序的算法是:
切换到分辨率为320×200/256色的VGA模式。320×200 = 64000(0xFA00)。256色的每个像素都是单字节的数据,所以缓冲区大小就是0xFA00字节。应用程序使用ES:DI寄存器对即可对像素进行寻址。
VGA图形缓冲区的段地址要存储在ES寄存器里,所以ES寄存器的值必须是0xA000。但是向ES寄存器传递0xA000的指令至少要占用4个字节(PUSH 0A000h/POP ES)。有关6位MS-DOS系统的内存模型,可参见第94章的详细介绍。
假设BX寄存器的值为零、程序段前缀(Program Segment Prefix)位于第0号地址处,那么2字节的LES AX,[BX] 指令就会在AX寄存器里存储0x20CD、在ES寄存器里存储0x9FFF。也就是说,这个程序会在图形缓冲区之前输出16个像素(字节)。但是因为该程序的运行平台是MS-DOS,而MS-DOS没有实现内存保护技术,所以这种问题不会引发程序崩溃。不过,屏幕右侧多出了一个16像素宽的红色条带,整个图像向左移动了16个像素。这就是在程序里节省2字节空间的代价。
单个循环处理全部像素。在处理图像问题时,一般的程序都会使用两个循环:一个循环遍历x坐标、另一个循环遍历y坐标。不过,在VGA图形缓存区里对某个像素进行定位时,双循环的程序必须通过乘法运算才能寻址。为此,这个程序的作者决定使用单循环的数据结构,通过除法运算获取当前点的坐标。经过转换以后,程序坐标的取值范围是:x~[−256,63],y~[−100,99]。正因如此,整个图像的中心点偏右。虽然直接把x的值减去160就可以使x的取值范围变成[−160,159],从而校准图像中心,但是“SUB DX,160”的指令占用4个字节,而作者的“DEC DH”(相当于SUB DX,0x100)只用两个字节。图像中心偏右算是节省文件空间的另一个代价吧。
检测当前点是否属于曼德博集合。检测算法和前文相同。
以CX寄存器作为循环计数器进行LOOP循环。作者没有明确地给循环计数器赋值,而是让CX寄存器继续沿用第35行的数值(320)。或许,数值越大越精确吧。这同样可以节省文件空间。
使用IMUL指令。因为操作符是有符号数,所以乘法指令不是MUL。同理,为了把原点坐标0,0调整到屏幕中心区域,在转换坐标时(除法)使用到的是SAR指令(带符号右移),而不是SHR指令。
简化逃逸检测的算法。常规算法检测的是坐标,即一对坐标值。而作者分三次检测了溢出问题:两个求平方操作和一次加法运算的溢出。事实正是如此,我们使用的是16位寄存器,寄存器内的数值不会超过[−32768,+32767]。在计算机进行有符号数的乘法运算时,只要坐标的值大于32767,那么该点的复数就不会属于曼德博集合。
后面再次出现的SAR指令(相当于除以64)设置了64级素灰度。调色版的灰度值越大,对比度就越高、图像就越清晰;反之,灰度越集中,对比度就越低、图像就越模糊。
程序执行到MandelBreak标签的情况分为两种:一种情况是循环结束时CX = 0这就说明,该点属于曼德博集合;另一种情况是程序发生溢出、此时CX存储着非零值。这个程序把CX的低8位(即CL)写到图形缓冲区里。在默认调色板中,0代表黑色,所以属于曼德博集合的坐标点都会显示纯黑色像素。当然,我们确实能够在绘图之前把调色板设置得更个性化一些,不过那种程序就不会只有64字节了!
这个程序采用的是无限循环,不能正常退出。如果还要判断循环终止或者进行互动响应,那么程序文件还要更大。
此外,这个程序的优化技巧同样值得一提:
使用单字节的CWD指令清空DX寄存器。相比之下,“XOR DX,DX”占用2个字节,而“MOV DX,0”则占用3个字节。
使用单字节的“XCHG AX,CX”完成双字节“MOV AX,CX”的操作。毕竟后面不再使用AX寄存器,所以交换数据完全不会造成问题。
因为程序没有对DI寄存器(图形缓冲区的位置)进行初始化赋值,所以它在启动时的初始值为0xFFFE(寄存器初始值由操作系统决定)。不过这无伤大雅,这个程序只要求DI的值保持在[0,0xFFFF]之间,软件用户也看不到屏幕区域之外的点(在320×200×256的图像缓冲区里,最后一个像素的地址是0xF9FF)。也就是说,因为这个程序前后衔接得严丝合缝,所以不管DI寄存器也没问题。否则,编程人员还要添加把DI寄存器置零、以及检测图形缓冲区结束边界的指令。
指令清单83.6 My“fixed”version
1 org 100h
2 mov al, 13h
3 int 10h
4
5 ; set palette
6 mov dx, 3c8h
7 mov al, 0
8 out dx, al
9 mov cx, 100h
10 inc dx
11 l00:
12 mov al, cl
13 shl ax, 2
14 out dx, al ; red
15 out dx, al ; green
16 out dx, al ; blue
17 loop l00
18
19 push 0a000h
20 pop es
21
22 xor di, di
23
24 FillLoop:
25 cwd
26 mov ax, di
27 mov cx, 320
28 div cx
29 sub ax, 100
30 sub dx, 160
31
32 xor bx, bx
33 xor si, si
34
35 MandelLoop:
36 mov bp, si
37 imul si, bx
38 add si, si
39 imul bx, bx
40 jo MandelBreak
41 imul bp, bp
42 jo MandelBreak
43 add bx, bp
44 jo MandelBreak
45 sub bx, bp
46 sub bx, bp
47
48 sar bx, 6
49 add bx, dx
50 sar si, 6
51 add si, ax
52
53 loop MandelLoop
54
55 MandelBreak:
56 xchg ax, cx
57 stosb
58 cmp di, 0FA00h
59 jb FillLoop
60
61 ; wait for keypress
62 xor ax, ax
63 int 16h
64 ; set text video mode
65 mov ax, 3
66 int 10h
67 ; exit
68 int 20h
笔者修正了上一个程序的缺陷:使用了平滑过渡的灰度调色板;把整个图形全部输出到了图形缓冲区里(第19、20行));使图像中心与屏幕中心重合(第30行);绘图结束后等待键盘敲击再退出程序(第58~68行)。不过,整个程序大了近一倍:它由54条指令构成,文件大小也增长到了105字节。该程序的绘图如图83.5所示。
图83.5 笔者的改进版程序的绘图
[1] http://trixter.oldskool.org/2012/12/17/maze-generation-in-thirteen-bytes/。
[2] 实际上,8位控制字相当于四条指令。因为另一个指令没有在本程序发挥作用,所以本文没有进行介绍。
[3] http://trixter.oldskool.org/2012/12/17/maze-generation-in-thirteen-bytes/。
[4] 请参照http://pferrie.host22.com/misc/10print.htm。
[5] 还被译作曼德布洛特集合,英文原文是Mandelbrot set。
[6] http://www.pouet.net/prod.php?which=53287。
[7] 微软较为全面地介绍了数学方面和编程方面的细节,建议读者阅读:https://msdn.microsoft. com/zh-cn/library/jj635753%28v=vs.85%29.aspx。
[8] 实际算法都令z0=0,因此前几步运算在计算z1。请参考前面推荐的微软文档。
[9] https://en.wikipedia.org/wiki/Mandelbrot_set。
[10] 修改版程序的下载地址是:http://beginners.re/examples/mandelbrot/dump_ iterations.exe。