Mandelbrot集(5)：硬件升级，算法优化：一切为了速度！ （yanlb2000, 2007.05.14, yanlb2000.blogcn.com） 上文，我简要分析了一下绘制Mandelbrot集的算法。提到，这个算法，主要可以归结为两个特点：浮点数密集型算法；适合编制为高度的并行运算。 以上分析的目的也是很明确的，就是要优化算法，提高绘制M集的速度。 这几年来，电脑绘制M集的速度的确有了飞速的发展。总结下来，大致可归结为几种类型方面的进步：有的纯粹是硬件运算能力的提高，有的是硬件能力和软件优化的配合，还有就是纯算法方面的优化。 ·当初我是在286上绘制M集图像的。80286是个16位的CPU，主频一般是12MHz到24MHz，只能进行整数运算，因为标配是没有数学协处理器的。而所有的浮点运算，都是通过软件来实现的。主频本来就低，又没有硬件级别的浮点数运算，所以，速度慢是必然的。记得当时画一个M集，640*480的分辨率，迭代限制为256次，居然需要约2个小时。所以，我们只能等待更快的CPU的出现。 ·此后的若干年，CPU严格遵循着摩尔定律发展着，386，486，Pentium, Pentium II, III, IV等相继问世，主频也一路高升，超过了几百MHz，达到了3GHz甚至更高。 而从80386DX开始，及以后所有型号x86系列，CPU就内置数学协处理器了。数学协处理器对M集算法的重要性是不言而喻的。 因此，单是主频的提升，以及硬件级的浮点数运算能力，就让该算法的实现速度大大提升了。 这些完全是硬件方面的进步。 ·1996年左右，Intel推出了带MMX扩展的Pentium处理器。MMX即是多媒体扩展，能够在同一个指令中，并行处理多个数据。这是Intel针对当时大量出现的图像、声音密集型应用而设计的。其实从技术上说，这应属于原来大型系统CPU中才有的并行运算技术SIMD(Single Instruction Multiple Data)，只是因为集成电路的迅猛发展，能够在民用CPU上实现了。 我当时就发觉，这对M集绘制算法是有很大的意义的。不过遗憾的是，MMX只对整形数提供了并行支持，不支持浮点数。但，这毕竟是一个方向，我相信对浮点数的SIMD的支持也会实现的。 果然，Intel随后在Pentium III等CPU中实现了SSE技术，能够提供对浮点数的并行运算。随着CPU的发展，更完善的SSE2， SSE3等指令集也相继出现。 M集算法中，每个点的运算过程都是完全相同而且互不干扰的，所以，这特别适宜于通过SSE指令来实现。当然，真正要发挥SSE的威力，还需要程序语言、编译器甚至操作系统等的支持和配合才行。 ·当奔腾升级到Pentium IV，主频上升到3GHz多接近4GHz的时候，无论是Intel还是AMD，都遇到了频率的瓶颈。CPU的主频已经很难象以前那样容易提升了。此后，CPU主要向着多核的方向发展，一个CPU封装中，有多个物理CPU核心。这样，电脑就可同时执行多个线程或进程，达到真正的并行处理。这种硬件进步，对M集算法也是有直接好处的。算法可以将绘制区域分块，交给多个线程，在多个CPU核心上并行处理。 当然，实现多核运算，也必须要编程语言、编译器、操作系统的支持。 这里，需要说明一下多核运算和上面SSE指令运算的异同。它们都是通过并行执行多个运算，来达到提升总体运算速度的目的。但在实现上是有区别的。 多核CPU简单地说就是一个电脑中就安装了多个物理的CPU核心，多个核心并行工作。根据算法特点，这些核心可以不同步，甚至执行其他程序指令。 而以上SSE指令则是在同一个CPU核心中，通过调用一条运算指令，就对多个数据同时进行运算。执行的是相同的一条指令。这，也正是将SSE归结为SIMD(Single Instruction Multiple Data，单指令多重数据)的原因。 而且，以上两种并行优化方式互不影响，可以同时叠加进行！ 上面说的是硬件方面的提升和优化。再来看看纯软件算法方面，是否有可改进的地方。 ·观察M集全图，并仔细考察算法，可以发现，当对M集内部的那些区域（集图像中心黑色葫芦形部分）进行计算的时候，每个绘图点所需的运算时间是最长的。因为，无论进行多少次迭代，都不能达到绝对值超出限定值的这个“溢出”条件，所以，总会进行“满负荷”的运算。而显然，最后的结果总是相同的，因为都是M集，所以，它们都属于同一种颜色（一般都画为黑色）。 所以，我们完全有理由对该区域进行优化。我们可以设定一个保守的M集区域。为简化算法，这可以是一个由几个矩形组合成的区域。每次运算的时候，先判断该绘图点是否属于该区域，如果是，则不进入迭代运算，而直接给出属于M集的结果。这样，就避免了大量的浮点运算。我那时在286上的Turbo C程序，就进行了这种优化，的确是大大缩短了绘图时间。 不过，设定这个优化区，是需要非常小心的，不能因为贪图更多的优化而靠边缘太近，否则一些本该有图案的地方会“黑”掉。 ·优化无止境，我们进一步讨论。 当画好一副M集图像之后，我们往往希望能够“深入”其中的一些局部，放大观察。而且，最好是类似“漫游”一样，实时动画显示，能够平滑地实现放大、移动等动画过程。这一定是个美妙的过程。显然，为了达到实时动画显示的效果，我们至少需要每秒十几或最好二十多个画面的刷新率(fps)。（作为参考，NTSC电视画面是每秒30帧，PAL电视画面是每秒25帧。）而即使是现在的CPU，要达到每秒算出几十个画面的M集，还是有困难的。所以，必须大幅度优化。 考察动画过程，相邻两个画面，其实大部分是相同的，只是区域增大或减少了一小部分。所以，可以考虑保留上个画面中的大部分计算结果，而只对放大或缩小的部分进行运算，将得到的结果，与上个画面的结果进行插值运算，修正之后，就可以得到一个新的画面了。 说详细一点，举个例子。比如，放大一个640*480的画面，假设其下一个画面，只保留了原画面中间的部分，上、下、左、右的各2行都将被“放大”到屏幕外。这样，原中间的636*476将被放大到640*480。所以，只需要计算新“扩展出”的大致4行4列的绘图区域的值，再将它们均匀插值到剩下的636*476中就可以了。 这其实不是一个“不失真”的优化算法，因为这与混沌的初始值敏感的特性是背道而驰的。一个绘图点的值，通过上述优化计算得到的结果，可能与真正计算的有出入，甚至较大的出入。但实际上，只要设计得好，这种误差是可以保持在一个很小级别内的。 但单单通过这种优化，我们就可以获得极大的速度提升，我们可以将运算量降低大致2个数量级，这是非常可观的！   最后，再谈谈绘图输出的问题。M集计算之后，还需要通过绘图设备画出来。最普通最方便的方法，当然就是在显示器上画出来（废话！）。但这里也是有讲究的。当绘制一个大画面或全屏幕画面的时候，像素点达到了几十万或几百万，这也不得不涉及到速度问题。如果使用传统的Windows GDI作图方式，那是非常低效的。画单个静态画面还勉强能胜任，但如果需要达到上面所说的实时动画绘图，GDI是达不到要求的。好在，自从Windows 9x之后，微软就一直在发展DirectX规范，包括能绕过低效的GDI，直接访问显示卡绘图能力的DirectDraw，DirectGraphic。这样，实时绘图输出也就不是问题了。  

