euler/python/e070.py

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)