高斯不是有个公式计算原点圆内整数坐标点个数的么，可以利用啊
与本贴问题相关的两个数列已经被【在线整数数列百科大全】收录了，
数列编号为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，但是我提交的时候没写上去，所以系统就自动填了提交者的名字。
等我拿到更多素材去更新这个数列的时候，看看能否把这个地方改掉。
KeyTo9_Fans 发表于 2017-11-16 22:58
与本贴问题相关的两个数列已经被【在线整数数列百科大全】收录了，
数列编号为A293288和A293289。
搜了下, Yi Yang总共贡献了14个数列.
---------------
插曲:最开始搜索的时候没有加引号,一不小心,搜到了 https://oeis.org/A102241. Yi Jing Algebra, 好神奇:http://www.yijing.co.uk ,:lol
这个问题可以推广到三维吗？即门格海绵（Menger sponge）与球的交集。我画了几幅图
补充内容 (2019-12-2 21:27):
http://www.lactamme.polytechnique.fr/Mosaic/images/MENG.C4.5.1.D/display.html 已经有人画过这样的图了
在3维的情况下，21#的担心成为现实，会有两个对顶小正方体分居球面内外，公共顶点刚好位于球面上。
也会有小正方体的中心正好位于球面上。\总是有解的，记解数为`t_n`，不知有无现成的数列。
本帖最后由 chyanog 于 2019-8-15 17:59 编辑
迭代次数，球外，球上，球内的方块的个数
{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 于 2019-12-2 19:35 编辑
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项需要十五六分钟。
#include <stdio.h>
#include <math.h>
#include <time.h>
#define MAX(x, y) (((x) > (y)) ? (x) : (y))
#define MIN(x, y) (((x) < (y)) ? (x) : (y))
unsigned __int64 in_count, out_count;
void f(double x, double y, double z, int m, int n)
{
double d = 2 * pow(3, -n);
if (n == 0)
{
double r = pow(3, -m);
double d1, d2, d3, d4, d5, d6, d7, d8;
d1 = pow(x + r, 2) + pow(y + r, 2) + pow(z + r, 2);
d2 = pow(x + r, 2) + pow(y + r, 2) + pow(z - r, 2);
d3 = pow(x + r, 2) + pow(y - r, 2) + pow(z + r, 2);
d4 = pow(x + r, 2) + pow(y - r, 2) + pow(z - r, 2);
d5 = pow(x - r, 2) + pow(y + r, 2) + pow(z + r, 2);
d6 = pow(x - r, 2) + pow(y + r, 2) + pow(z - r, 2);
d7 = pow(x - r, 2) + pow(y - r, 2) + pow(z + r, 2);
d8 = pow(x - r, 2) + pow(y - r, 2) + pow(z - r, 2);
if (MAX(MAX(MAX(MAX(MAX(MAX(MAX(d1, d2), d3), d4), d5), d6), d7), d8) < 1)
{
in_count++;
}
if (MIN(MIN(MIN(MIN(MIN(MIN(MIN(d1, d2), d3), d4), d5), d6), d7), d8) > 1)
{
out_count++;
}
}
else
{
f(x - d, y - d, z - d, m, n - 1);
f(x - d, y - d, z, m, n - 1);
f(x - d, y - d, z + d, m, n - 1);
f(x - d, y, z - d, m, n - 1);
f(x - d, y, z + d, m, n - 1);
f(x - d, y + d, z - d, m, n - 1);
f(x - d, y + d, z, m, n - 1);
f(x - d, y + d, z + d, m, n - 1);
f(x, y - d, z - d, m, n - 1);
f(x, y - d, z + d, m, n - 1);
f(x, y + d, z - d, m, n - 1);
f(x, y + d, z + d, m, n - 1);
f(x + d, y - d, z - d, m, n - 1);
f(x + d, y - d, z, m, n - 1);
f(x + d, y - d, z + d, m, n - 1);
f(x + d, y, z - d, m, n - 1);
f(x + d, y, z + d, m, n - 1);
f(x + d, y + d, z - d, m, n - 1);
f(x + d, y + d, z, m, n - 1);
f(x + d, y + d, z + d, m, n - 1);
}
}
int main()
{
clock_t t0 = clock();
for (int i = 1; i <= 6; i++)
{
in_count = 0;
out_count = 0;
f(0, 0, 0, i, i);
printf("(%d), %llu, %llu, %llu\n", i, out_count, (__int64)pow(20, i) - out_count - in_count, in_count);
}
printf("Elapsed time %0.4fs\n", (clock() - t0) / 1000.0);
return 0;
}
嗯？又被挖出来了。
需要我写用于三维的代码吗:lol
这个序列会需要更多数据
https://oeis.org/A329302