谢尔宾斯基地毯与圆的交集

wayne · 发表于 2017-9-28 16:22:41

我把我的代码贴出来，用到了 C++ + cln(https://www.ginac.de/CLN/)
圆上的格子用线性队列存储。元素个数以$8/3$的倍率指数递增，所以此处是性能瓶颈（内存和cpu），时空的复杂度均为 $O((8/3)^n),n $为迭代次数。

#include <iostream>
#include <queue>
#include <cln/rational_io.h>
using namespace std;
typedef cln::cl_RA T;
//typedef unsigned long long T;
typedef struct counter{T half;T one;} counter;
typedef struct Point2d{cln::cl_RA x;cln::cl_RA y;} Point2d;
typedef struct Context {
counter outer;
std::queue<Point2d> on;
counter inner;
cln::cl_RA sideLen;
} Context;
void walkMengerCircle(Context &c){
static int kernel[8][2] = {{-1, -1},{-1, 0},{-1, 1},{0, -1},{0, 1},{1, -1},{1, 0},{1, 1}};
cln::cl_RA newslength = c.sideLen/3;
c.sideLen = newslength;
///////////////////////////
counter out;
out.half = c.outer.half*2;
out.one = c.outer.one*8 + c.outer.half*3;
///////////////////////////
counter in;
in.half = out.half*3-1;
in.one = c.inner.one*8 + c.inner.half*3;
///////////////////////////
for(size_t j=0,s = c.on.size();j<s;++j){
Point2d v = c.on.front();
c.on.pop();
for(size_t i=0;i<8;++i){
Point2d center = {v.x+ kernel[i][0]*newslength,v.y+ kernel[i][1]*newslength};
Point2d lower = {center.x - newslength/2,center.y - newslength/2};
Point2d upper = {center.x + newslength/2,center.y + newslength/2};
cl_signean s1 = cln::compare("1/4",lower.x*lower.x + lower.y*lower.y);
cl_signean s2 = cln::compare("1/4",upper.x*upper.x + upper.y*upper.y);
//assuming that s2<s1 is always valid.
if(s2<0&&s1>0){
c.on.push(center);
} else if(s1<0){
out.one ++;
} else {
in.one ++;
}
}
}
///////////////////////////
c.on.push({"1/2"-newslength/2,newslength});
in.one +=2;
///////////////////////////
c.inner = in;
c.outer = out;
}
int main()
{
/*
{{{0, 8, 0}, 8},
{{4, 28, 32}, 64},
{{104, 76, 332}, 512},
{{984, 204, 2908}, 4096},
{{8288, 580, 23900}, 32768},
{{67744, 1556, 192844}, 262144},
{{545904, 4180, 1547068}, 2097152},
{{4377944, 11204, 12388068}, 16777216},
{{35052688, 29724, 99135316}, 134217728},
{{280499768, 79276, 793162780}, 1073741824},
{{2244206376, 212076, 6345516140}, 8589934592},
{{17954209288, 565692, 50764701756}, 68719476736},
{{143635171632, 1509332, 406119132924}, 549755813888},
{{1149085383304, 4026028, 3248957101772}, 4398046511104},
{{9192693786392, 10740796, 25991667561644}, 35184372088832},
{{73541578891856, 28646804, 207933369171996}, 281474976710656},
{{588332707461496, 76396620, 1663467029827132}, 2251799813685248},
{{4706661863240488, 203728972, 13307736442512524}, 18014398509481984},
{{37653295449008288, 543283204, 106461892083564380}, 144115188075855872},
{{301226365040331144, 1448779164, 851695138117736668}, 1152921504606846976}}
*/
// {{1,0},{{"4/9", "1/9"}, {"4/9", "2/9"}, {"4/9", "1/3"}},{2,3},"1/9"}
// {{2^(n-2),*},******************************************,{3*2^(n-2)-1,*},"1/3^n"}
Context context;
context.outer = {1,0};
context.on.push({"4/9", "1/9"});
context.on.push({"4/9", "2/9"});
context.on.push({"4/9", "3/9"});
context.inner = {2,3};
context.sideLen = "1/9";
for(int depth=3;depth<=15;++depth){
walkMengerCircle(context);
cout <<depth<<":{"<<8*context.outer.one+4*context.outer.half<<","
<<8*context.on.size()+4<<","
<< 8*context.inner.one+ 4*context.inner.half<<"},"<<context.sideLen;
// cout <<depth<<":{"<<context.outer.half<<","<<context.outer.one<<"},"
// <<context.on.size()<<",{"
// <<context.inner.half<<","<<context.inner.one<<"},"<<context.sideLen;
cout <<endl;
}
return 0;
}

复制代码

征得KeyTo9_Fans的同意后，也把他的代码贴出来：

#include<cstdio>
int i,l=20;
unsigned __int64 ou[96],in[96],s,r;
void f(unsigned __int64 a,unsigned __int64 b,unsigned __int64 c,unsigned __int64 d,unsigned __int64 r,int h)
{
if(b*b+d*d<=r)in[h]+=2;
else
if(a*a+c*c>=r)ou[h]+=2;
else
if(h++<l)
r*=9,
f(a+a+a,a+a+b,c+c+c,c+c+d,r,h),
f(a+a+b,a+b+b,c+c+c,c+c+d,r,h),
f(a+b+b,b+b+b,c+c+c,c+c+d,r,h),
f(a+a+a,a+a+b,c+c+d,c+d+d,r,h),
f(a+b+b,b+b+b,c+c+d,c+d+d,r,h),
f(a+a+a,a+a+b,c+d+d,d+d+d,r,h),
f(a+a+b,a+b+b,c+d+d,d+d+d,r,h),
f(a+b+b,b+b+b,c+d+d,d+d+d,r,h);
}
int main()
{
ou[2]=1,in[2]=3;
f(3,5,7,9,81,2),f(5,7,7,9,81,2);
for(r=9,s=16,i=2;i<=l;i++,s<<=3,r*=3)
{
in[i]+=5;
f(1,3,r-2,r,r*r,i);
if(i>2)in[i]+=in[i-1]<<3,ou[i]+=ou[i-1]<<3;
printf("{%I64u,%I64u,%I64u},%I64u,{%.16lf,%.16lf,%.16lf}\n",ou[i]<<2,s-in[i]-ou[i]<<2,in[i]<<2,s<<2,1.0*ou[i]/s,1.0*(s-in[i]-ou[i])/s,1.0*in[i]/s);
}
return 0;
}

复制代码

他只取$1/8$的圆，直接计算第$2$层到第$l$层内外部的完整正方形个数，没有保存中间结果，全程采用$64bit$的整数运算，时间复杂度为$O((8/3)^n)$，空间复杂度为$O(n)$。
将上述代码稍作修改，进行更多层迭代：

#include<cstdio>
int i,l=32;
unsigned __int64 ou[96],in[96],s,r,z=1ULL<<63;
void f(unsigned __int64 a,unsigned __int64 b,unsigned __int64 c,unsigned __int64 d,unsigned __int64 r,int h)
{
if(r-b*b-d*d<z)in[h]+=2;
else
if(a*a+c*c-r<z)ou[h]+=2;
else
if(h++<l)
r*=9,
f(a+a+a,a+a+b,c+c+c,c+c+d,r,h),
f(a+a+b,a+b+b,c+c+c,c+c+d,r,h),
f(a+b+b,b+b+b,c+c+c,c+c+d,r,h),
f(a+a+a,a+a+b,c+c+d,c+d+d,r,h),
f(a+b+b,b+b+b,c+c+d,c+d+d,r,h),
f(a+a+a,a+a+b,c+d+d,d+d+d,r,h),
f(a+a+b,a+b+b,c+d+d,d+d+d,r,h),
f(a+b+b,b+b+b,c+d+d,d+d+d,r,h);
}
int main()
{
ou[2]=1,in[2]=3;
f(3,5,7,9,81,2),puts("*"),f(5,7,7,9,81,2),puts("**");
for(r=9,s=16,i=2;i<=l;i++,s<<=3,r*=3)
{
in[i]+=5;
f(1,3,r-2,r,r*r,i);
printf("%I64u %I64u\n",in[i],ou[i]);
}
return 0;
}

复制代码

虽然坐标的平方和会超出$64bit$整数范围，但是与圆半径的平方进行比较的时候，两者的差异不会超过$2^63$，因此只看平方运算结果的最后$64bit$足矣

于是跑32层迭代也是妥妥的，根本无需进行更高精度的运算，只是需要

18天罢了。

wayne · 发表于 2017-9-29 17:13:00

仔细琢磨起来，这个维度怎么来度量有点名堂，突然间觉得我们在这编程失去了参考意义。
首先比较平凡的结论是，经过$n$次迭代后，小格子的边长是 $1/3^n$，黑白格子沿着$x$轴，$y$轴均匀分布，其个数是 $3^n$, 二维平面的面积也是被黑白格子均匀的填充着，个数是$9^n$
但如果计算黑格子，非平凡的问题就来了：

在二维平面内的$9^n$个黑白格子的集体中，
1）$x$轴，$y$轴穿过的黑格子数目是 $2^n$, 沿着对角线直线穿过的黑格子数目也是 $2^n$,其他直线不是。
2）$x^2+y^2=1$的圆周非均匀分布的穿过的黑格子个数是 $(8/3)^n$数量级
--------------------------------------------------------------------------------------------------
3）整个二维的平面非各向同性的包含 $8^n$个格子。
4）在$x^2+y^2<1$的圆内且在大正方形区域内非均匀的包含的黑格子个数是 $8^n$数量级
5）在$x^2+y^2>1$的圆外且在大正方形区域内非均匀的包含的黑格子个数是 $8^n$数量级

==========================================================
6）当n趋于无穷的时候，地毯的二维空间的面积为0
7）当n趋于无穷的时候，地毯的一维空间的点的个数为0
==========================================================

问题就出在我们正在用二维的格子，或者一维的点来度量 $log(8)/log(3)$维的空间面积。是否本身就不恰当？是否应该选择一个 $log(8)/log(3)$维的基本元素来组成这个空间的“面积”才比较合适？连续性在这里到底是啥玩意？

楼主Fans在这里将二维空间的圆跟 $log(8)/log(3)$维的地毯混在一块到底该作何理解？

zeroieme · 发表于 2017-9-29 19:06:01

wayne 发表于 2017-9-29 17:13
仔细琢磨起来，这个维度怎么来度量有点名堂，突然间觉得我们在这编程失去了参考意义。
首先比较平凡的结论 ...

这个白格子是白格子。黑格子却不是黑格子，仍然是空洞的。姑且算“灰格子”吧。

这些灰格子不是二维还是$Log_3(8)$维几何体。所以前面的编程还是意义的，算是数值上进行勒贝格积分吧。

zeroieme · 发表于 2017-9-29 22:03:48

点就是0维的点。
平面是空洞的，但有连续线段，就是每个谢尔宾斯基地毯的纵横 1/3,2/3分割线。

KeyTo9_Fans · 发表于 2017-10-5 00:09:30

wayne 发表于 2017-9-28 13:29
......
修改Fans的代码，……，2天12小时27分完成30次迭代，预计在2017年10月23日完成32次迭代。

等不及$32$次迭代了，先处理一下前$30$项数据：

{0,8,0}
{4,28,32}
{104,76,332}
{984,204,2908}
{8288,580,23900}
{67744,1556,192844}
{545904,4180,1547068}
{4377944,11204,12388068}
{35052688,29724,99135316}
{280499768,79276,793162780}
{2244206376,212076,6345516140}
{17954209288,565692,50764701756}
{143635171632,1509332,406119132924}
{1149085383304,4026028,3248957101772}
{9192693786392,10740796,25991667561644}
{73541578891856,28646804,207933369171996}
{588332707461496,76396620,1663467029827132}
{4706661863240488,203728972,13307736442512524}
{37653295449008288,543283204,106461892083564380}
{301226365040331144,1448779164,851695138117736668}
{2409810924185402784,3863345612,6813561108806027412}
{19278487403784147304,10302538780,54508488880751520380}
{154227899257743649824,27473690092,436067911073488311796}
{1233823194135208730288,73263231116,3488543288661173252292}
{9870585553277031899992,195369181668,27908346309484760627908}
{78964684426737227956608,520985280228,223266770476399080439708}
{631717475415287101141792,1389296277316,1786134163812581951993244}
{5053739803326001562892000,3704793953044,14289073310504360438453772}
{40429918426617891895941808,9879454623900,114312586484044763011824820}
{323439347412969480294205216,26345219429236,914500691872384449385489772}

复制代码

首先按照hujunhua在$27$楼求面积上下界的代数平均数的方法，求得以下面积序列：

0.5
0.71875
0.72265625
0.73486328125
0.73822021484375
0.73860931396484375
0.73869609832763671875
0.73872029781341552734375
0.73872639238834381103515625
0.73872731812298297882080078125
0.73872764804400503635406494140625
0.73872775249765254557132720947266
0.73872777209544437937438488006592
0.73872777529368249815888702869415
0.73872777568459468966466374695301
0.7387277758233707913859689142555
0.73872777585101889741281411261298
0.73872777586064863886150533289765
0.73872777586196634869164512338102
0.73872777586238128934292834770758
0.73872777586245612292843720769753
0.73872777586247222368306432349616
0.73872777586247703259742587895975
0.73872777586247779770407592540328
0.73872777586247796382687745843604
0.73872777586247800131432562497993
0.7387277758624780077222150705618
0.73872777586247800949936495034918
0.73872777586247800985928169095094
0.73872777586247800992560980652542

复制代码

然后对上述面积序列的相邻两项作差，得到：

0.5
0.21875
0.00390625
0.01220703125
0.00335693359375
0.00038909912109375
0.00008678436279296875
0.00002419948577880859375
0.00000609457492828369140625
0.00000092573463916778564453125
0.00000032992102205753326416015625
0.00000010445364750921726226806641
0.00000001959779183380305767059326
0.00000000319823811878450214862823
0.00000000039091219150577671825886
0.00000000013877610172130516730249
0.00000000002764810602684519835748
0.00000000000962974144869122028467
0.00000000000131770983013979048337
0.00000000000041494065128322432656
0.00000000000007483358550885998995
0.00000000000001610075462711579863
0.00000000000000480891436155546359
0.00000000000000076510665004644353
0.00000000000000016612280153303276
0.00000000000000003748744816654389
0.00000000000000000640788944558187
0.00000000000000000177714987978738
0.00000000000000000035991674060176
0.00000000000000000006632811557448

复制代码

再对上述序列相邻两项作商，得到：

1 2.2857142857142857142857142857143
2 56
3 0.32
4 3.6363636363636363636363636363636
5 8.6274509803921568627450980392157
6 4.4835164835164835164835164835165
7 3.5862068965517241379310344827586
8 3.9706601466992665036674816625917
9 6.5835010060362173038229376257545
10 2.8059280169371912491178546224418
11 3.1585399832822513234884367771037
12 5.3298682012251717096714318820994
13 6.1276837764822977392293561428623
14 8.1814744801512287334592972660495
15 2.8168552557472735651016350792179
16 5.0193710045295383725782481751307
17 2.871116132676942942459197582619
18 7.3079377784330865393824639461471
19 3.1756585575906052333539618816487
20 5.5448452517900141689032828041249
21 4.6478309397269083352206042808509
22 3.3481059167599611362212550782029
23 6.285286320885529757449321963342
24 4.6056690772476863673157818475219
25 4.4314246409893269354527425736722
26 5.8502020805603697358084394450449
27 3.6057113237677602868585857216454
28 4.9376694088085151486851873689654
29 5.4263073431900628788285733338051

复制代码

然后根据上述序列作散点图，得到：

最后根据散点的变化趋势预测后续数据，并逆推得到所求面积的极限为（$20$个有效数字）：

$0.738727775862478009946(3)$

最后的$(3)$表示最后两位数字$46$的

$68%$置信区间是$46\pm03$，
$95%$置信区间是$46\pm06$，
$99.7%$置信区间是$46\pm09$，
$99.99%$置信区间是$46\pm12$，
$99.9999%$置信区间是$46\pm15$，
……

zeroieme · 发表于 2017-10-12 19:01:44

抛砖引玉说一下我的极限推导思路。就是仿照古典积分从矩形逼近曲线的过程。
把地毯切割，作出累积曲线。切割地毯1维宽度每缩小1/3，累积的谢尔宾斯基地毯缩小为1/8。

未切割的累积曲线就是

未命名a-1.png

切割成3段有这2种曲线
未命名b-1.png

切割成9段有这几种曲线
未命名c-1.png

每种曲线，是切割宽度3^(-n)和高度y坐标的函数。结合终止累积位置x的数值就得到所需区域的谢尔宾斯基地毯近似数值。
当n->∞就能得到精确值。

Ickiverar · 发表于 2017-10-17 16:14:28

使用以下思路可以避免 int64 之间的乘法，而且大大降低了操作数的大小（操作数最大约3^n，对n=32，远远小于maxint64），且消灭了对 r 的传参。
思路：每次迭代，正方形坐标就会从 x，y，变成 3x+2，3y+2，3x，3y，3x-2，3y-2的一些组合。而我们关心的只有 fh=(x+1)^2-(y+1)^2-r^2 和 fl=(x-1)^2+(y-1)^2-r^2 这两个值，如果 fh<=0，正方形在圆内，如果 fl>=0，正方形在圆外。
那么我们可以直接建立 fh 和 fl 的迭代关系，而不需要每次都由平方差计算 fh 和 fl。
一般来说，新的 fh 是 9倍原来的fh 加上x和y的一些线性组合。fl 类似。

#include<stdio.h>
#include<stdint.h>
uint64_t in[64]={0},out[64]={0};
int maxl=24;
//fh:=(x+1)^2-(y+1)^2-r^2
//fl:=(x-1)^2-(y-1)^2-r^2
void ffull(int l,int64_t x,int64_t y,int64_t fh,int64_t fl){
if(++l==maxl){
in[l]+=32;
out[l]+=32;
return;
};
x*=3;
y*=3;
fh*=9;
fl*=9;
int64_t nfh,nfl;
#define proc(nx,ny,dfh,dfl) \
do{ \
nfh=fh-4*(dfh),nfl=fl+4*(dfl); \
if(nfl>=0)out[l]+=8; \
else if(nfh<=0)in[l]+=8; \
else ffull(l,nx,ny,nfh,nfl); \
}while(0)
proc(x+2,y+2,0,2*x+2*y-4);
proc(x+2,y,y+2,2*x+y-4);
proc(x+2,y-2,2*y+2,2*x-2);
proc(x,y+2,x+2,x+2*y-4);
proc(x,y-2,x+2*y+4,x-2);
proc(x-2,y+2,2*x+2,2*y-2);
proc(x-2,y,2*x+y+4,y-2);
proc(x-2,y-2,2*x+2*y+4,0);
#undef proc
}
void fhalf(int l,int64_t x,int64_t fh,int64_t fl){
if(++l==maxl){
in[l]+=32;
return;
};
x*=3;
fh*=9;
fl=9*fl+4*(2*x-4);
in[l]+=20;
ffull(l,x+2,2,fh,fl);
fhalf(l,x+2,fh-8,fl);
}
int main(){
in[2]=32,out[2]=4;
ffull(2,8,6,49,-7);
printf("*");
ffull(2,8,4,25,-23);
printf("*");
ffull(2,8,2,9,-31);
printf("*");
fhalf(2,8,1,-31);
double inres=0,outres=0;
for(int i=2;i<=maxl;++i){
inres=8*inres+in[i];
outres=8*outres+out[i];
printf("%18.16lf\n",inres/(outres+inres));
}
return 0;
}

复制代码

在我的电脑上，单线程24层不到一分钟。
我再试试并行化。

Ickiverar · 发表于 2017-10-19 15:58:56

32层计算完毕。
计算开始于 2017年10月17日22:30:49.46
计算结束于 2017年10月19日13:04:02.45
耗时 38:33:12.99
结果：

{0, 8, 0}
{4, 28, 32}
{104, 76, 332}
{984, 204, 2908}
{8288, 580, 23900}
{67744, 1556, 192844}
{545904, 4180, 1547068}
{4377944, 11204, 12388068}
{35052688, 29724, 99135316}
{280499768, 79276, 793162780}
{2244206376, 212076, 6345516140}
{17954209288, 565692, 50764701756}
{143635171632, 1509332, 406119132924}
{1149085383304, 4026028, 3248957101772}
{9192693786392, 10740796, 25991667561644}
{73541578891856, 28646804, 207933369171996}
{588332707461496, 76396620, 1663467029827132}
{4706661863240488, 203728972, 13307736442512524}
{37653295449008288, 543283204, 106461892083564380}
{301226365040331144, 1448779164, 851695138117736668}
{2409810924185402784, 3863345612, 6813561108806027412}
{19278487403784147304, 10302538780, 54508488880751520380}
{154227899257743649824, 27473690092, 436067911073488311796}
{1233823194135208730288, 73263231116, 3488543288661173252292}
{9870585553277031899992, 195369181668, 27908346309484760627908}
{78964684426737227956608, 520985280228, 223266770476399080439708}
{631717475415287101141792, 1389296277316, 1786134163812581951993244}
{5053739803326001562892000, 3704793953044, 14289073310504360438453772}
{40429918426617891895941808, 9879454623900, 114312586484044763011824820}
{323439347412969480294205216, 26345219429236, 914500691872384449385489772}
{2587514779303826096159159104, 70253935029820, 7316005534979145849098804868}
{20700118234430796112810254760, 187343834543092, 58528044279833354136899152484}

复制代码

KeyTo9_Fans · 发表于 2017-10-19 21:46:04

耗时 38:33:12.99 也就是 1.6064天。

竟然真的是1.6天

#####

根据楼上给出的$31$层和$32$层的结果，求得面积中项：

$31:\ 0.73872777586247800993618025568671$
$32:\ 0.7387277758624780099397883894226$

相邻两项作差：

$31:\ 0.00000000000000000001057044916129$
$32:\ 0.00000000000000000000360813373589$

相邻两项作商：

$30:\ 6.2748625495858708913211618483576$
$31:\ 2.9296167866911510306989426440087$

散点图：

从图中可以看到，

最新的两个数据点并不在我画的两条曲线框定的范围里，

说明我在$35$楼预测的极限与真实的极限偏差较大，

真实的极限不在我给出的$68%$置信区间里，

而是在我给出的$95%$置信区间的边缘，

差一点就超出$95%$置信区间的范围了。

#####

根据最新的$2$个数据重新预测极限，结果为（$21$个有效数字）：

$0.7387277758624780099409(3)$

KeyTo9_Fans · 发表于 2017-11-3 08:38:55

哇！这个贴子变成精品贴了？好开心呀~

#####

拿到$18$项数据之后，我就把这个数列提交到【在线整数数列百科大全】网站上去了。

可惜自提交至今的一个多月里，虽然做了多次修改，但还是没有通过审核。

因此这个数列还没有正式成为【在线整数数列百科大全】里的词条。

#####

而这两个数列提交上去之后则是秒过，差别也太大了。

账号		自动登录	找回密码
密码			欢迎注册

[原创] 谢尔宾斯基地毯与圆的交集

点评

评分

点评

点评

评分

点评