今天又是因为可以用py而高兴的一天。

继续咱的sierpinski地毯计划。

二，随机算法

在二十年前，磁盘容量以MB还是KB计的时候，分形解决计图的问题确实有很大的优势。存至多十来个数就好了。我要在二十年前，敢用这种几百×几百的矩阵直接存数据，那肯定被计算机打死，呃，被老板打死。。。电脑怎么烧掉的都不知道。。。

这种意义上说，随机算法才是传统算法。。。之前那个数学硬解的玩意，我是不信十五年前的计算机能支持的。。。（说的好像那时候我有电脑用似的。。。呃好像还真有。。。）

说到分形，不得不提迭代函数系。

康托三分集和谢尔宾斯基垫片，谢尔宾斯基地毯都是由若干个压缩映射确定的。

[1]沙震，阮火军：分形与拟合[M]. 杭州：浙江大学出版社，2005.3

啥是完备度量空间(X;d)上的双曲迭代函数系呢？每个映射都是压缩映射。

我们今天用的是带概率的迭代函数系。不是有几个压缩映射么？给一个加权，通过抛硬币选择走那条路。视为一次迭代。

如果在随机迭代出来的点序列中，去掉前若干项（万一我从中间开始的，岂不是凉凉月色为你。。。），那对于充分长的一段点列，它与K的hausdorff距离小于e的概率大于1-e。

换句话说，如果我迭代了足够多次，那么我应该能跑出来，至于跑不出来的情形，我觉得你应该换个随机数生成法。。。不对，不如买彩票。。。买好彩票就能换台电脑了，反复迭代就可以买一台超算。

唔，咱应该不用解释啥是完备，啥是度量空间，啥是压缩映射吧？

至于想问hausdorff距离的同学，来蹭课吧。

接着我们进入算法环节。

这个sierpinski地毯呢，有八个压缩映射确定。

正常的流程是写九个函数。

def f1():

def f2():

def frandom():

然而机智的我早就看穿一切。

为啥要写8+1个？我写9+1个不就好了？中间第5个是空函数嘛。

我先写了一个映射，将random的范围(0,1)映射到(0,4/9)U(5/9,1)。

这样我×3取整就能得到是第几排了，去掉整数部分，再×3取整就知道是第几列了。

我们刚刚是不是取到了两个0，1，2的数？

嘿嘿嘿！（发出了政委的笑声）

def next(x,y):

    t=random.random()

    if t<=0.5:

        t=t/4.5*4

    elif t>0.5:

        t=(t-0.5)/4.5*4+5/9

    i=int(t*3)

    j=int((t-int(t*3)/3)*9)

    xx=x/3+i/3

    yy=y/3+j/3

    return xx,yy

主函数是这样的：

import cv2

import numpy as np

import ct

import random

times=7

k=3**times

img=ct.cut(times)

x=random.random()

y=random.random()

a=np.zeros((k,k),dtype=int)

trytime=20000000

for i in range(trytime):

    [x,y]=ct.next(x,y)

    if i>300:

        a[int(k*x)][int(k*y)]=1

img=ct.print(a,times,img)

cv2.imshow(“out”,img)

name=str(“bigrandom”+str(int(trytime))+”.jpg”)

cv2.imwrite(name,img)

我决定再强调一遍前xxx项要扔掉！

我们要写什么样的算法？要写禁得起比较的算法，所以这矩阵多大呢？

3的5次方*3的5次方。

上面分别是迭代80w,10w,2w,5k,1k次后的结果。

下面的分别是迭代2000w,300w,65w,10w次的结果。

为啥是这些数，我手里的图片查看器会反走样，要把图片放大才会看到点。。。在点很稀疏的时候，这里，我不知道，如果你看到一坨黑的，那有可能，就是被反走样了。也可能就是因为点少的可怜。。。

至于为啥今天和昨天的颜色相反。这是因为昨天我是用np.zeros然后用1标记白的（去掉的）；今天我还这么来，用np.zeros初始化，用1标记黑的。你想一致，很简单，你用np.ones初始化，然后用0标记黑的。

“可以，但没必要。”

下次不出意外的话是sierpinski垫片。

能用py就很开心。我写两个用了十分钟。

如果是c++我可能环境都没配置好。不对，我有可能就把这事晾一边了。

昨天晚上在准备羽毛球裁判的讲座。总觉得有啥忘了的。。。（划掉）

我自闭了。这图片咋调大小。。。woc,我的公式也被吃了。