How to quickly calculate the number of reduced own fractions?

Question

Currently, I am solving a mathematical problem in which I need to calculate the number of reduced proper fractions with a numerator and denominator of more than 1,000,000 (10 ^ 6).

I have code that works great for small numbers; the given example (value = 8) gives the correct (given) answer 21.

But this code seems very slow for large numbers for reasons I don't know. I read a bunch of similar questions about SO, but could not find anything useful. I looked at this one in more detail, but it really didn't help me. My code runs at an acceptable speed for values ​​up to 1000, then it becomes super-super-slow.

import math def is_prime(n): if n == 2: return True if n % 2 == 0 or n <= 1: return False sqr = int(math.sqrt(n)) + 1 for divisor in range(3, sqr, 2): if n % divisor == 0: return False return True def get_divisors(n): liste = [1] if n % 2 == 0: liste.append(2) for divisor in range(3, n+1): if n % divisor == 0: liste.append(divisor) return liste def intersect(a, b): return list(set(a) & set(b)) until = 1000 primes = list() for i in range(int(until)): if i != 1 and i != 0: if is_prime(i): primes.append(i) pos = 0 for i in range(1, until+1): if i%50 == 0: print(i) if is_prime(i): pos += (i-1) else: di = get_divisors(i) di.remove(1) for j in range(1, i): dj = get_divisors(j) if intersect(di, dj)==[]: pos+=1 print(pos) 

I want to know which parts of my program slow down and how to fix these problems.

Answer 1

We do not need simple factorization or inclusion-exception to solve this problem. We can modify the Euler sieve to dynamically calculate the Euler totalizer using its product formula, leading to the O(n) algorithm.

For each number in our accumulative list, we also store its totality. We use the fact that for the tuple (n, phi(n)) , where phi is the Euler function in magnitude:

  if m = p * n, for prime p: if p does not divide n then phi(m) = (p - 1) * phi(n) else phi(m) = p * phi(n) 

(Note that primes are generated as part of this algorithm by simply increasing the pointer to the next number not yet listed.)

For example n = 12 :

 list = [(2,1)] total = 1 prime: 2 2*2: (2 * 2, 2 * 1) total = total + 2 list: [(2,1), (4,2)] 2*2*2: (2 * 4, 2 * 2) total = total + 4 list: [(2,1), (4,2), (8,4)] prime: 3 total = total + 2 list: [(2,1), (3,2), (4,2), (8,4)] 3*3: (3 * 3, 3 * 2) total = total + 6 list: [(2,1), (3,2), (4,2), (8,4), (9,6)] 3*2: (3 * 2, 1 * (3-1)) total = total + 2 list: [(2,1), (3,2), (4,2), (6,2), (8,4), (9,6)] 3*4: (3 * 4, 2 * (3-1)) total = total + 4 list: [(2,1), (3,2), (4,2), (6,2), (8,4), (9,6), (12,4)] prime: 5 total = total + 4 list: [(2,1), (3,2), (4,2), (5,4), (6,2), (8,4), (9,6), (12,4)] 5*2: (5 * 2, 1 * (5-1)) total = total + 4 list: [(2,1), (3,2), (4,2), (5,4), (6,2), (8,4), (9,6), (10,4), (12,4)] prime: 7 total = total + 6 list: [(2,1), (3,2), (4,2), (5,4), (6,2), (7,6), (8,4), (9,6), (10,4), (12,4)] prime: 11 total = total + 10 list: [(2,1), (3,2), (4,2), (5,4), (6,2), (7,6), (8,4), (9,6), (10,4), (11,10), (12,4)] 

Python code (this code is actually an adaptation of the code from btilly, optimizing my idea):

 n = 8 total = 0 a = [None for i in range(0, n+1)] s = [] p = 1 while (p < n): p = p + 1 if a[p] is None: print("\n\nPrime: " + str(p)) a[p] = p - 1 total = total + a[p] s.append(p) limit = n / p new_s = [] for i in s: j = i while j <= limit: new_s.append(j) print j*p, a[j] a[j * p] = a[j] * p if (not j % p) else a[j] * (p-1) total = total + a[j * p] j = j * p s = new_s print("\n\nAnswer: " + str(total) + " => " + str([(k,d) for k,d in enumerate(a)])) 

Answer 2

If the prime numbers themselves are fast enough (I would recommend the Eratosthenes sieve, which is better than your usual primary testing for mass prime generation), divisor generation is not

 for divisor in range(3, n+1): if n % divisor == 0: liste.append(divisor) 

