- 积分
- 15
- 注册时间
- 2008-7-20
- 仿真币
-
- 最后登录
- 1970-1-1
|
发表于 2009-9-8 14:01:29
|
显示全部楼层
来自 黑龙江哈尔滨
1# FreddyMusic
Mandelbrot 编译? 分形? 多少次迭代用了2秒?绘图时多少个点? 机器大概配置?
凑个热闹
-
- with(DocumentTools):
- kernelopts( gcfreq=10^7 ):
- # Computes number of iterations it takes for each point in a strip of the complex plane to escape (if at all)
- # from the Mandelbrot set. Stops when the absolute value of the iterate is greater than the bailout value,
- # or when the maximum number of iterations is exceeded. Operates inplace on imageArray.
- MandelLoop := proc( X, Y, imageArray, i_low, i_high, j_low, j_high, iter, bailout )
- local i, j, Xc, Yc, Xtemp, Ytemp, Xold, Yold, k;
- option hfloat;
- for i from i_low to i_high do
- for j from j_low to j_high do
- Xtemp := X[i];
- Ytemp := Y[j];
- Xc := Xtemp;
- Yc := Ytemp;
- k := 0;
- while k < iter do
- Xold := Xtemp;
- Yold := Ytemp;
- Xtemp := Xold^2-Yold^2+Xc;
- Ytemp := 2*Xold*Yold+Yc;
- if Xtemp^2+Ytemp^2 >= bailout then
- imageArray[i,j] := k;
- break;
- end if;
- k := k+1;
- end do
- end do;
- end do;
- end proc:
- #If we're using Maple 13 or higher then multithread calculations, otherwise singlethread the calcualtions
- if parse(convert( kernelopts(version), string )[7..8]) >=13 then
- # This procedure recursively divides the complex plane into two strips, and creates a Task for each strip.
- # These Tasks are sent to seperate cores for processesing (or enter a queue if the cores are busy)
- # Strips with 5 rows or less are sent to MandelLoop().
- MandelTask := proc ( X, Y, imageArray, i_low, i_high, j_low, j_high, iter, bailout )
- local i, j, i_mid;
- if ( i_high - i_low > 5 ) then
- i_mid := floor( (i_high-i_low)/2 ) + i_low;
- Threads:-Task:-Continue( null,
- Task=[MandelTask, X, Y, imageArray, i_low, i_mid, j_low, j_high, iter, bailout],
- Task=[MandelTask, X, Y, imageArray, i_mid+1, i_high, j_low, j_high, iter, bailout ] );
- else
- MandelLoop( X, Y, imageArray, i_low, i_high, j_low, j_high, iter, bailout );
- end if;
- end proc:
- # Main entry function.
- # numPoints - number of discretization points along real and imaginary axis on complex plane
- # iter - number of iterations for each discretized point on the complex plane
- # X1 and X2 - left and right point on the real axis
- # Y1 and Y2 - top and bottom point on imaginary axis
- # bailout - value of iterate above which a point is no longer considered to be in the Mandelbrot set
- Mandelbrot := proc (numPoints::integer,iter::integer, X1::float, X2::float, Y1::float, Y2::float, bailout::float)
- local X::Vector, Y::Vector, imageArray::Array, i::int, j::int:
- X := Vector(numPoints, i->X1+(X2-X1)*(i-1)/(numPoints-1) , datatype = float[8]);
- Y := Vector(numPoints, i->Y1+(Y2-Y1)*(i-1)/(numPoints-1) , datatype = float[8]);
- imageArray := Array(1 .. numPoints, 1 .. numPoints, datatype = float[8]);
- Threads:-Task:-Start( MandelTask, X, Y, imageArray, 1, numPoints, 1, numPoints, iter, bailout );
- return imageArray;
- end proc:
- DocumentTools[SetProperty](threadedLabel,caption,cat("Multi-threading across ", kernelopts(numcpus), " CPUs"));
- else
- Mandelbrot := proc ( numPoints::integer,iter::integer, X1::float, X2::float, Y1::float, Y2::float, bailout::float)
- local X::Vector, Y::Vector, imageArray::Array:
- X := Vector(numPoints, i->X1+(X2-X1)*(i-1)/(numPoints-1) , datatype = float[8]);
- Y := Vector(numPoints, i->Y1+(Y2-Y1)*(i-1)/(numPoints-1) , datatype = float[8]);
- imageArray := Array(1 .. numPoints, 1 .. numPoints, datatype = float[8]);
- MandelLoop( X, Y, imageArray, 1, numPoints, 1, numPoints, iter, bailout );
- return imageArray;
- end proc:
- DocumentTools[SetProperty](threadedLabel,caption,"Single-threaded mode");
- end if;
- #############################################
- #Set Default Values
- defaultNumPoints:=100:
- defaultIterations:=50:
- defaultxMandStart:=-2.0:
- defaultxMandEnd:=0.7:
- defaultyMandStart:=-1.35:
- defaultyMandEnd:=1.35:
- defaultBailout:=10.0:
- defaultNumPoints:=100:
- defaultIterations:=50:
- defaultSmoothing:=false:
- defaultGrayScaleColoring:=true:
- SetProperty(xMandStartLabel, caption, defaultxMandStart);
- SetProperty(xMandEndLabel, caption, defaultxMandEnd);
- SetProperty(yMandStartLabel , caption, defaultyMandStart);
- SetProperty(yMandEndLabel, caption, defaultyMandEnd);
- SetProperty(bailoutBox, value, defaultBailout);
- SetProperty(pointsBox, value, defaultNumPoints);
- SetProperty(iterationsBox, value, defaultIterations);
- SetProperty(smoothingCB,value,defaultSmoothing);
- SetProperty(grayScaleRB,value, defaultGrayScaleColoring);
- ############################################
- #############################################
- #Initialize the plot
- bailout:=defaultBailout:
- iterations:=defaultIterations:
- numPoints:=defaultNumPoints:
- smoothing:=defaultSmooting:
- if defaultGrayScaleColoring=true then coloring:=none: else coloring:=HUE: end if:
- xMandStart:=defaultxMandStart:
- xMandEnd:=defaultxMandEnd:
- yMandStart:=defaultyMandStart:
- yMandEnd:=defaultyMandEnd:
- points := Mandelbrot(numPoints, iterations, xMandStart, xMandEnd, yMandStart, yMandEnd, bailout);
- MandelbrotPlot:= plots[listdensityplot](points,style=patchnogrid,smooth=defaultSmoothing,colorstyle=coloring,axes=none);
- Do(%Plot1 = MandelbrotPlot);
- #############################################
- hasBeenDragged:=0:
复制代码 |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
|