一种用于互质数模连续范围的快速算法/公式

Question

在我的项目中，有一个问题。为了简化问题，这里将问题阐述。有两个正互质的整数：a和b，其中a < b。列出了从1到b-1的a的倍数，然后是对b的模运算。 a mod b，2*a mod b，3*a mod b，...，(b-1)*a mod b 现在，还有另一个整数，称为n (1 <= n < b)。通过列表中的前n个数字，我们必须找出比m（1 <= m < b）小的数字有多少个。这可以通过蛮力法完成，从而给出一个O(n)。

例如： a=6, b=13, n=8, m=6 列表如下： 6, 12, 5, 11, 4, 10, 3, 9, 2, 8, 1, 7 这是从1到12的数字的排列，因为如果我们包括另一个数字0，任何两个互质数的模运算都会产生数字的排列。如果我们取a = 2，b = 13，那么列表将是2, 4, 6, 8, 10, 12, 1, 3, 5, 7, 9, 11，这会产生一种模式。但是如果a和b非常大（在我的项目中它们可以达到10^20），那么我不知道如何推导出这样大的数字的模式。

现在回到示例，我们取列表中的前n = 8个数字，得到： 6, 12, 5, 11, 4, 10, 3, 9 将less-than运算符与m = 6应用，它给出了小于m的数字总数为3，如下所述： 0, 0, 1, 0, 1, 0, 1, 0 其中0表示不小于m，1表示小于m。

由于上述算法的时间复杂度为O(n)，这对于范围在[0, 10^20]内的情况是不可接受的。因此，社区能否给出提示/线索/技巧，使我能够达到O(log n)甚至更好的O(1)解决方案？

- Nonymous NT

Answer 1

（警告：我对乘数范围不是[0，n)有些紧张，所以我进行了调整。很容易进行补偿。）

我将用经过测试的Python代码进行草图，并实现一个时间复杂度为 O(log max {a, b})的算法。首先，这里是一些实用函数和一个简单的实现。

from fractions import gcd
from random import randrange


def coprime(a, b):
    return gcd(a, b) == 1


def floordiv(a, b):
    return a // b


def ceildiv(a, b):
    return floordiv(a + b - 1, b)


def count1(a, b, n, m):
    assert 1 <= a < b
    assert coprime(a, b)
    assert 0 <= n < b + 1
    assert 0 <= m < b + 1
    return sum(k * a % b < m for k in range(n))

现在，我们该如何加快速度呢？第一个改进是将乘数分成不相交的范围，使得在每个范围内，对应的a的倍数在两个b的倍数之间。知道最低值和最高值后，我们可以通过向上取整除法计算小于m的倍数的数量。

def count2(a, b, n, m):
    assert 1 <= a < b
    assert coprime(a, b)
    assert 0 <= n < b + 1
    assert 0 <= m < b + 1
    count = 0
    first = 0
    while 0 < n:
        count += min(ceildiv(m - first, a), n)
        k = ceildiv(b - first, a)
        n -= k
        first = first + k * a - b
    return count

这不够快。第二个改进是用递归调用替换大部分while循环。在下面的代码中，j 是“完整”迭代次数，意味着存在一种环绕方式。term3 负责剩余的迭代，使用类似于 count2 的逻辑。

每个完整迭代都会在阈值 m 下贡献 floor(m / a) 或 floor(m / a) + 1 个残留物。我们是否获得了 + 1 取决于该迭代的 first 值。 first 从 0 开始，并在 while 循环中的每次迭代中通过 a - (b % a) 模 a 更改。当它低于某个阈值时，我们会获得 + 1，并且可以通过递归调用计算此计数。

def count3(a, b, n, m):
    assert 1 <= a < b
    assert coprime(a, b)
    assert 0 <= n < b + 1
    assert 0 <= m < b + 1
    if 1 == a:
        return min(n, m)
    j = floordiv(n * a, b)
    term1 = j * floordiv(m, a)
    term2 = count3(a - b % a, a, j, m % a)
    last = n * a % b
    first = last % a
    term3 = min(ceildiv(m - first, a), (last - first) // a)
    return term1 + term2 + term3

运行时间可以类比于欧几里得求最大公约数算法。

下面是一些测试代码，以证明我的正确性声明。在测试性能之前请记得删除断言。

def test(p, f1, f2):
    assert 3 <= p
    for t in range(100):
        while True:
            b = randrange(2, p)
            a = randrange(1, b)
            if coprime(a, b):
                break
        for n in range(b + 1):
            for m in range(b + 1):
                args = (a, b, n, m)
                print(args)
                assert f1(*args) == f2(*args)


if __name__ == '__main__':
    test(25, count1, count2)
    test(25, count1, count3)

- David Eisenstat