Project Euler 127 - abc-hits

Question

Problem (Rephrased from here):

• The radical of \$n\$, \$rad(n)\$, is the product of distinct prime factors of \$n\$. For example, \$504 = 2^3 × 3^2 × 7\$, so \$rad(504) = 2 × 3 × 7 = 42\$.

We shall define the triplet of positive integers \$(a, b, c)\$ to be an abc-hit if:

• \$GCD(a, b) = GCD(a, c) = GCD(b, c) = 1\$
• \$a < b\$
• \$a + b = c\$
• \$rad(abc) < c\$

If \$(a,b,c)\$ is an abc-hit such that \$c<120000\$, Find sum of all \$c\$.

My code solves this problem in just under 2.5 minutes. Ideally I would like it to run in under a minute if possible. Originally I calculated the radicals using prime factorization, switching to a prime sieve reduced the running time by more than half. I can't (or don't know how to) identify any further bottlenecks that might slow down the program.

My code:

import math

def problem_one_hundred_and_twenty_seven() -> int:
    answer = 0
    limit = 120000
    # limit = 1000

    # Modified sieve of Eratosthenes
    rad = [0] + [1] * limit
    for i in range(2, len(rad)):
        if rad[i] == 1:
            for j in range(i, len(rad), i):
                rad[j] *= i

    for c in range(limit):
        for a in range(1, c // 2):
            if rad[a] * rad[c] >= c:
                continue
            b = c - a
            if b <= a:
                continue
            # If gcd(a, b) = 1 then gcd(c,b) = gcd(a+b,b) = gcd(a,b) = gcd(a,a+b) = gcd(a,c) = 1.
            # If a, b, c are all co-prime then rad(a * b * c) = rad(a) * rad(b) * rad(c).
            # If gcd(a, b) != 1 then they share a prime factor, => gcd(rad(a), rad(b)) != 1.
            # so gcd(rad(a), rad(b)) = 1 => gcd(a, b) = 1
            if math.gcd(a, b) == 1 and rad[a] * rad[b] * rad[c] < c:
                # print(a, b, c)
                answer += c

    return answer
    

if __name__ == "__main__":
    start_time = time.time()
    print(problem_one_hundred_and_twenty_seven())
    end_time = time.time()
    print(f"Found in {end_time - start_time} seconds")

Answer 1

Given that: \$GCD(a,b)=GCD(a,c)=GCD(b,c)=1\$

Since \$a\$, \$b\$ and \$c\$ have no common divisors then it also follows that: \$GCD(rad(a),rad(b))=GCD(rad(a),rad(c))=GCD(rad(b),rad(c))=1\$

And: \$rad(abc) = rad(a)*rad(b)*rad(c)\$.

---

Since \$GCD(a,b)=GCD(a,c)=GCD(b,c)=1\$ then when \$a\$, \$b\$ and \$c\$ are all even then \$GCD(a,b)>=2\$ so you can skip all the even \$a\$ values if \$c\$ is also even which eliminates a quarter of the checks that you need to do.

You can simplify the final part of the code to:

    for c in range(2, limit):
        rc = rad[c]
        if rc == c:
            continue
        step = 1 if c%2 == 1 else 2
        for a in range(1, (c + 1) // 2, step):
            ra = rad[a]
            rb = rad[c - a]
            if ra * rb * rc < c and math.gcd(ra, rb) == 1:
                # print(a, c - a, c)
                answer += c

Answer 2

Possible Bug

Whenever c is odd, for a in range(1, c // 2) misses an a, b combo.

For example: if c == 21, then for a in range(1, 21 // 2) becomes for a in range(1, 10), which means it never tries a == 10, b == 11.

Useless test

Since the range(...) enforces a < c // 2, because b = c - a, then b <= a will never be true, so that continue will never be executed.

Answer 3

Because c gets quite large, there is a small advantage to altering your short cut test to avoid a multiplication inside the loop on a (cost for this is a single divide outside). There are a lot of cutoffs on rad[a] is too large (and there might be even more if you apply this heuristic on b instead).

for c in range(limit):
     for a in range(1, c // 2):
         if rad[a] * rad[c] >= c:
             continue

Can be made marginally faster by changing it to the equivalent form (we know that c, rad[c] > 0)

for c in range(limit):
     alim = c // rad[c];
     for a in range(1, (c + 1) // 2):
         if rad[a]  >= alim:
             continue