Continue reading about Mandelbrot集(5)：硬件升级，算法优化：一切为了速度！

Mandelbrot集(4)：算法分析   （yanlb2000, 2007.05.11, yanlb2000.blogcn.com）   继续讨论Mandelbrot集。 今天，来进一步分析一下绘制M集的算法。分析，是为了更多地了解掌握此算法。而且，在此基础上，更可探讨如何优化算法。   绘制M集，其实就是找出屏幕上每个绘图点对应的复数C，然后计算其对应的颜色就可以了。而这个颜色，又是由该绘图点所对应的复数所需要的迭代运算次数来决定的。 前面已经提到，这个迭代式是： Z := Z*Z + C 其中C是需要考察（计算的）复数点，对于每个绘图点来说，C就是是常数。而变量Z的初始值是0。 我们需要考察每次迭代之后Z的绝对值，看是否大于某常数（可以取为2)。   先分析绘图不同的绘图区域的特点。 对于远离原点的点，比如在以原点为圆心，半径为2的园之外的点，无须迭代，其绝对值就大于2了。不妨称这个区域为NM。 对于接近原点的点，大部分来说，无论迭代多少次，Z的绝对值都不会大于2。它们就是M集上的点了。这里也是最耗费时间的区域，如何设置迭代次数的上限，关系到整体的运算量。这个区域就直接称呼为M吧。 而对于介于以上二者之间的点，或者说是处于M集边缘的那些区域，则需要细细考察。如果一定次数之后绝对值超过2,它们仍不属于M集，但其颜色取决于对应的迭代次数，不同的迭代次数，将绘制出变化多端的图案。当然也有可能超过指定次数后绝对值仍不超过2,那就认为（但不一定全部是！）它们属于M集，仍绘制成黑色。我们将这个区域称为EM。 以上，区域NM和M是"平凡"的，因为结果很明显，要么属于M集，要么不属于M集（NM区域）。复杂的在于第三种情况，EM区域。这里，就算相邻"很近"的两个点，它们对应的运算结果就可能是很不同的（在最终绘图结果上，将表现为不同的像素颜色）。而这，也正是混沌现象的基本特征之一。很微小的初始条件的差别，或者说是对系统初始状态的很小的一点"扰动"，就可能使系统的后续发展最后朝完全不同的方向进行，得到截然不同的结果。这也就是所谓的"蝴蝶效应"。 但同时，任何两个点之间，即使相邻很近，它们之间也是没有"联系"的。每个点最后的结果，只与自己对应的复数值有关，不受任何其他点的影响。这样，我们就可以同时进行多个点的运算，而不用担心点与点之间有因果、先后关系。这是一个非常重要的特性，使我们可以在电脑上实现并行运算！   继续分析。 对于每个点，都需要进行多次的迭代运算。显然，在这个算法中，最多用到的，就是复数的平方、加法、取绝对值等，而这些运算在电脑上最终都可以归结为浮点数的加法和乘法。判定复数绝对值不超过2,等价于判定其平方不超过4。 所以，这是一个典型的运算密集型的算法，程序运行的大部分时间，将花费在这些浮点数运算上。所以，如何提高浮点数运算速度，将对优化算法有决定性作用。当然，从另一方面来说，如何用最少的计算量，达到同样的效果，也是需要考虑的地方。   总结一下，该算法主要有两大特点： 1, 典型的浮点数计算密集型算法。 出于这一点，提高浮点数处理速度，降低运算量，是优化算法的两大内容。 2, 非常适合于并行运算。 可以从CPU硬件、指令集、编译器、源代码等多方面考虑提高并行度。   以上两点，为优化该算法指明了方向。  

