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

Ickiverar

用于三维的c++代码，在我的电脑上单线程算到11层需要3秒。
要算更多层，只需要更改const int maxl=层数;

1. #include<stdint.h>
   
2. #include<stdio.h>
   
3. 
   
4. const int maxl=11;
   
5. int64_t in[64]={0},on[64]={0},out[64]={0};
   
6. int64_t g;
   
7. 
   
8. #define xm(proc,x,y,z)                  \
   
9.     proc(x-2,y-2,z-2,2*x+2*y+2*z+6,0);  \
   
10.     proc(x-2,y-2,z,2*x+2*y+z+6,z-2);    \
    
11.     proc(x-2,y-2,z+2,2*x+2*y+4,2*z-2);  \
    
12.     proc(x-2,y,z-2,2*x+y+2*z+6,y-2);    \
    
13.     proc(x-2,y,z+2,2*x+y+4,y+2*z-4);    \
    
14.     proc(x-2,y+2,z-2,2*x+2*z+4,2*y-2);  \
    
15.     proc(x-2,y+2,z,2*x+z+4,2*y+z-4);    \
    
16.     proc(x-2,y+2,z+2,2*x+2,2*y+2*z-4);
    
17. 
    
18. #define x0(proc,x,y,z)                  \
    
19.     proc(x,y-2,z-2,x+2*y+2*z+6,x-2);    \
    
20.     proc(x,y-2,z+2,x+2*y+4,x+2*z-4);    \
    
21.     proc(x,y+2,z-2,x+2*z+4,x+2*y-4);    \
    
22.     proc(x,y+2,z+2,x+2,x+2*y+2*z-6);   
    
23. 
    
24. #define xp(proc,x,y,z)                  \
    
25.     proc(x+2,y-2,z-2,2*y+2*z+4,2*x-2);  \
    
26.     proc(x+2,y-2,z,2*y+z+4,2*x+z-4);    \
    
27.     proc(x+2,y-2,z+2,2*y+2,2*x+2*z-4);  \
    
28.     proc(x+2,y,z-2,y+2*z+4,2*x+y-4);    \
    
29.     proc(x+2,y,z+2,y+2,2*x+y+2*z-6);    \
    
30.     proc(x+2,y+2,z-2,2*z+2,2*x+2*y-4);  \
    
31.     proc(x+2,y+2,z,z+2,2*x+2*y+z-6);    \
    
32.     proc(x+2,y+2,z+2,0,2*x+2*y+2*z-6);
    
33. 
    
34. void ffull(int l,int64_t x,int64_t y,int64_t z,int64_t fh,int64_t fl){
    
35.     ++l;
    
36.     x*=3;
    
37.     y*=3;
    
38.     z*=3;
    
39.     fh*=9;
    
40.     fl*=9;
    
41.     int64_t nfh,nfl;
    
42. #define proc(nx,ny,nz,dfh,dfl)          \
    
43. do{                                     \
    
44.     nfh=fh-4*(dfh),nfl=fl+4*(dfl);      \
    
45.     if(nfl>=0)out[l]+=g;                \
    
46.     else if(nfh<=0)in[l]+=g;            \
    
47.     else {                              \
    
48.         on[l]+=g;                       \
    
49.         if(l<maxl)                      \
    
50.             ffull(l,nx,ny,nz,nfh,nfl);  \
    
51.     }                                   \
    
52. }while(0)
    
53. 
    
54.     xm(proc,x,y,z);
    
55.     x0(proc,x,y,z);
    
56.     xp(proc,x,y,z);
    
57. 
    
58. }
    
59. 
    
