euler/python/e072.py

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)