曼德勃羅集
幾何圖形之一
曼德勃羅特集是人類有史以來做出的最奇異,最瑰麗的幾何圖形,曾被稱為“上帝的指紋”。這個點集均出自公式:Zn+1=(Zn)^2+C,對於非線性迭代公式Zn+1=(Zn)^2+C,所有使得無限迭代后的結果能保持有限數值的複數C的集合,構成曼德勃羅集。
這是一個迭代公式,式中的變數都是複數。這是一個大千世界,從他出發可以產生無窮無盡美麗圖案,他是曼德勃羅特教授在二十世紀七十年代發現的。
圖1:曼德勃羅集
圖形是由美國數學家曼徳勃羅特教授於1975年夏天一個寂靜的夜晚,在冥思苦想之餘翻看兒子的拉丁文字典時想到的,起拉丁文的原意是“產生無規則的碎片”。請看如下的圖形產生過程,其中后一個圖均是前一個圖的某一局部放大:
圖2
如下是產生上圖的出發點
圖3
出發點:
實部Real 0.2537269133080432,
虛部Imag.0.000365995381749671135
The width of that screen is9.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表示虛部最大
值.
to
to
2.曼德勃羅特集的計算
顯示曼德勃羅特集是處理位映射圖像的一個例子。首先要對圖像進行計算,且計算量很
大。曼德勃羅特集是複數平面中的點集,當對一個函數迭代計算時,這些點將處於准穩定狀
態(quasi-stable),即將會增加或減少,但不會超過某一限度。通常該函數為:
式中是複數的第次迭代,c是確定該點在複平面中位置的複數值.z的初始
值為0.迭代將一直進行下去,直到z的幅值(向量長度,這裡為22ba+)大於2或者迭代次
數已經達到某種任意的規定的限度.
化簡計算.
用zreal表示z的實部,zimag表示z的虛部。則
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,表明沒有任務,程序終止