Continue reading about Mandelbrot集(4)：算法分析

Mandelbrot集(3)：我来画画看   （yanlb2000, 2007.04.28, yanlb2000.blogcn.com）   前面两篇文章，我简单介绍了Mandelbrot集，包括定义，如何画，大概是什么样子的。 这个M集，是混沌和分形的代表性图形，有着无穷精细复杂的结构。而且，用电脑来画出来，本身也不复杂。 所以，我很久前，就尝试用电脑编程来画M集的图像了。多年来，也算写过好几个版本了。这里就总结一下吧，也算是一次有意思的回顾。   1, VAX / VMS, VT340, Fortran 我最早画M集，还是十几年以前了。那时候还是学生，因为课程需要，在系里机房有上机的机会。但那时候的机器，是一台DEC VAX主机，带着20、30个终端！操作系统是VMS。大部分终端是纯字符模式的VT220,只有少数是带有图形功能的VT340(如果还没记错的话)。使用的编程语言是Fortran，而作图指令则是VT340终端扩展的ESC指令序列。 那时候上机机会本来就不多，更难“抢”到有限的图形终端。所以，那时候只是大致画出了M集的轮廓，很不完善。后来也马上不了了之了。   2, 286 / 486, DOS, TC 后来有机会用上了286电脑，上机机会也多了，能够好好画画这个神秘的图像了。那时候的操作系统还是DOS，编程用的开发语言是Turbo C 2.0(TC，非常经典的C语言开发环境)。使用的作图子系统是TC带的BGI库(Borland Graphics Interface)。 286电脑，用的CPU当然是80286,主频大约20MHz左右，而且是不带数学协处理器的，所有浮点运算都是软件实现的。所以，当时画一个M集的全图，640*480的分辨率，竟然需要2个小时左右。那时候，为了看一副图像，是很需要耐心的！ 再后来，用上了486电脑，当然是带有数学协处理器的。这时候，画一副图像的速度得以大大提高，大概需要一两分钟吧。   3,  Windows GDI, Delphi 终于进入Windows时代了。我选择了Borland的Delphi(使用Pascal语言)来画M集图像。作图指令方面，调用的是VCL控件的Canvas对象，实际上就是在Windows GDI上作图。电脑性能提高相对以前提高不少了，画一个M集图像大概需要十几秒几十秒左右。   4, AutoCAD, Lisp 用AutoCAD画M集是不常见的，至少我没听说过。早期的AutoCAD的二次开发语言一直是Lisp语言，称为AutoLisp。我用AutoCAD来画M集，一方面是想画个立体的M集图像，另一方面，是因为我对Lisp语言的喜爱。 Lisp语言是与Fortran差不多“古老”的语言，擅长表处理，号称人工智能语言。其实，我一直认为Fortran、C、Pascal、Basic等语言，其实都很相像。而Lisp就与众不同！ 诚然，排除硬件上的关系，这也是我所有相应软件中画M集最慢的。AutoLisp本来就是在AutoCAD下运行的解释型语言，运行是偏慢的。而且我画出来的，实际上都是3D对象。我使用了有高度的PLine对象，根据颜色的不同（也就是迭代次数的不同），使用不同的高度，类似于地图上的等高线和实际地形的关系。这样，画得慢是自然的。但画好之后，就不是简单的平面点阵图了，而是AutoCAD中真正的3D对象，可以渲染，可以变换多种角度和视角来观看。还是很有意思的。   5, Windows DirectX, VC++ 为了突破低效的Windows GDI，达到更快的作图速度，甚至是动画式的实时渲染，也为了练习使用DirectX，我最近使用MS VC++和DirectX SDK重写了M集作图的程序。

