数学研发论坛

 找回密码
 欢迎注册
楼主: KeyTo9_Fans

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

[复制链接]
发表于 2017-9-28 16:22:41 | 显示全部楼层
我把我的代码贴出来,用到了 C++  + cln(https://www.ginac.de/CLN/)
圆上的格子用线性队列存储。元素个数以$8/3$的倍率指数递增,所以此处是性能瓶颈(内存和cpu),时空的复杂度均为 $O((8/3)^n),n $为迭代次数。
  1. #include <iostream>
  2. #include <queue>
  3. #include <cln/rational_io.h>
  4. using namespace std;

  5. typedef cln::cl_RA T;
  6. //typedef unsigned long long T;
  7. typedef struct counter{T half;T one;} counter;
  8. typedef struct Point2d{cln::cl_RA x;cln::cl_RA y;} Point2d;
  9. typedef struct Context {
  10.     counter outer;
  11.     std::queue<Point2d> on;
  12.     counter inner;
  13.     cln::cl_RA sideLen;
  14. } Context;

  15. void walkMengerCircle(Context &c){
  16.     static int kernel[8][2] = {{-1, -1},{-1, 0},{-1, 1},{0, -1},{0, 1},{1, -1},{1, 0},{1, 1}};

  17.     cln::cl_RA newslength = c.sideLen/3;
  18.     c.sideLen = newslength;
  19.     ///////////////////////////
  20.     counter out;
  21.     out.half = c.outer.half*2;
  22.     out.one  = c.outer.one*8 + c.outer.half*3;
  23.     ///////////////////////////
  24.     counter in;
  25.     in.half = out.half*3-1;
  26.     in.one  = c.inner.one*8 + c.inner.half*3;
  27.     ///////////////////////////
  28.     for(size_t j=0,s = c.on.size();j<s;++j){
  29.         Point2d v = c.on.front();
  30.         c.on.pop();
  31.         for(size_t i=0;i<8;++i){
  32.             Point2d center = {v.x+ kernel[i][0]*newslength,v.y+ kernel[i][1]*newslength};
  33.             Point2d lower = {center.x - newslength/2,center.y - newslength/2};
  34.             Point2d upper = {center.x + newslength/2,center.y + newslength/2};
  35.             cl_signean s1 = cln::compare("1/4",lower.x*lower.x + lower.y*lower.y);
  36.             cl_signean s2 = cln::compare("1/4",upper.x*upper.x + upper.y*upper.y);
  37.             //assuming that s2<s1 is always valid.
  38.             if(s2<0&&s1>0){
  39.                 c.on.push(center);
  40.             } else if(s1<0){
  41.                 out.one ++;
  42.             } else {
  43.                 in.one ++;
  44.             }
  45.         }
  46.     }
  47.     ///////////////////////////
  48.     c.on.push({"1/2"-newslength/2,newslength});
  49.     in.one +=2;
  50.     ///////////////////////////
  51.     c.inner = in;
  52.     c.outer = out;
  53. }

  54. int main()
  55. {
  56.     /*
  57.      {{{0, 8, 0}, 8},
  58.       {{4, 28, 32}, 64},
  59.       {{104, 76, 332}, 512},
  60.       {{984, 204, 2908}, 4096},
  61.       {{8288, 580, 23900}, 32768},
  62.       {{67744, 1556, 192844}, 262144},
  63.       {{545904, 4180, 1547068}, 2097152},
  64.       {{4377944, 11204, 12388068}, 16777216},
  65.       {{35052688, 29724, 99135316}, 134217728},
  66.       {{280499768, 79276, 793162780}, 1073741824},
  67.       {{2244206376, 212076, 6345516140}, 8589934592},
  68.       {{17954209288, 565692, 50764701756}, 68719476736},
  69.       {{143635171632, 1509332, 406119132924}, 549755813888},
  70.       {{1149085383304, 4026028, 3248957101772}, 4398046511104},
  71.       {{9192693786392, 10740796, 25991667561644}, 35184372088832},
  72.       {{73541578891856, 28646804, 207933369171996}, 281474976710656},
  73.       {{588332707461496, 76396620, 1663467029827132}, 2251799813685248},
  74.       {{4706661863240488, 203728972, 13307736442512524}, 18014398509481984},
  75.       {{37653295449008288, 543283204, 106461892083564380}, 144115188075855872},
  76.       {{301226365040331144, 1448779164, 851695138117736668}, 1152921504606846976}}
  77.     */

  78.     // {{1,0},{{"4/9", "1/9"}, {"4/9", "2/9"}, {"4/9", "1/3"}},{2,3},"1/9"}
  79.     // {{2^(n-2),*},******************************************,{3*2^(n-2)-1,*},"1/3^n"}

  80.     Context context;
  81.     context.outer = {1,0};
  82.     context.on.push({"4/9", "1/9"});
  83.     context.on.push({"4/9", "2/9"});
  84.     context.on.push({"4/9", "3/9"});
  85.     context.inner = {2,3};
  86.     context.sideLen = "1/9";

  87.     for(int depth=3;depth<=15;++depth){
  88.         walkMengerCircle(context);
  89.         cout <<depth<<":{"<<8*context.outer.one+4*context.outer.half<<","
  90.              <<8*context.on.size()+4<<","
  91.              << 8*context.inner.one+ 4*context.inner.half<<"},"<<context.sideLen;
  92. //        cout <<depth<<":{"<<context.outer.half<<","<<context.outer.one<<"},"
  93. //             <<context.on.size()<<",{"
  94. //             <<context.inner.half<<","<<context.inner.one<<"},"<<context.sideLen;
  95.         cout <<endl;
  96.     }

  97.     return 0;
  98. }
复制代码

征得KeyTo9_Fans的同意后,也把他的代码贴出来:
  1. #include<cstdio>

  2. int i,l=20;
  3. unsigned __int64 ou[96],in[96],s,r;

  4. void f(unsigned __int64 a,unsigned __int64 b,unsigned __int64 c,unsigned __int64 d,unsigned __int64 r,int h)
  5. {
  6.         if(b*b+d*d<=r)in[h]+=2;
  7.         else
  8.                 if(a*a+c*c>=r)ou[h]+=2;
  9.                 else
  10.                         if(h++<l)
  11.                                 r*=9,
  12.                                 f(a+a+a,a+a+b,c+c+c,c+c+d,r,h),
  13.                                 f(a+a+b,a+b+b,c+c+c,c+c+d,r,h),
  14.                                 f(a+b+b,b+b+b,c+c+c,c+c+d,r,h),
  15.                                 f(a+a+a,a+a+b,c+c+d,c+d+d,r,h),
  16.                                 f(a+b+b,b+b+b,c+c+d,c+d+d,r,h),
  17.                                 f(a+a+a,a+a+b,c+d+d,d+d+d,r,h),
  18.                                 f(a+a+b,a+b+b,c+d+d,d+d+d,r,h),
  19.                                 f(a+b+b,b+b+b,c+d+d,d+d+d,r,h);
  20. }

  21. int main()
  22. {
  23.         ou[2]=1,in[2]=3;
  24.         f(3,5,7,9,81,2),f(5,7,7,9,81,2);
  25.         for(r=9,s=16,i=2;i<=l;i++,s<<=3,r*=3)
  26.         {
  27.                 in[i]+=5;
  28.                 f(1,3,r-2,r,r*r,i);
  29.                 if(i>2)in[i]+=in[i-1]<<3,ou[i]+=ou[i-1]<<3;
  30.                 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);
  31.         }
  32.         return 0;
  33. }
