计算半径为 R 且尺寸为 D 的球体内的整数点答案

【问题标题】：Counting integer points inside a sphere of radius R and dimension D计算半径为 R 且尺寸为 D 的球体内的整数点
【发布时间】：2016-09-21 01:03:15
【问题描述】：

我正在尝试编写一个有效的算法来计算半径为 R 和维度 D 的球体内的点数。球体始终位于原点。假设我们有一个维度为 2（圆）、半径为 5 的球体。

我的策略是在第一象限内生成所有可能的点，因此对于上面的示例，我们知道 (1,2) 在圆中，因此该点的所有 + / - 组合也必须是简单的维度平方。因此，对于在 n 维球体的单个象限中找到的每个点，我们将 2 ^ 维添加到总数中。

我不确定是否有更有效的解决方案来解决这个问题，但这是我迄今为止在实施方面所拥有的。

int count_lattice_points(const double radius, const int dimension) {
int R = static_cast<int>(radius);

int count  = 0;

std::vector<int> points;
std::vector<int> point;

for(int i = 0; i <= R; i++)
    points.push_back(i);

do {
    for(int i = 0; i < dimension - 1; i++)
        point.push_back(points.at(i));

    if(isPointWithinSphere(point, radius)) count += std::pow(2,dimension);
    point.clear();

}while(std::next_permutation(points.begin(), points.end()));

return count + 3;
}

在这种情况下我可以解决或改进什么？

【问题讨论】：

你是假设它的原点总是在 (0, 0) 还是 (0, 0, 0)？如果没有这种限制，您当前的推理似乎无法正常工作。
@JamesAdkison 是的，球体总是在原点
代码 sn -p 的输出没有给出二维的预期结果，见The On-line Encyclopedia of Integer Sequences。我看到一些问题。一个点的坐标应该是dimension，但是第二个for 循环中的i < dimension - 1 会导致测试点中的坐标太少。当某些坐标为零时，添加 2^dimension 是不正确的，因为零的 +/- 只是零。它应该是 2^（正坐标的数量）。另外，如果只是置换值，我不明白如何计算像 (1,1) 这样的点。

标签： c++ algorithm geometry computational-geometry

【解决方案1】：

对于 2D 情况，这是 Gauss's circle problem. 一个可能的公式：

N(r) = 1 + 4 * r + 4 * Sum[i=1..r]{Floor(Sqrt(r^2-i^2))}

（中心点+四个象限，4*r为轴上的点，其他为象限内）。

请注意，对于 2D 情况，没有已知的简单封闭数学表达式。

一般而言，您对象限、八分圆等的想法是正确的，但检查所有点太昂贵了。

人们可能会找到从 0 到 r^2 从 1..D 组成所有正方形的方法数整数平方（(4)式的扩展）.

请注意，组合数学有助于加快计算速度。例如，找到方法的数量就足够了由 D 个自然平方生成 X^2，并乘以 2^D（不同的符号组合）；找出从 D-1 个自然平方生成 X^2 的方法数，然后乘以 D*2^(D-1)（不同的符号组合 + D 个零加数的位）等

D=2、R=3 的示例

addends: 0,1,4,9
possible sum     compositions    number of variants        
0               0+0             1
1               0+1,1+0         2*2=4
2               1+1             4      
4               0+4,4+0         2*2=4
5               1+4,4+1         2*4=8  
8               4+4             4
9               0+9,9+0         2*2=4
-------------------------------------
                                29

【讨论】：

你能举一个相同半径和尺寸的具体例子吗？
r=5 的示例太长，添加了 r=3
我只有最后一个请求，你能知道伪代码/实现在 R 和 D 方面是什么样的吗？
您需要将整数划分为已知部分（正方形）的算法，用于重复排列的公式。我会使用动态编程构建一个带有可能分区的数组 [0..R^2] 以避免多次计算..
可以在monsiterdex.wordpress.com/2013/04/05/… 找到类似的方法，包括源代码。该方法包括找到半径的分区，然后为球体中的每个分区计算可以通过置换坐标和翻转非零坐标的符号在球体中表示的方式的数量。

【解决方案2】：

类似于that described by MBo 的方法，包括源代码，可以在以下位置找到 https://monsiterdex.wordpress.com/2013/04/05/integer-lattice-in-n-dimensional-sphere-count-of-points-with-integer-coordinates-using-parallel-programming-part-i/.

该方法包括找到半径的分区，然后为球体中的每个分区计算可以通过置换坐标和翻转非零坐标的符号在球体中表示的方式数。

【讨论】：

@BhargavRao，我删除这个答案会更好吗？我知道对 MBo 的回答发表评论会更好，但我不能这样做。
不，别说了。它看起来好多了atm。我想一定有人看到你最初的回答投了反对票。
这并没有提供问题的答案。要批评或要求作者澄清，请在他们的帖子下方留下评论。 - From Review
@SebastianLenartowicz 在撰写本文时我无法发表评论。其次，链接详细说明了一个解决方案，比Mbo 的答案更彻底

【解决方案3】：

我在这里展示了我的 2D 算法（带有一些源代码和一个丑陋但方便的插图）： https://stackoverflow.com/a/42373448/5298879

它比 MBo's 在其中一个季度计算圆的原点和边缘之间的点快大约 3.4 倍。

您只需想象一个内接正方形，然后只计算该正方形外的八分之一。

public static int gaussCircleProblem(int radius) {
    int allPoints=0; //holds the sum of points
    double y=0; //will hold the precise y coordinate of a point on the circle edge for a given x coordinate.
    long inscribedSquare=(long) Math.sqrt(radius*radius/2); //the length of the side of an inscribed square in the upper right quarter of the circle
    int x=(int)inscribedSquare; //will hold x coordinate - starts on the edge of the inscribed square
    while(x<=radius){
        allPoints+=(long) y; //returns floor of y, which is initially 0
        x++; //because we need to start behind the inscribed square and move outwards from there
        y=Math.sqrt(radius*radius-x*x); // Pythagorean equation - returns how many points there are vertically between the X axis and the edge of the circle for given x
    }
    allPoints*=8; //because we were counting points in the right half of the upper right corner of that circle, so we had just one-eightth
    allPoints+=(4*inscribedSquare*inscribedSquare); //how many points there are in the inscribed square
    allPoints+=(4*radius+1); //the loop and the inscribed square calculations did not touch the points on the axis and in the center
    return allPoints;
}

【讨论】：