Continue reading about Mandelbrot集(3)：我来画画看

Mandelbrot集(2)：如何画，"长"什么样？ （yanlb2000, 2007.04.27, yanlb2000.blogcn.com） 上次，我简单介绍了Mandelbrot集的定义。这个定义和相应的运算，都是抽象的数学概念和过程。现在，我们要将M集这个数学概念形象地画出来，看看到底"长"什么模样。特别是，我们要好好看看M集的边界！ 我们要将复数平面上的图像画到我们的电脑屏幕上。这么做。首先定义一个复数平面区间，并将这个区间映射到屏幕作图区域上。这样，屏幕上每个作图点，就对应了一个复数。要决定这幅"画"怎么画，其实也很简单，只要定义每个作图点的颜色就行了。屏幕上的图画，不都是由彩色的点阵组成的么？ 如上，对于屏幕上的一个绘图点，我们可以算出其对应的复数值。然后，根据前文的迭代公式，进行运算。 如果在有限次数之内（比如256次），运算结果都不超过给定的常数值（比如2)，那么就认为该点属于M集。其对应颜色是黑色。 如果在上述有限次数之内，运算结果超出给定常数值，那么就认为该点不属于M集。给其赋予另外的颜色。而且，为了使画面颜色丰富，也为了形象地表示出该点是经过了多少次运算的，所以，根据运算的次数，来决定一个相应的颜色。 由上述算法，对作图区域上每个点都可以算出其颜色值了。接下来的任务，就是把这个算法编制成程序来实现了。 实际上，该算法并不复杂。所以，网上随便搜一下，就有大把的画Mandelbrot集的软件。 下面，还是实际看一下，这个M集是什么样子的吧。 先来看看Mandelbrot集的全貌。     是不是象个葫芦？其中，红线为x轴和y轴，两轴交叉处即为原点(0, 0)。图像最左端是(-2, 0)。 图中，靠近原点的一片区域是黑色。这表示，这些区域对应的任意一个复数，经过制定次数的迭代运算之后，结果序列的绝对值都没能超过指定的常数(一般还是选该常数R＝2)。所以，暂且认为它们是或可能是属于M集的。它们被描绘成黑色。 同时地，图像最外围是也黑色。这些区域表示，其中的点对应的任意复数C，其绝对值（到原点距离）早就超过2了。甚至不用做迭代，就可知它们不属于M集。它们也被描绘成黑色。 从最外围逐渐往原点靠近，情况变得微妙起来。有些点，可能其本身绝对值没有超过2。但是，经过若干次迭代运算之后，绝对值就超过2了。这些点显然也不属于M集。为了使图像丰富多彩，更为了表明该区域的点的特性，所以，根据它们需要迭代运算的次数n，来决定这个点的颜色。具体的颜色当然可以任意，但为了表达这种渐变的特性，所以让相邻大小的次数，其色彩值也是渐变的。（如上所述，颜色值是由运算次数来决定的。） 可以看到随着向原点的接近，需要迭代的次数也逐渐增加，色彩也相应渐渐变化。 所以，在这个图像中，最外围的黑色区域，以及所有有颜色的区域，都不属于M集。而每点不同的颜色，表明该点需要迭代运算的次数。只有中间那片黑色，才是M集的范围。 显然，这个区域的形状，出乎人们最初的预计。比如，我原来的第一判断是，这应该是个简单的园。然而，事实上，这个葫芦状的区域有点古怪，复杂得超出我们的预料。 通过计算和图形显示可发现，这个区域的边界非常地"不确定"。在黑色连接颜色区域(M集过渡到非M集)的地方，却也是最复杂、最精彩的地方。这里，没有一条明确的界线，将M集和非M集区分开来。对于这些区域的每个点，只能是实际地迭代运算一番，才能得出结论说，该点不属于M集，或者是可能属于M集。 所以，这个区域，显然是个开区间。 这个区域，其行为是"混沌"的。 这个边界，具备所谓的"分形"特性。 这个边界，是处处连续但处处不可微的。   好，我们不防将这个图像的某些区域放大，来一次M集上的旅行吧。也正好领略一下其复杂和精细的结构。 下面的第一幅图像，是前面那个大图的红色小框区域的放大。接下来的每个图像上的红色小框区域，代表了下一个将要放大的区域。  我们是不是放大了很多倍了？其实，只要我们愿意，这个“漫游”的游戏可以无限制地进行下去。无论放大多少倍，M集的边界永远是那么复杂。无限精细的结构，这是分形图形的特征之一。 当然，事实上，由于任何计算机的运算精度都是有限的，所以，现实中的任何程序其实是不能无限放大下去的。但这仅仅是物理实现上的限制。理论上，这个图像就是无限复杂的。    好，再来一次漫游。还是从原来那个全图开始，但放大另一个局部。   在某个层次，你是否发现，这个局部图像与原来的那个全图是相似的？其实，在不同放大比例的不同的局部，存在着无数个与全体类似的局部图形。 局部相似于整体，这是分形的另一个重要特征。    

