这是一种可能的方法，灵感来自您问题中的评论。（对于 tl;dr 版本，请跳至point_to_line此答案底部的定义：仅给出第一象限的映射。扩展到整个平面作为一项不太困难的练习。）

你的问题说：

特别是，没有计算光线和正 x 轴之间的每个格点（除非你知道一种聪明的方法来做到这一点而不将它们全部列举出来）。

这里是一个算法来做到这一点不计数枚举点; 它的效率类似于寻找最大公约数的欧几里得算法。我不确定它在多大程度上算作“易于计算”或“聪明”。

假设我们给定一个点(p, q)整数坐标，既p和q正（使点位于在第一象限）。我们不妨也假设q < p，因此该点(p, q)位于 x 轴y = 0和对角线之间y = x：如果我们可以解决位于对角线下方的第一象限的一半的问题，我们可以利用对称性来解决它一般。

写下和M大小的界限，以便在您的示例中我们想要。pqM = 2^31

那么严格在三角形内的格点数为：

• x 轴 y = 0
• y = (q/p)x从原点开始并通过的射线(p, q)，以及
• 垂直线 x = M

是x在(0, M)of 中的整数范围内的总和?qx/p? - 1。

为方便起见，我将删除-1和包含0在总和的范围内；这些变化都是微不足道的。现在，我们所需要的核心功能是评价的总和的能力?qx/p?作为x在在区间的整数范围[0, M)。在此过程中，我们可能还希望能够计算出一个密切相关的总和：在?qx/p?相同范围内的总和x（结果证明将这两个值一起计算是有意义的）。

出于测试目的，这里是我们感兴趣的函数的缓慢、朴素但显然正确的版本，这里是用 Python 编写的：

