优化列表理解以查找互质数对

Question

给定 A,B 打印 (a,b) 对的数量，使得 GCD(a,b)=1 和 1<=a<=A and 1<=b<=B.

这是我的回答：

return len([(x,y) for x in range(1,A+1) for y in range(1,B+1) if gcd(x,y) == 1])

我的答案对于小范围来说效果很好，但如果范围增加则需要足够的时间。 比如

• 1 <= A <= 10^5
• 1 <= B <= 10^5

有更好的写法或者可以优化吗？

Answer 1

由于您需要计算每个

pair

数字是否为gcd == 1，因此您应该预先计算所有素因子集。这样，您稍后可以通过检查两个数字素因子集的交集来快速确定两个数字是否互质。我们可以通过“类似筛子”的方法快速完成此操作。

factors = [set() for n in range(N)]
factored = collections.defaultdict(set)
for n in range(2, N):
    if not factors[n]:           # no factors yet -> n is prime
        for m in range(n, N, n): # all multiples of n up to N
            factors[m].add(n)
            factored[n].add(m)

在此之后，

factors[n]

将保存一组

n

（废话）的所有素因数，以及

factored[n]

所有由n分解的数字。这现在会派上用场，因为否则我们仍然需要检查最多 10,000 x 10,000 对数字，这在 Python 中仍然相当慢。但是结合使用 factors 和

factored

集，我们现在可以通过消除与

n

共享质因数的数字来快速找到给定数字的所有互质数。

for n in range(1, N):
    coprimes = set(range(1, N))  # start with all the numbers in the range
    for f in factors[n]:         # eliminate numbers that share a prime factor
        coprimes -= factored[f]
    print "%d is coprime with %r others" % (n, len(coprimes))

对于

N == 100

，结果对我来说似乎是合理的，对于 

N == 10000

，在我的计算机上大约需要 10 秒。这可能需要一些工作来适应您的实际问题，但我想这是一个好的开始。

根据

wikipedia

Answer 2

from collections import deque

def coprimes():
    tree = deque([[2, 1], [3, 1]])
    while True:
        m, n = tree.popleft()
        yield m, n
        tree.append([2 * m - n, m])
        tree.append([2 * m + n, m])
        tree.append([m + 2 * n, n])

这不是最快的算法，但却是最容易理解的。另请参阅

http://en.wikipedia.org/wiki/Farey_sequence

我知道这是一个老问题，但由于我需要类似的东西来解决另一个问题，并且无法足够快地找到东西，所以这是一个新的答案。需要

840ms

Answer 3

(x, y)

s.t.

x < y <= 10,000

。这个想法是首先获取素数列表（使用primesieve

）。然后，从该列表中，我们使用 算术基本定理

) 构造直到 

N 的所有整数。我在这里

提出了论点：我们只需要迭代所有可能的乘积

prod(p_i**k_i)，其中p_i是i^th

素数，

k_i

是任何非负整数。

下面的迭代器

gen_prod_primes()

可以方便地返回与每个整数一起使用的不同素数。因此，我们可以通过使用所有其他素数（在生成第一个整数时未使用的素数）来生成该整数的所有互素。把它们放在一起：

from primesieve import primes

def rem_p(q, p, p_distinct):
    q0 = q
    while q % p == 0:
        q //= p
    if q != q0:
        if p_distinct[-1] != p:
            raise ValueError(f'rem({q}, {p}, ...{p_distinct[-4:]}): p expected at end of p_distinct if q % p == 0')
        p_distinct = p_distinct[:-1]
    return q, p_distinct

def add_p(q, p, p_distinct):
    if len(p_distinct) == 0 or p_distinct[-1] != p:
        p_distinct += (p, )
    q *= p
    return q, p_distinct

def gen_prod_primes(p, upto=None):
    if upto is None:
        upto = p[-1]
    if upto >= p[-1]:
        p = p + [upto + 1]  # sentinel
    
    q = 1
    i = 0
    p_distinct = tuple()
    
    while True:
        while q * p[i] <= upto:
            i += 1
        while q * p[i] > upto:
            yield q, p_distinct
            if i <= 0:
                return
            q, p_distinct = rem_p(q, p[i], p_distinct)
            i -= 1
        q, p_distinct = add_p(q, p[i], p_distinct)

def coprimes(upto):
    p = list(primes(upto))
    ps = set(p)
    for y, p_distinct in gen_prod_primes(p, upto):
        for x, _ in gen_prod_primes(list(ps - set(p_distinct)), upto=y):
            yield x, y

现在：

>>> sorted(coprimes(6))
[(1, 1),
 (1, 2),
 (1, 3),
 (1, 4),
 (1, 5),
 (1, 6),
 (2, 3),
 (2, 5),
 (3, 4),
 (3, 5),
 (4, 5),
 (5, 6)]

%timeit list(coprimes(10_000))
840 ms ± 804 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

# but time grows fast:
%timeit list(coprimes(20_000))
6.57 s ± 9.97 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit list(coprimes(30_000))
13.1 s ± 23.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

>>> len(list(coprimes(10_000)))
1246336

>>> len(list(coprimes(20_000)))
9834865

>>> len(list(coprimes(30_000)))
19845763

备注：

迭代器产生无序互质。结果是所有互质数的一半；只需反转每个元组即可得到另一半。