简介

曼德勃罗特集的计算与编程

简介

曼德勃罗特集是人类有史以来做出的最奇异,最瑰丽的几何图形.曾被称为“上帝的指纹”。 这个点集均出自公式:Zn+1=(Zn)^2+C,

这是一个迭代公式,式中的变量都是复数.这是一个大千世界,从他出发可以产生无穷无尽美丽图案,他是曼德勃罗特教授在二十世纪七十年代发现的.你看上图中,有的地方象日冕,有的地方象燃烧的火焰,只要你计算的点足够多,不管你把图案放大多少倍,都能显示出更加复杂的局部.这些局部既与整体不同,又有某种相似的地方,好像着梦幻般的图案具有无穷无尽的细节和自相似性.曼德勃罗特教授称此为"魔鬼的聚合物".为此,曼德勃罗特在1988年获得了"科学为艺术大奖".

图形是由美国数学家曼徳勃罗特教授于1975年夏天一个寂静的夜晚，在冥思苦想之余翻看儿子的拉丁文字典是想到的，起拉丁文的原意是“产生无规则的碎片”

请看如下的图形产生过程,其中后一个图均是前一个图的某一局部放大:

如下是产生上图的出发点

出发点:实部 Real 0.2537269133080432 , 虚部 Imag. 0.000365995381749671135 The width of that screen is 9.45e-17

曼德勃罗特集的计算与编程

曼德勃罗特集是易并行计算的一个典型例子.采用分治技术,并行算法设计时分为静态

任务分配和动态任务分配(可用work-pool or processor farm).

1.测试用例

其中底部的数据(real_min, imag_min) to (real_max, imag_max)表示复平面窗口,real_min表示

实部最小值, imag_min表示虚部最小值,real_max表示实部最大值, imag_max表示虚部最大

值.

(-0.84950, 0.21000) to (-0.84860, 0.21090) (0.32000, -0.45000) to (0.50000, 0.05000) (0.26304, 0.00233) to (0.26329, 0.00267)

(-0.63500, 0.68000) to (-0.62500, 0.69000) (-0.46510, -0.56500) to (-0.46470, -0.56460) (-1.50000, -1.00000) to (0.50000, 1.00000)

2.曼德勃罗特集的计算

显示曼德勃罗特集是处理位映射图像的一个例子.首先要对图像进行计算,且计算量很

大.曼德勃罗特集是复数平面中的点集,当对一个函数迭代计算时,这些点将处于准稳定状

态(quasi-stable),即将会增加或减少,但不会超过某一限度.通常该函数为:

zk+1 = zk

2 + c

式中zk+1是复数z = a + bi的第k+1次迭代,c是确定该点在复平面中位置的复数值.z的初始

值为0.迭代将一直进行下去,直到z的幅值(向量长度,这里为22ba+)大于2或者迭代次

数已经达到某种任意的规定的限度.

化简计算.

z2 = a2 - b2 + 2abi

用zreal表示z的实部,zimag表示z的虚部.则

zreal = a2 - b2 + creal

zimag = 2ab + cimag

3.顺序代码

structure complex {

float real;

float image;

};

对一点值的计算并返回迭代次数的例程:

int cal_pixel(complex c)

{ int count, max;

complex z;

float temp, lengthsq;

max = 256;

z.real = 0;

z.imag = 0;

count = 0;

do {

temp = z.real * z.real - z.imag * z.imag + c.real;

z.imag = 2 * z.real * z.imag + c.imag;

z.real = temp;

lengthsqc = z.real * z.real + z.imag * z.imag;

count++;

} while ((lengthsq < 4.0) && (count < max)); //直到z的幅值大于2或者迭代次数达到max

return count;

}

因此,所有给定的曼德勃罗特点必将处在以原点为中心,半径为2的圆中.

计算和显示这些点的代码需要对坐标系统进行一定的缩放来与显示区域的坐标系统相匹配.

假设显示高度为disp_height,宽度为disp_width.而点在显示区域中的位置为(x, y).如果显

示复数平面的这个窗口具有最小值(real_min, imag_min)和最大值(real_max, imag_max),则每

个点需用以下系数加以缩放.

c.real = real_min + x * (real_max - real_min) / disp_width;

c.imag = imag_min + y * (imag_max - imag_min) / disp_height;

设置scale_real = (real_max - real_min) / disp_width;

scale_imag = (imag_max - imag_min) / disp_height;

缩放代码:

for (x = 0; x < disp_width; x++)

for (y = 0; y < disp_height; y++) {

c.real = real_min + ((float) x * scale_real);

c.imag = imag_min + ((float) y * scale_imag);

color = cal_pixel(c);

display(x, y, color);

}

4.并行代码

1)静态任务分配

假定显示区域为640 × 480,在一个进程中要计算10行.即将10 × 640像素变为一组,共48

个进程.为代码如下:

主进程Master

for ( i = 0, row = 0; i < 48; i++, row = row + 10)

send (&row, pi)

for (i = 0; i < (480 * 640); i++)

recv(&c, &color, PANY);

display(c, color);

}

从进程Slave (process i)

recv(&row, Pmaster);

for (x = 0; x < disp_width; x++)

for (y = row; y < (row + 10); y++) {

c.real = real_min + ((float) x * scale_real);

c.imag = imag_min + ((float) y * scale_imag);

color = cal_pixel(c);

send(&c, &color, Pmaster);

}

改进:成组发送数据.减少通信启动次数.先将结果保存在数组中,然后以一个消息将整个

数组发送给主进程.主进程可用一个通配符以任意顺序接收来自从进程的消息.

2)动态任务分配—工作池/处理器农庄(work pool / processor farm)

动态负载平衡可以用工作池方法实现.在我们的问题中,像素集(更确切应该是坐标集)构成

了任务.任务数是固定的,要计算的像素数在计算前是已知的.各个处理器从工作池中请求

下一个像素子集的坐标.

主进程

count = 0;

row = 0;

for (k = 0; k< procno; k++) {

send(&row, pk, data_tag);

count++;

row++;

}

do {

recv(&slave, &r, color, PANY, result_tag);

count--;

if (row 0)

从进程

recv(y, Pmaster, ANYTAG, source_tag); //接收Pmaster发送的第y行的点

while(source_tag == data_tag) { //判断是否还有消息需要处理

c.imag = imag_min + ((float) y * scale_imag);

for (x = 0; x < disp_width; x++) {

c.real = real_min + ((float) x * scale_real);

color = cal_pixel(c);

}

send(&i, &y, color, Pmaster, result_tag); //将所计算的第y行点的color发给Pmaster

recv(y, &r, Pmaster, source_tag); //接收下一任务

} //如果退出while循环, 则一定是source_tag = termiate_tag, 表明没有任务,程序终止