from lib_prime import primes from functools import lru_cache def relative_primes_count_factors(n, factors): """ https://math.stackexchange.com/questions/1178847/relative-primes """ d = 1 for f in factors: n *= (f - 1) d *= f return n // d def prime_factors_unique(n, primes_list, primes_set): if n in primes_set: return [n] fs = [] rest = n for p in primes_list: if rest == 1: return fs if rest % p == 0: fs.append(p) while rest % p == 0: rest = rest // p def euler_072_brutal_brute_force(): d_max = 10**6 primes_list = primes(10**6) primes_set = set(primes_list) s = 0 for n in range(2, d_max + 1): factors = prime_factors_unique(n, primes_list, primes_set) rel_primes_count = relative_primes_count_factors(n, factors) s += rel_primes_count if n % 100000 == 0: print(n) # progress... haha. print(s) return s primes_list = primes(10**6) primes_set = set(primes_list) def euler_totient(p, a): return (p - 1) * p ** (a - 1) @lru_cache(maxsize=10**6) def relative_primes(n): if n == 1: return 1 if n in primes_set: return euler_totient(n, 1) for p in primes_list: if n % p == 0: a = 0 while n % p == 0: n = n // p a += 1 return euler_totient(p, a) * relative_primes(n) def euler_072(): s = 0 d_max = 10**6 for n in range(2, d_max + 1): rel_primes_count = relative_primes(n) s += rel_primes_count return s if __name__ == "__main__": print("e072.py: " + str(euler_072())) assert(euler_072() == 303963552391)