60. void fhalf(int l,int64_t y,int64_t z,int64_t fh,int64_t fl){
    
61.     const int64_t x=0;
    
62.     ++l;
    
63.     y*=3;
    
64.     z*=3;
    
65.     fh*=9;
    
66.     fl*=9;
    
67.     int64_t nfh,nfl;
    
68. #define proch(nx,ny,nz,dfh,dfl)         \
    
69. do{                                     \
    
70.     nfh=fh-4*(dfh),nfl=fl+4*(dfl);      \
    
71.     if(nfl>=0)out[l]+=g/2;              \
    
72.     else if(nfh<=0)in[l]+=g/2;          \
    
73.     else {                              \
    
74.         on[l]+=g/2;                     \
    
75.         if(l<maxl)                      \
    
76.             fhalf(l,ny,nz,nfh,nfl);     \
    
77.     }                                   \
    
78. }while(0)
    
79. 
    
80.     x0(proch,x,y,z);
    
81.     xp(proc,x,y,z);
    
82. }
    
83. 
    
84. int main(){
    
85.     on[0]=1;
    
86.     on[1]=20;
    
87.     g=8;
    
88.     ffull(1,2,2,2,18,-6);
    
89.     g=24;
    
90.     fhalf(1,  2,2,10,-6);
    
91.     for(int i=1;i<=maxl;++i){
    
92.         /*
    
93.             Note: To avoid overflow,
    
94.             in[i] & out[i] are not real counts of level i cubes
    
95.             unless after executing following lines :
    
96. 
    
97.         in[i]+=20*in[i-1];
    
98.         out[i]+=20*out[i-1];
    
99. 
    
100.             This can be done in Mathematica
     
101.         */
     
102.         printf("%llu %llu %llu\n",in[i],on[i],out[i]);
     
103.     }
     
104. }
     

复制代码

mathe

我试着运行了一下，和chyang的结果不同
0 20 0
44 216 140
1440 1224 1656
8400 7968 8112
50112 54456 54792
356232 367944 364944
2458824 2444520 2455536
16239768 16284576 16366056
108395712 108525168 108770640
723385080 723469800 723648480
4823029008 4822582872 4823784120

mathe

0 20 0
44 216 140
2320 1224 4456
54800 7968 97232
1146112 54456 1999432
23278472 367944 40353584
468028264 2444520 809527216
9376805048 16284576 16206910376
187644496672 108525168 324246978160
3753613318520 723469800 6485663211680
75077089399408 4822582872 129718088017720
1501573939443704 32150080560 2594393910475736
30031693123408792 214334797104 51888092541794104

Ickiverar

Ickiverar 发表于 2019-12-9 18:29
用于三维的c++代码，在我的电脑上单线程算到11层需要3秒。
要算更多层，只需要更改const int maxl=层数;
...

8线程并行写完了，20秒13层，预计3天能算到18层。开算。

mathe

能够分别计算结果关于模2^64和2^64-1的结果吗？如果这样，最后就不怕溢出，最后通过中国剩余定理可以得出更大范围的结果了。
毕竟还是指数递增的，差值好像公差接近$\frac{20}3$,18项应该没有问题，但是计算更多项就有问题了。

Ickiverar

18层的数据
内部 球面上 外部

1. 0 20 0
   
2. 44 216 140
   
3. 2320 1224 4456
   
4. 54800 7968 97232
   
5. 1146112 54456 1999432
   
6. 23278472 367944 40353584
   
7. 468028264 2444520 809527216
   
8. 9376805048 16284576 16206910376
   
9. 187644496672 108525168 324246978160
   
10. 3753613318520 723469800 6485663211680
    
11. 75077089399408 4822582872 129718088017720
    
12. 1501573939443704 32150080560 2594393910475736
    
13. 30031693123408792 214334797104 51888092541794104
    
14. 600635291349699440 1428893539656 1037763279756760904
    
15. 12012715352891236816 9525966857976 20755275121141905208
    
16. 240254370564131706488 63506419362240 415105565929448931272
    
17. 4805087834658423387496 423376135436184 8302111741965441176320
    
18. 96101759515675744360544 2822507323770168 166042237661816931869288

复制代码

mathe

上面的数据已经更新在 https://oeis.org/A329302