def floor_sum_slow(p, q, M):
    """
    Sum of floor(q * x / p) for 0 <= x < M.

    Assumes p positive, q and M nonnegative.
    """
    return sum(q * x // p for x in range(M))


def ceil_sum_slow(p, q, M):
    """
    Sum of ceil(q * x / p) for 0 <= x < M.

    Assumes p positive, q and M nonnegative.
    """
    return sum((q * x + p - 1) // p for x in range(M))

和一个例子使用：

>>> floor_sum_slow(51, 43, 2**28)  # takes several seconds to complete
30377220771239253
>>> ceil_sum_slow(140552068, 161600507, 2**28)
41424305916577422

可以更快地评估这些总和。第一个关键观察是如果q >= p，那么我们可以应用欧几里德“除法算法”并写出q = ap + r一些整数a和r。总和然后简化：该ap部分贡献了一个因子a * M * (M - 1) // 2，我们从计算减少floor_sum(p, q, M)到计算floor_sum(p, r, M)。类似地， 的计算ceil_sum(p, q, M)减少到 的计算ceil_sum(p, q % p, M)。

第二个关键观察是我们可以用 表示floor_sum(p, q, M)，的上限ceil_sum(q, p, N)在哪里。为此，我们考虑矩形，并使用线将该矩形分成两个三角形。矩形内位于直线上或直线下方的格点数为，而位于直线上方的矩形内格点数为。由于晶格点的矩形中的总数是，我们可以推断出的值从的反之亦然，和副。N(q/p)M[0, M) x (0, (q/p)M)y = (q/p)xfloor_sum(p, q, M)ceil_sum(q, p, N)(N - 1)Mfloor_sum(p, q, M)ceil_sum(q, p, N)

结合这两个想法，并仔细研究细节，我们最终得到了一对相互递归的函数，如下所示：

def floor_sum(p, q, M):
    """
    Sum of floor(q * x / p) for 0 <= x < M.

    Assumes p positive, q and M nonnegative.
    """
    a = q // p
    r = q % p
    if r == 0:
        return a * M * (M - 1) // 2
    else:
        N = (M * r + p - 1) // p
        return a * M * (M - 1) // 2 + (N - 1) * M - ceil_sum(r, p, N)


def ceil_sum(p, q, M):
    """
    Sum of ceil(q * x / p) for 0 <= x < M.

    Assumes p positive, q and M nonnegative.
    """
    a = q // p
    r = q % p
    if r == 0:
        return a * M * (M - 1) // 2
    else:
        N = (M * r + p - 1) // p
        return a * M * (M - 1) // 2 + N * (M - 1) - floor_sum(r, p, N)

执行与之前相同的计算，我们得到完全相同的结果，但这次结果是即时的：

>>> floor_sum(51, 43, 2**28)
30377220771239253
>>> ceil_sum(140552068, 161600507, 2**28)
41424305916577422

一些实验应该让您相信floor_sum和floor_sum_slow函数在所有情况下都会给出相同的结果，对于ceil_sum和 也是如此ceil_sum_slow。

这是一个使用floor_sum并ceil_sum为第一象限提供适当映射的函数。我没能抗拒让它成为一个完整的双射的诱惑，按照它们出现在每条射线上的顺序枚举点，但是你可以通过简单地将两个分支中的+ gcd(p, q)术语替换为 来解决这个问题+ 1。

from math import gcd

def point_to_line(p, q, M):
    """
    Bijection from [0, M) x [0, M) to [0, M^2), preserving
    the 'angle' ordering.
    """
    if p == q == 0:
        return 0
    elif q <= p:
        return ceil_sum(p, q, M) + gcd(p, q)
    else:
        return M * (M - 1) - floor_sum(q, p, M) + gcd(p, q)

扩展到整个平面应该很简单，尽管由于二进制补码表示中的负范围和正范围之间的不对称性而有点混乱。

这是 case 的可视化演示M = 7，使用以下代码打印：

M = 7
for q in reversed(range(M)):
    for p in range(M):
        print(" {:02d}".format(point_to_line(p, q, M)), end="")
    print()

结果：

 48 42 39 36 32 28 27
 47 41 37 33 29 26 21
 46 40 35 30 25 20 18
 45 38 31 24 19 16 15
 44 34 23 17 14 12 11
 43 22 13 10 09 08 07
 00 01 02 03 04 05 06

• Why do you describe this as only a partial answer? It seems like it's the kind of underlying algorithm OP asked for. Am I missing something?

这不符合您对“简单”功能的要求，也不符合“合理高效”的要求。但原则上它会起作用，并且它可能会让人了解问题的难度。为了简单起见，让我们只考虑 0 < y ? x，因为完整的问题可以通过将完整的 2D 平面分成八个八分圆并以基本相同的方式将每个八分圆映射到其自己的整数范围来解决。

A point (x₁, y₁) is "anticlockwise" of (x₂, y₂) if and only if the slope y₁/x₁ is greater than the slope y₂/x₂. To map the slopes to integers in an order-preserving way, we can consider the sequence of all distinct fractions whose numerators and denominators are within range (i.e. up to 2³¹), in ascending numerical order. Note that each fraction's numerical value is between 0 and 1 since we are just considering one octant of the plane.

This sequence of fractions is finite, so each fraction has an index at which it occurs in the sequence; so to map a point (x, y) to an integer, first reduce the fraction y/x to its simplest form (e.g. using Euclid's algorithm to find the GCD to divide by), then compute that fraction's index in the sequence.

It turns out this sequence is called a Farey sequence; specifically, it's the Farey sequence of order 2³¹. Unfortunately, computing the index of a given fraction in this sequence turns out to be neither easy nor reasonably efficient. According to the paper
Computing Order Statistics in the Farey Sequence by Corina E. P?tra?cu and Mihai P?tra?cu, there is a somewhat complicated algorithm to compute the rank (i.e. index) of a fraction in O(n) time, where n in your case is 2³¹, and there is unlikely to be an algorithm in time polynomial in log n because the algorithm can be used to factorise integers.

All of that said, there might be a much easier solution to your problem, because I've started from the assumption of wanting to map these fractions to integers as densely as possible (i.e. no "unused" integers in the target range), whereas in your question you wrote that the number of distinct fractions is about 60% of the available range of size 2⁶⁴. Intuitively, that amount of leeway doesn't seem like a lot to me, so I think the problem is probably quite difficult and you may need to settle for a solution that uses a larger output range, or a smaller input range. At the very least, by writing this answer I might save somebody else the effort of investigating whether this approach is feasible.