簡介
這是一個疊代公式,式中的變數都是複數.這是一個大千世界,從他出發可以產生無窮無盡美麗圖案,他是曼德勃羅特教授在二十世紀七十年代發現的.
曼德勃羅集你看上圖中,有的地方像日冕,有的地方像燃燒的火焰,只要你計算的點足夠多,不管你把圖案放大多少倍,都能顯示出更加複雜的局部.這些局部既與整體不同,又有某種相似的地方,好像著夢幻般的圖案具有無窮無盡的細節和自相似性.曼德勃羅特教授稱此為"魔鬼的聚合物".為此,曼德勃羅特在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, 表明沒有任務,程式終止
