112 lines
3.4 KiB
Python
112 lines
3.4 KiB
Python
from lib_misc import gcd
|
|
from lib_misc import is_permutation
|
|
|
|
|
|
class Primes(object):
|
|
def __init__(self, n_max):
|
|
import bitarray
|
|
self.n_max = n_max
|
|
b = bitarray.bitarray(n_max)
|
|
b.setall(True)
|
|
n = 1
|
|
b[n - 1] = False
|
|
while n * n <= n_max:
|
|
if b[n - 1] is True:
|
|
for i in range(n + n, n_max + 1, n):
|
|
b[i - 1] = False
|
|
n += 1
|
|
self.b = b
|
|
|
|
def iter_down(self):
|
|
for i in range(self.n_max, 0, -1):
|
|
if self.b[i - 1]:
|
|
yield i
|
|
raise StopIteration
|
|
|
|
def iter_up(self):
|
|
for i in range(1, self.n_max + 1):
|
|
if self.b[i - 1]:
|
|
yield i
|
|
raise StopIteration
|
|
|
|
def iter_range(self, n_min, n_max):
|
|
for i in range(n_min, n_max + 1):
|
|
if self.b[i - 1]:
|
|
yield i
|
|
raise StopIteration
|
|
|
|
def is_prime(self, n):
|
|
if n > self.n_max:
|
|
raise Exception("n greater than n_max")
|
|
return self.b[n - 1]
|
|
|
|
|
|
def relative_primes_count_naiv(n):
|
|
return len([i for i in range(1, n) if gcd(n, i) == 1])
|
|
|
|
|
|
def relative_primes_count_factors(n, fs):
|
|
"""
|
|
XXX: this method has a bug that first occurs for n = 60. Use totient
|
|
function from e072.py instead! For some curious reason it is good enough
|
|
for this problem. Probably because we only deal with two factors ever.
|
|
"""
|
|
from itertools import combinations
|
|
rel_primes_count = n - 1 # n itself is not a relative prime
|
|
for f in fs:
|
|
rel_primes_count -= (n // f - 1)
|
|
for f_1, f_2 in combinations(fs, 2):
|
|
f = f_1 * f_2
|
|
rel_primes_count += (n // f - 1)
|
|
return rel_primes_count
|
|
|
|
|
|
def get_phi(n):
|
|
r = relative_primes_count_naiv(n)
|
|
return n / r
|
|
|
|
|
|
def get_phi_factors(n, fs):
|
|
r = relative_primes_count_factors(n, fs)
|
|
return n / r
|
|
|
|
|
|
def euler_070():
|
|
"""
|
|
I struggled harder than I should have with this problem. I realized
|
|
quickly that a prime can't be the solution because a prime minus one
|
|
cannot be a permutation of itself. I then figured that the solution is
|
|
probably a number with two prime factors. I implemented an algorithm, but
|
|
it did not yield the right solution.
|
|
|
|
I tried a couple of things like squaring one prime factor which does not
|
|
yield a solution at all and three factors. Finally, I came up with a
|
|
faster algorithm to get the number of relative primes faster. With that
|
|
procedure I was then able to bruteforce the problem in 10 minutes.
|
|
|
|
When analyzing the solution I saw that it actually consists of two primes
|
|
which means my orginal algorithm had a bug. After reimplenting it was able
|
|
to find the solution in under 30 seconds. We could further optimize this
|
|
by making the search range for the two factors smaller. """
|
|
n = 10**7
|
|
ps = Primes(n // 1000)
|
|
phi_min = 1000
|
|
n_phi_min = 0
|
|
|
|
for p_1 in ps.iter_down():
|
|
for p_2 in ps.iter_range(1000, n // 1000):
|
|
n_new = p_1 * p_2
|
|
if n_new > n:
|
|
break
|
|
rel_primes_n_new = relative_primes_count_factors(n_new, [p_1, p_2])
|
|
phi_new = n_new / rel_primes_n_new
|
|
if phi_new < phi_min and is_permutation(n_new, rel_primes_n_new):
|
|
phi_min = phi_new
|
|
n_phi_min = n_new
|
|
return n_phi_min
|
|
|
|
|
|
if __name__ == "__main__":
|
|
print("e070.py: " + str(euler_070()))
|
|
assert(euler_070() == 8319823)
|