diff --git a/project_euler/problem_127/__init__.py b/project_euler/problem_127/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/project_euler/problem_127/sol1.py b/project_euler/problem_127/sol1.py new file mode 100644 index 000000000000..17e888438355 --- /dev/null +++ b/project_euler/problem_127/sol1.py @@ -0,0 +1,109 @@ +""" +Project Euler Problem 127: https://projecteuler.net/problem=127 + +abc-hits + +It takes about 10 minutes to run. +'Brute-force' solution that uses the following simplifications: +- if gcd(a, b) = 1 then gcd(a, c) = 1 and gcd(b, c) = 1 +- rad(a*b*c) = rad(a) * rad(b) * rad(c), for gcd(a, b) = 1 +- if a is even, b cannot b even for gcd(a, b) = 1 to be true. + +>>> solution(1000) +12523 +""" + +from numpy import sqrt + +N = 120000 + + +def generate_primes(n: int) -> list[bool]: + """ + Generates primes boolean array up to n. + + >>> generate_primes(2) + [False, False, True] + >>> generate_primes(5) + [False, False, True, True, False, True] + """ + primes = [True] * (n + 1) + primes[0] = primes[1] = False + for i in range(2, int(sqrt(n + 1)) + 1): + if primes[i]: + j = i * i + while j <= n: + primes[j] = False + j += i + return primes + + +def rad(n: int, primes_list: list[int]) -> int: + """ + Calculated rad - product of unique prime factors for n, using prime numbers + list primes_list. + + >>> rad(1, [1]) + 1 + >>> rad(12, [2, 3]) + 6 + """ + f = 1 + for p in primes_list: + if p > n: + break + if n % p == 0: + f *= p + return f + + +def gcd(a: int, b: int) -> int: + """ + Calculates greatest common divisor of a and b. + + >>> gcd(1, 10) + 1 + >>> gcd(14, 48) + 2 + """ + while b: + a, b = b, a % b + return a + + +def solution(c_less: int = 120000) -> int: + """ + Calculates all primes, rads, and then loops over a, b checking the conditions. + + >>> solution(10) + 9 + >>> solution(100) + 316 + """ + primes_bool = generate_primes(c_less) + primes_list = [] + for i in range(2, len(primes_bool)): + if primes_bool[i]: + primes_list += [i] + + rads = [1] * (c_less + 1) + for i in range(c_less + 1): + rads[i] = rad(i, primes_list) + + sum_c = 0 + for a in range(1, c_less): + rad_a = rads[a] + if a % 2 == 1: + r = range(1, min(a, c_less - a)) + else: + r = range(1, min(a, c_less - a), 2) + for b in r: + c = a + b + if rad_a * rads[b] * rads[c] < c and gcd(rad_a, rads[b]) == 1: + sum_c += c + + return sum_c + + +if __name__ == "__main__": + print(f"{solution() = }")