Continue reading about Mandelbrot集(2)：如何画，"长"什么样？

Mandelbrot集(1)：定义   （yanlb2000, 2007.04.26, yanlb2000.blogcn.com）   Mandelbrot集(Mandelbrot Set，有人译为曼德尔布罗特集，我觉得挺拗口，还是直接写为Mandelbrote集,或简称M集)是分形(Fractal)和混沌(Chaos)领域的最著名最典型的话题。那个类似葫芦一样的分形图案，几乎成了混沌和分形的象征。 我很早就对分形感兴趣了。所以，这里，我也准备写在博客上，既是总结，也是给大家介绍一下。 分形和混沌是比较新兴的学科，其涉及的领域也很多。我这里，主要还是介绍一下最著名的Mandelbrt集。 对于复数平面上的一个点C，以及作为变量的复数Z（初始值为0），定义一个迭代运算过程： Z := Z^2 + C， 或者写为 Zn+1 = Zn ^ 2 + C 这个过程将产生一系列的 Z0,  Z1,  Z2, ...。 我们要考察的是Zn离原点的距离，即abs(Zn)。如果无论经过多少次运算，Zn 的绝对值都是有限的，就说C点属于Mandelbrot集。如果经过有限次数的运算之后，其绝对值可超过任意给定的值，（比如说R=2)，那么就说C点不属于Mandelbrot集。 如果你对复数平面不熟悉，那也可以理解为一个普通的二维平面，其上每个点都有x和y方向的坐标，这个点就代表了一个复数 C = x + y*i。x和y就是这个复数的实部和虚部。 显然，Mandelbrot集是个抽象的数学概念。经过一些试验计算，我们可以发现，那些离原点比较远的点，肯定是不属于M集的，不断地平方再加自己，可以迅速超过任何给定数值。而那些离原点很近的点，则无论如何都不可能变大。因此，那些紧靠原点的点，显然属于M集。而远离原点的，肯定不属于M集。 现在问题来了，在M集和非M集之间，分界线在哪里？有没有一个明确的界线，将M集和非M集划分开来？或者说，能明确地给出这个界线的定义（解析式）吗？又或者说，有没有一个明确而简便的方法，使我们可以对任意一个复数C，给出其是否属于M集的结论？ 很遗憾，这个问题远非相像的那么简单。没有这样的简单方法，没有这样明确的界线。 对于一个给定复数C，除了实际验算，无法给出明确的答案。 而就算根据前面的定义实际验算了，结论也是复杂的。 如果经过一定次数的迭代运算，Z的绝对值超出了设定的常数R。那么很好，这个C不属于M集。 但也有可能，就算经过1000次运算，其绝对值还是很小。那么，就可以说C属于M集了吗？不一定！有可能，在接下来的1001次或更以后，就可以发现Z超出了R！ 按理，上述迭代过程是个非常确定性的过程，而且很简单。所以，对于任意一个给定的C，其是否属于M集，应该是确定的。但实际上，对于某些C值，我们有可能无论经过多少次迭代，都无法给出结论，而我们又不能说，这个C就不属于M集了，因为说不定增加迭代次数，就发现超出R了。结果不明，这就是“混沌”了。

Continue reading about Mandelbrot集(1)：定义