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

倪举鹏

高斯不是有个公式计算原点圆内整数坐标点个数的么，可以利用啊

KeyTo9_Fans

与本贴问题相关的两个数列已经被【在线整数数列百科大全】收录了，

数列编号为A293288和A293289。

前者是圆内小正方形个数，后者是圆上小正方形个数。

这两个数列都更新到了第$32$项，都是直接使用了Ickiverar在$38$楼提供的数据，

感谢Ickiverar、Chyanog和Wayne提供的宝贵素材：

LINKS:
Yi Yang: Table of n, a(n) for n = 0..32
Chyanog: Illustration of initial terms
Wayne, Ickiverar: The C++ code generating the sequence

其中【Table of n, a(n) for n = 0..32】的贡献者应该是Ickiverar，但是我提交的时候没写上去，所以系统就自动填了提交者的名字。

等我拿到更多素材去更新这个数列的时候，看看能否把这个地方改掉。

wayne

KeyTo9_Fans 发表于 2017-11-16 22:58
与本贴问题相关的两个数列已经被【在线整数数列百科大全】收录了，

数列编号为A293288和A293289。

搜了下, Yi Yang总共贡献了14个数列.
---------------
插曲:  最开始搜索的时候没有加引号,一不小心,搜到了 https://oeis.org/A102241.          Yi Jing Algebra, 好神奇:http://www.yijing.co.uk ,

chyanog

这个问题可以推广到三维吗？即门格海绵（Menger sponge）与球的交集。我画了几幅图


补充内容 (2019-12-2 21:27):
http://www.lactamme.polytechniqu ... .5.1.D/display.html 已经有人画过这样的图了

hujunhua

在3维的情况下，21#的担心成为现实，会有两个对顶小正方体分居球面内外，公共顶点刚好位于球面上。
也会有小正方体的中心正好位于球面上。\[x^2+y^2+z^2=9^n,xyz≠0\]总是有解的，记解数为`t_n`，不知有无现成的数列。

chyanog

迭代次数，球外，球上，球内的方块的个数
{1, 0, 20, 0}
{2, 140, 216, 44}
{3, 4456, 1224, 2320}
{4, 97232, 7968, 54800}
{5, 1999432, 54456, 1146112}
{6, 40353584, 367944, 23278472}

chyanog

chyanog 发表于 2019-8-15 15:50
迭代次数，球外，球上，球内的方块的个数
{1, 0, 20, 0}
{2, 140, 216, 44}

更新到第8项
(1), 0, 20, 0
(2), 140, 216, 44
(3), 4456, 1224, 2320
(4), 97232, 7968, 54800
(5), 1999432, 54456, 1146112
(6), 40353584, 367944, 23278472
(7), 809527216, 2444520, 468028264
(8), 16206910376, 16284576, 9376805048

代码没有根据对称性作优化，时间复杂度O(20^n)，还有很大优化空间，计算到第8项需要十五六分钟。

1. #include <stdio.h>
   
2. #include <math.h>
   
3. #include <time.h>
   
4. 
   
5. #define MAX(x, y) (((x) > (y)) ? (x) : (y))
   
6. #define MIN(x, y) (((x) < (y)) ? (x) : (y))
   
7. 
   
8. unsigned __int64 in_count, out_count;
   
9. 
   
10. void f(double x, double y, double z, int m, int n)
    
11. {
    
12.         double d = 2 * pow(3, -n);
    
13.         if (n == 0)
    
14.         {
    
15.                 double r = pow(3, -m);
    
16.                 double d1, d2, d3, d4, d5, d6, d7, d8;
    
17.                 d1 = pow(x + r, 2) + pow(y + r, 2) + pow(z + r, 2);
    
18.                 d2 = pow(x + r, 2) + pow(y + r, 2) + pow(z - r, 2);
    
19.                 d3 = pow(x + r, 2) + pow(y - r, 2) + pow(z + r, 2);
    
20.                 d4 = pow(x + r, 2) + pow(y - r, 2) + pow(z - r, 2);
    
21.                 d5 = pow(x - r, 2) + pow(y + r, 2) + pow(z + r, 2);
    
22.                 d6 = pow(x - r, 2) + pow(y + r, 2) + pow(z - r, 2);
    
23.                 d7 = pow(x - r, 2) + pow(y - r, 2) + pow(z + r, 2);
    
24.                 d8 = pow(x - r, 2) + pow(y - r, 2) + pow(z - r, 2);
    
25.                 if (MAX(MAX(MAX(MAX(MAX(MAX(MAX(d1, d2), d3), d4), d5), d6), d7), d8) < 1)
    
26.                 {
    
27.                         in_count++;
    
28.                 }
    
29.                 if (MIN(MIN(MIN(MIN(MIN(MIN(MIN(d1, d2), d3), d4), d5), d6), d7), d8) > 1)
    
30.                 {
    
31.                         out_count++;
    
32.                 }
    
33.         }
    
34.         else
    
35.         {
    
36.                 f(x - d, y - d, z - d, m, n - 1);
    
37.                 f(x - d, y - d, z, m, n - 1);
    
38.                 f(x - d, y - d, z + d, m, n - 1);
    
39.                 f(x - d, y, z - d, m, n - 1);
    
40.                 f(x - d, y, z + d, m, n - 1);
    
41.                 f(x - d, y + d, z - d, m, n - 1);
    
42.                 f(x - d, y + d, z, m, n - 1);
    
43.                 f(x - d, y + d, z + d, m, n - 1);
    
44.                 f(x, y - d, z - d, m, n - 1);
    
45.                 f(x, y - d, z + d, m, n - 1);
    
46.                 f(x, y + d, z - d, m, n - 1);
    
47.                 f(x, y + d, z + d, m, n - 1);
    
48.                 f(x + d, y - d, z - d, m, n - 1);
    
49.                 f(x + d, y - d, z, m, n - 1);
    
50.                 f(x + d, y - d, z + d, m, n - 1);
    
51.                 f(x + d, y, z - d, m, n - 1);
    
52.                 f(x + d, y, z + d, m, n - 1);
    
53.                 f(x + d, y + d, z - d, m, n - 1);
    
54.                 f(x + d, y + d, z, m, n - 1);
    
55.                 f(x + d, y + d, z + d, m, n - 1);
    
56.         }
    
57. }
    
58. 
    
59. int main()
    
60. {
    
61.         clock_t t0 = clock();
    
62.         for (int i = 1; i <= 6; i++)
    
63.         {
    
64.                 in_count = 0;
    
65.                 out_count = 0;
    
66.                 f(0, 0, 0, i, i);
    
67.                 printf("(%d), %llu, %llu, %llu\n", i, out_count, (__int64)pow(20, i) - out_count - in_count, in_count);
    
68.         }
    
69.         printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);
    
70.         return 0;
    
71. }
    

复制代码

Ickiverar

嗯？又被挖出来了。
需要我写用于三维的代码吗

mathe

这个序列会需要更多数据

mathe

https://oeis.org/A329302