CODEWITHSUNDEEP
pythonlisttupleslist-comprehensiongreatest-common-divisor
LearningNinja
LearningNinja

Reputation: 447

Optimizing list comprehension to find pairs of co-prime numbers

Given A,B print the number of pairs (a,b) such that GCD(a,b)=1 and 1<=a<=A and 1<=b<=B.

Here is my answer:

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

My answer works fine for small ranges but takes enough time if the range is increased. such as

is there a better way to write this or can this be optimized?

Upvotes: 3

Views: 5363

Answers (3)

Pierre D
Pierre D

Reputation: 26211

I know that's an oldie, but since I needed something like that for another problem, and couldn't find something fast enough, here is a new answer. It takes 840ms to list all (unordered) coprimes (x, y) s.t. x < y <= 10,000.

The idea is to first get a list of primes (using primesieve). Then, from that list, we construct all integers up to N using the Fundamental theorem of arithmetic). I presented the argument here: we just need to iterate through all the possible products prod(p_i**k_i) where p_i is the i^th prime number and k_i is any non-negative integer.

The iterator gen_prod_primes() below conveniently returns the distinct primes that were used along with each integer. So we can just generate all coprimes of that integer by using all the other primes (the ones not used in making the first integer).

Putting it all together:

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

Now:

>>> 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

Notes:

  • the iterator yields coprimes unordered.
  • the result is half all co-primes; just reverse each tuple to get the other half.

Upvotes: 1

tobias_k
tobias_k

Reputation: 82899

Since you need to calculate whether gcd == 1 for each pair of numbers, you should pre-calculate all the prime factor sets. This way, you can later very quickly determine whether two numbers are coprime by checking the intersection of their prime factor sets. We can do this fast in a sieve-like approach.

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)

After this, factors[n] will hold a set of all the prime factors of n (duh), and factored[n] all the numbers that are factored by n. This will come in handy now, since otherwise we would still need to check up to 10,000 x 10,000 pairs of numbers, which can still be rather slow in Python. But using the factors and factored sets in combination, we can now quickly find all the co-primes for a given number by eliminating the numbers that share a prime factor with 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))

For N == 100 the result looks plausible to me, and for N == 10000 it takes about 10 seconds on my computer. This might require some work to fit your actual problem, but I guess it's a good start.

Upvotes: 3

gog
gog

Reputation: 11347

According to wikipedia this should generate all coprimes:

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])

This is not the fastest algorithm, but the easiest to understand. See also http://en.wikipedia.org/wiki/Farey_sequence

Upvotes: 3

Related Questions