Solved 72 in Python.

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:

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)

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)