this loop has O(n) complexity when it can have O(n**0.5) complexity like this:

 liste = set() # okay, the var name isn't optimal now :) for divisor in range(3, int((n+1)**0.5)+1): if n % divisor == 0: other = n // divisor liste.add(divisor) liste.add(other) 

When you find the divider "low" (<= sqrt(n) ), other is just an extra "high" divider. This makes the cycle much shorter.

since you are going to perform the intersection by converting to set , why not create a set in the first place? The advantage is that in a perfect square you do not add the same value twice (you can check, but you are not sure that it will be faster, since the collision is rare)

now everything is set like this:

 if intersect(di, dj)==[]: pos+=1 

becomes:

 if not di.isdisjoint(dj): pos += 1 

• no set creation at each iteration
• no set to find out if the sets have common values. Just use set.isdisjoint

Last thing: As noted in the comments, when testing for primitiveness in the main loop, you call is_prime() again instead of storing your generated primes in set and testing with in .

Answer 3

For starters, stop doing simple factorization. Instead, use the Euclidean algorithm to check if the numerator and denominator has GCD 1.

If you really need prime numbers, find generation acceleration algorithms like the Sieve of Eratosthenes .

Finally, think about it: generate prime numbers (efficiently); factor, each denominator is how you are doing now. However, instead of going through all possible numerators, how can you generate real numerators for your needs?

Answer 4

A few simple optimizations are possible. You can try them and see if the performance sufficient to work in your use case improves.

Firstly, you can improve the way you search for primes to until . Instead of checking every number in the range, regardless of all numbers up to its square root, you could remember the prime numbers you have already found so far and only check them to the square root of the new candidate. If the number is odd, it is divided by some prime number less than or equal to its square root. Since the number of characters is less than n , less than n than n , as n gets larger, you can avoid many divisor checks by reusing the partial results that you already created at each step. This idea is similar to a sieve or Eratosthenes, not writing down all the numbers and not deleting them, but simply writing down those that are just the main ones when you go.

Secondly, when you check if fractions are in the least conditions, you do not need to look for any common divisors; just check for simple common factors. What for? if the fraction is not in the lower terms, there is a common denominator, which is divided by some prime number. Therefore, we can do without getting all the factors and focus only on getting a simple factorization. This can be accelerated only by checking the prime numbers that we received at the first stage as candidate factors, and when we find a match by dividing the found factor at each step. Thus, we can also avoid many checks. Something like that:

 p = 1 while number > 1 do if number % primes[p] == 0 then print "Another prime factor is " + primes[p] number = number / primes[p] else then p = p + 1 

The second of these optimizations is likely to give you the greatest benefit for your dollar.

Answer 5

Perform simple factorizations of all possible denominators. From each simple factorization, we calculate the denominator totator ( https://en.wikipedia.org/wiki/Euler%27s_totient_function ). This is the number of fractions with this denominator. Add all those and you're done.

Answer 6

This problem must be sought from a different angle. Let's look at the problem.

For the number n = a^k1 * b ^ k2 * c ^ k3 , where a, b, c ... is a prime number, the answer to question n - all number that is less than n and divides by either a or b or c ...

So, assuming we know that a, b, c ... how to calculate all number that is less than n and divides by a or b or c ? Let's look at the inclusion-exclusion theorem and this question. The number of positive integers in [1,1e18] that cannot be separated by any integers in [2,10] . Now we can have a reasonable algorithm, which over time will be O (2 ^ m), where m is the number of primes. Assuming that n <10 ^ 6, we see that m <8 (a very rough estimate: 2 * 3 * 5 * 7 * 11 * 13 * 17 * 23> 10 ^ 6)

Finally, now to get a list of primes for each number, we can use Sieve of Eratosthenes

Pseudocode:

 max = 1000001 // First, get all the prime number factors for each number list_prime = int[max][] for i in range (2, max): if list_prime[i] = []: list_prime[i].append(i) for j = i; j < max; j += i list_prime[j].append(i) // Second, calculate the number of value less than i and not share any factor with i result = 0 for i in range (2, max - 1): number_of_prime_factors = length(list_prime[i]) number_of_value = i for j = 1; j < 1 << number_of_prime_factors; j++: value = 0; for k in range (0, number_of_prime_factors): if ((1 << k) & j) != 0 : value += i / list_prime[i][k] if number of bit in j % 2 == 0: number_of_value -= value else: number_of_value += value result += number_of_value return result 

The time complexity should be O (n * 2 ^ m * m) with m <8 and n <= 10 ^ 6 </p>

Note With the Totient function, the time complexity should be further reduced to O (n * m), thanks @Matt Timmermans for his answer