复制代码

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

  2. int i,l=32;
  3. unsigned __int64 ou[96],in[96],s,r,z=1ULL<<63;

  4. void f(unsigned __int64 a,unsigned __int64 b,unsigned __int64 c,unsigned __int64 d,unsigned __int64 r,int h)
  5. {
  6.         if(r-b*b-d*d<z)in[h]+=2;
  7.         else
  8.                 if(a*a+c*c-r<z)ou[h]+=2;
  9.                 else
  10.                         if(h++<l)
  11.                                 r*=9,
  12.                                 f(a+a+a,a+a+b,c+c+c,c+c+d,r,h),
  13.                                 f(a+a+b,a+b+b,c+c+c,c+c+d,r,h),
  14.                                 f(a+b+b,b+b+b,c+c+c,c+c+d,r,h),
  15.                                 f(a+a+a,a+a+b,c+c+d,c+d+d,r,h),
  16.                                 f(a+b+b,b+b+b,c+c+d,c+d+d,r,h),
  17.                                 f(a+a+a,a+a+b,c+d+d,d+d+d,r,h),
  18.                                 f(a+a+b,a+b+b,c+d+d,d+d+d,r,h),
  19.                                 f(a+b+b,b+b+b,c+d+d,d+d+d,r,h);
  20. }

  21. int main()
  22. {
  23.         ou[2]=1,in[2]=3;
  24.         f(3,5,7,9,81,2),puts("*"),f(5,7,7,9,81,2),puts("**");
  25.         for(r=9,s=16,i=2;i<=l;i++,s<<=3,r*=3)
  26.         {
  27.                 in[i]+=5;
  28.                 f(1,3,r-2,r,r*r,i);
  29.                 printf("%I64u %I64u\n",in[i],ou[i]);
  30.         }
  31.         return 0;
  32. }
