diff --git a/python/e070.py b/python/e070.py index 644ca48..b07e842 100644 --- a/python/e070.py +++ b/python/e070.py @@ -46,6 +46,11 @@ def relative_primes_count_naiv(n): 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: diff --git a/python/e071.py b/python/e071.py index cbfc096..8af46a8 100644 --- a/python/e071.py +++ b/python/e071.py @@ -29,6 +29,22 @@ def euler_071(): return n_min +def euler_071_farey(): + """ + https://www.cut-the-knot.org/blue/Farey.shtml + n_min must be located between 2/5 and 3/7 for F_7 + To calculate F_d where d <= 10**6 we continously calculate the + median between 3/7 and the new fraction. + """ + n, d = (3, 7) + n_D, d_D = (2, 5) + while d_D + d <= 10**6: + n_D = n_D + n + d_D = d_D + d + return n_D + + if __name__ == "__main__": + print(euler_071_farey()) print("e071.py: " + str(euler_071())) assert(euler_071() == 428570) diff --git a/python/e072.py b/python/e072.py index e9e7cae..e7da668 100644 --- a/python/e072.py +++ b/python/e072.py @@ -1,8 +1,105 @@ +from collections import namedtuple +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 + + +Fraction = namedtuple("Fraction", ["n", "d"]) + + +# Procedures that were not actually helpful +def get_farey_series_length(f_min, f_max, d_max): + if f_min.d + f_max.d > d_max: + return 0 + f_med = Fraction(f_min.n + f_max.n, f_min.d + f_max.d) + l_min = get_farey_series_length(f_min, f_med, d_max) + l_max = get_farey_series_length(f_med, f_max, d_max) + return l_min + 1 + l_max + + +def get_farey_series(f_min, f_max, d_max): + if f_min.d + f_max.d > d_max: + return [] + f_med = Fraction(f_min.n + f_max.n, f_min.d + f_max.d) + l_min = get_farey_series(f_min, f_med, d_max) + l_max = get_farey_series(f_med, f_max, d_max) + return l_min + [f_med] + l_max + + +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(): - return 0 + 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() == 0) + assert(euler_072() == 303963552391)