5
\$\begingroup\$

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")
\$\endgroup\$
2
  • \$\begingroup\$ Note that your rephrasing of the problem excludes the information that the sum is 12523 if we set the limit at 1000 instead of 120000. This information is extremely useful to test your code! \$\endgroup\$
    – Stef
    Commented Jun 21 at 15:16
  • 1
    \$\begingroup\$ Sharing solutions to Project Euler problems beyond the first 100 goes against the terms of their website: projecteuler.net/about \$\endgroup\$
    – ViHdzP
    Commented Jun 22 at 3:15

3 Answers 3

3
\$\begingroup\$

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
\$\endgroup\$
8
\$\begingroup\$

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.

\$\endgroup\$
0
1
\$\begingroup\$

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
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.