复制代码

虽然坐标的平方和会超出$64bit$整数范围,但是与圆半径的平方进行比较的时候,两者的差异不会超过$2^63$,因此只看平方运算结果的最后$64bit$足矣
于是跑32层迭代也是妥妥的,根本无需进行更高精度的运算,只是需要18天罢了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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)$维 的地毯 混在一块到底该作何理解?


毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-29 19:06:01 | 显示全部楼层
wayne 发表于 2017-9-29 17:13
仔细琢磨起来,这个维度怎么来度量有点名堂,突然间觉得我们在这编程失去了参考意义。
首先比较平凡的结论 ...

这个白格子是白格子。黑格子却不是黑格子,仍然是空洞的。姑且算“灰格子”吧。
这些灰格子不是二维还是\(Log_3(8)\)维几何体。所以前面的编程还是意义的,算是数值上进行勒贝格积分吧。

点评

灰格子最初是全连通的,在二维空间,无穷迭代后,崩塌成二维的不连续。退化至$log_3(8)$维空间。接下来我们就束手无策,不知道能干啥了,^_^  发表于 2017-9-30 13:08
灰点之间是连续的,倒是不同白格子之间是不连通的  发表于 2017-9-30 11:18
每个点都是 灰点,没有连续的概念  发表于 2017-9-29 21:05
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-9-29 22:03:48 | 显示全部楼层
点就是0维的点。
平面是空洞的,但有连续线段,就是每个谢尔宾斯基地毯的纵横 1/3,2/3分割线。

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-10-5 00:09:30 | 显示全部楼层
wayne 发表于 2017-9-28 13:29
......
修改Fans的代码,……,2天12小时27分完成30次迭代,预计在2017年10月23日完成32次迭代。


等不及$32$次迭代了,先处理一下前$30$项数据:
  1. {0,8,0}
  2. {4,28,32}
  3. {104,76,332}
  4. {984,204,2908}
  5. {8288,580,23900}
  6. {67744,1556,192844}
  7. {545904,4180,1547068}
  8. {4377944,11204,12388068}
  9. {35052688,29724,99135316}
  10. {280499768,79276,793162780}
  11. {2244206376,212076,6345516140}
  12. {17954209288,565692,50764701756}
  13. {143635171632,1509332,406119132924}
  14. {1149085383304,4026028,3248957101772}
  15. {9192693786392,10740796,25991667561644}
  16. {73541578891856,28646804,207933369171996}
  17. {588332707461496,76396620,1663467029827132}
  18. {4706661863240488,203728972,13307736442512524}
  19. {37653295449008288,543283204,106461892083564380}
  20. {301226365040331144,1448779164,851695138117736668}
  21. {2409810924185402784,3863345612,6813561108806027412}
  22. {19278487403784147304,10302538780,54508488880751520380}
  23. {154227899257743649824,27473690092,436067911073488311796}
  24. {1233823194135208730288,73263231116,3488543288661173252292}
  25. {9870585553277031899992,195369181668,27908346309484760627908}
  26. {78964684426737227956608,520985280228,223266770476399080439708}
  27. {631717475415287101141792,1389296277316,1786134163812581951993244}
  28. {5053739803326001562892000,3704793953044,14289073310504360438453772}
  29. {40429918426617891895941808,9879454623900,114312586484044763011824820}
  30. {323439347412969480294205216,26345219429236,914500691872384449385489772}
