Solve problem 70 in Python.

python/e070.py

@@ -0,0 +1,106 @@
+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):
+    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)