复制代码

首先按照hujunhua在$27$楼求面积上下界的代数平均数的方法,求得以下面积序列:
  1. 0.5
  2. 0.71875
  3. 0.72265625
  4. 0.73486328125
  5. 0.73822021484375
  6. 0.73860931396484375
  7. 0.73869609832763671875
  8. 0.73872029781341552734375
  9. 0.73872639238834381103515625
  10. 0.73872731812298297882080078125
  11. 0.73872764804400503635406494140625
  12. 0.73872775249765254557132720947266
  13. 0.73872777209544437937438488006592
  14. 0.73872777529368249815888702869415
  15. 0.73872777568459468966466374695301
  16. 0.7387277758233707913859689142555
  17. 0.73872777585101889741281411261298
  18. 0.73872777586064863886150533289765
  19. 0.73872777586196634869164512338102
  20. 0.73872777586238128934292834770758
  21. 0.73872777586245612292843720769753
  22. 0.73872777586247222368306432349616
  23. 0.73872777586247703259742587895975
  24. 0.73872777586247779770407592540328
  25. 0.73872777586247796382687745843604
  26. 0.73872777586247800131432562497993
  27. 0.7387277758624780077222150705618
  28. 0.73872777586247800949936495034918
  29. 0.73872777586247800985928169095094
  30. 0.73872777586247800992560980652542
复制代码

然后对上述面积序列的相邻两项作差,得到:
  1. 0.5
  2. 0.21875
  3. 0.00390625
  4. 0.01220703125
  5. 0.00335693359375
  6. 0.00038909912109375
  7. 0.00008678436279296875
  8. 0.00002419948577880859375
  9. 0.00000609457492828369140625
  10. 0.00000092573463916778564453125
  11. 0.00000032992102205753326416015625
  12. 0.00000010445364750921726226806641
  13. 0.00000001959779183380305767059326
  14. 0.00000000319823811878450214862823
  15. 0.00000000039091219150577671825886
  16. 0.00000000013877610172130516730249
  17. 0.00000000002764810602684519835748
  18. 0.00000000000962974144869122028467
  19. 0.00000000000131770983013979048337
  20. 0.00000000000041494065128322432656
  21. 0.00000000000007483358550885998995
  22. 0.00000000000001610075462711579863
  23. 0.00000000000000480891436155546359
  24. 0.00000000000000076510665004644353
  25. 0.00000000000000016612280153303276
  26. 0.00000000000000003748744816654389
  27. 0.00000000000000000640788944558187
  28. 0.00000000000000000177714987978738
  29. 0.00000000000000000035991674060176
  30. 0.00000000000000000006632811557448
复制代码

再对上述序列相邻两项作商,得到:
  1. 1 2.2857142857142857142857142857143
  2. 2 56
  3. 3 0.32
  4. 4 3.6363636363636363636363636363636
  5. 5 8.6274509803921568627450980392157
  6. 6 4.4835164835164835164835164835165
  7. 7 3.5862068965517241379310344827586
  8. 8 3.9706601466992665036674816625917
  9. 9 6.5835010060362173038229376257545
  10. 10 2.8059280169371912491178546224418
  11. 11 3.1585399832822513234884367771037
  12. 12 5.3298682012251717096714318820994
  13. 13 6.1276837764822977392293561428623
  14. 14 8.1814744801512287334592972660495
  15. 15 2.8168552557472735651016350792179
  16. 16 5.0193710045295383725782481751307
  17. 17 2.871116132676942942459197582619
  18. 18 7.3079377784330865393824639461471
  19. 19 3.1756585575906052333539618816487
  20. 20 5.5448452517900141689032828041249
  21. 21 4.6478309397269083352206042808509
  22. 22 3.3481059167599611362212550782029
  23. 23 6.285286320885529757449321963342
  24. 24 4.6056690772476863673157818475219
  25. 25 4.4314246409893269354527425736722
  26. 26 5.8502020805603697358084394450449
  27. 27 3.6057113237677602868585857216454
  28. 28 4.9376694088085151486851873689654
  29. 29 5.4263073431900628788285733338051
复制代码

然后根据上述序列作散点图,得到:

data.png

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

$0.738727775862478009946(3)$

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

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

评分

参与人数 1威望 +3 金币 +3 贡献 +3 经验 +3 鲜花 +3 收起 理由
wayne + 3 + 3 + 3 + 3 + 3

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-10-12 19:01:44 | 显示全部楼层
抛砖引玉说一下我的极限推导思路。就是仿照古典积分从矩形逼近曲线的过程。
把地毯切割,作出累积曲线。切割地毯1维宽度每缩小1/3,累积的 谢尔宾斯基地毯缩小为1/8。

curveS.png


未切割的累积曲线就是

未命名a-1.png

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

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

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



毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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 类似。
  1. #include<stdio.h>
  2. #include<stdint.h>
  3. uint64_t in[64]={0},out[64]={0};
  4. int maxl=24;
  5. //fh:=(x+1)^2-(y+1)^2-r^2
  6. //fl:=(x-1)^2-(y-1)^2-r^2
  7. void ffull(int l,int64_t x,int64_t y,int64_t fh,int64_t fl){
  8.         if(++l==maxl){
  9.                 in[l]+=32;
  10.                 out[l]+=32;
  11.                 return;
  12.         };
  13.         x*=3;
  14.         y*=3;
  15.         fh*=9;
  16.         fl*=9;
  17.         int64_t nfh,nfl;
  18. #define proc(nx,ny,dfh,dfl)          \
  19. do{                                  \
  20.         nfh=fh-4*(dfh),nfl=fl+4*(dfl);   \
  21.         if(nfl>=0)out[l]+=8;             \
  22.         else if(nfh<=0)in[l]+=8;         \
  23.         else ffull(l,nx,ny,nfh,nfl);     \
  24. }while(0)

  25.         proc(x+2,y+2,0,2*x+2*y-4);
  26.         proc(x+2,y,y+2,2*x+y-4);
  27.         proc(x+2,y-2,2*y+2,2*x-2);
  28.         proc(x,y+2,x+2,x+2*y-4);
  29.         proc(x,y-2,x+2*y+4,x-2);
  30.         proc(x-2,y+2,2*x+2,2*y-2);
  31.         proc(x-2,y,2*x+y+4,y-2);
  32.         proc(x-2,y-2,2*x+2*y+4,0);
  33. #undef proc
  34. }

  35. void fhalf(int l,int64_t x,int64_t fh,int64_t fl){
  36.         if(++l==maxl){
  37.                 in[l]+=32;
  38.                 return;
  39.         };
  40.         x*=3;
  41.         fh*=9;
  42.         fl=9*fl+4*(2*x-4);
  43.         in[l]+=20;
  44.         ffull(l,x+2,2,fh,fl);
  45.         fhalf(l,x+2,fh-8,fl);
  46. }

  47. int main(){
  48.         in[2]=32,out[2]=4;
  49.         ffull(2,8,6,49,-7);
  50.         printf("*");
  51.         ffull(2,8,4,25,-23);
  52.         printf("*");
  53.         ffull(2,8,2,9,-31);
  54.         printf("*");
  55.         fhalf(2,8,1,-31);
  56.         double inres=0,outres=0;
  57.         for(int i=2;i<=maxl;++i){
  58.                 inres=8*inres+in[i];
  59.                 outres=8*outres+out[i];
  60.                 printf("%18.16lf\n",inres/(outres+inres));
  61.         }
  62.         return 0;
  63. }
复制代码


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

点评

额……我发现我把注释和解释里 fh 和 fl 的定义写错了…… 应该是 x^2+y^2-r^2 的形式,我经常把第一个加号写成减号……不过代码是对的。而且正在跑的代码和贴出来的有区别,最终结果将是35楼的数据表格式。  发表于 2017-10-19 00:05
我也发现我原来预计错了,不过我做出的新预计是2天半左右。我从 10月17日22时多 正式开始计算,仍然在等待结果中……  发表于 2017-10-19 00:00
如果24层只需55秒,那么32层只需1.6天  发表于 2017-10-18 22:56
我搞错了层数……事实上我应该说23层单线程不到一分钟的。 我的电脑有四线程,所以我进行了暴力三线程并行。现在24层大约耗时55秒。 我即将开始进行32层的计算,预计需要4到5天。  发表于 2017-10-17 18:05
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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
结果:
  1. {0, 8, 0}
  2. {4, 28, 32}
  3. {104, 76, 332}
  4. {984, 204, 2908}
  5. {8288, 580, 23900}
  6. {67744, 1556, 192844}
  7. {545904, 4180, 1547068}
  8. {4377944, 11204, 12388068}
  9. {35052688, 29724, 99135316}
  10. {280499768, 79276, 793162780}
  11. {2244206376, 212076, 6345516140}
  12. {17954209288, 565692, 50764701756}
  13. {143635171632, 1509332, 406119132924}
  14. {1149085383304, 4026028, 3248957101772}
  15. {9192693786392, 10740796, 25991667561644}
  16. {73541578891856, 28646804, 207933369171996}
  17. {588332707461496, 76396620, 1663467029827132}
  18. {4706661863240488, 203728972, 13307736442512524}
  19. {37653295449008288, 543283204, 106461892083564380}
  20. {301226365040331144, 1448779164, 851695138117736668}
  21. {2409810924185402784, 3863345612, 6813561108806027412}
  22. {19278487403784147304, 10302538780, 54508488880751520380}
  23. {154227899257743649824, 27473690092, 436067911073488311796}
  24. {1233823194135208730288, 73263231116, 3488543288661173252292}
  25. {9870585553277031899992, 195369181668, 27908346309484760627908}
  26. {78964684426737227956608, 520985280228, 223266770476399080439708}
  27. {631717475415287101141792, 1389296277316, 1786134163812581951993244}
  28. {5053739803326001562892000, 3704793953044, 14289073310504360438453772}
  29. {40429918426617891895941808, 9879454623900, 114312586484044763011824820}
  30. {323439347412969480294205216, 26345219429236, 914500691872384449385489772}
  31. {2587514779303826096159159104, 70253935029820, 7316005534979145849098804868}
  32. {20700118234430796112810254760, 187343834543092, 58528044279833354136899152484}
复制代码

点评

耗时果然是1.6天耶~你的程序太棒了!  发表于 2017-10-19 22:14

评分

参与人数 2威望 +12 金币 +12 贡献 +12 经验 +12 鲜花 +12 收起 理由
wayne + 6 + 6 + 6 + 6 + 6 赞一个!
KeyTo9_Fans + 6 + 6 + 6 + 6 + 6 赞一个!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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$

散点图:

data.png

从图中可以看到,

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

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

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

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

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

#####

根据最新的$2$个数据重新预测极限,结果为($21$个有效数字):

$0.7387277758624780099409(3)$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-11-3 08:38:55 | 显示全部楼层
哇!这个贴子变成精品贴了?好开心呀~

#####

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

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

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

#####

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

点评

"Number of level $n$ squares on a Sierpinski carpet that are enclosed in a circle with the same center and diameter",翻译过来是"谢尔宾斯基地毯在同径同心圆内的$n$级小正方形个数",已经说得很清   发表于 2017-11-3 12:41
评委们表示看懂了这个数列的意义了吗?  发表于 2017-11-3 10:35
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-10-18 20:31 , Processed in 0.071425 second(s), 18 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表