Problem 127: abc-hits

Project Euler Problem 127 Solution

EulerSolve provides an optimized solution for Project Euler Problem 127, abc-hits, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary An abc-hit is a triple $(a,b,c)$ such that $a \lt b$, $a+b=c$, $\gcd(a,b)=1$, and $\operatorname{rad}(abc) \lt c$. Here $\operatorname{rad}(n)$ means the product of the distinct prime divisors of $n$. The task is to add all values of $c$ for abc-hits with $c \lt 120{,}000$. A naive search would try every $a$ for every $c$, which is far too much work. The implementations succeed because they exploit the exact multiplicative structure of $\operatorname{rad}(abc)$, and they scan candidate values in increasing order of radical so that most candidates are rejected very early. Mathematical Approach The key point is that once $c$ is fixed, the whole problem is controlled by the radicals of $a$, $b=c-a$, and $c$. That turns the search from a generic triple loop into a heavily pruned one-parameter scan. Radicals and Pairwise Coprimality Write $$\operatorname{rad}(n)=\prod_{p\mid n} p.$$ Because $c=a+b$, the condition $\gcd(a,b)=1$ automatically forces the other two gcds to be 1 as well: $$\gcd(a,c)=\gcd(a,a+b)=\gcd(a,b)=1,\qquad \gcd(b,c)=\gcd(b,a+b)=\gcd(a,b)=1.$$ So in every abc-hit the three numbers are pairwise coprime, and therefore $$\operatorname{rad}(abc)=\operatorname{rad}(a)\operatorname{rad}(b)\operatorname{rad}(c).$$ This identity is what makes the pruning in the code mathematically sound....

Detailed mathematical approach

Problem Summary

An abc-hit is a triple $(a,b,c)$ such that $a \lt b$, $a+b=c$, $\gcd(a,b)=1$, and $\operatorname{rad}(abc) \lt c$. Here $\operatorname{rad}(n)$ means the product of the distinct prime divisors of $n$. The task is to add all values of $c$ for abc-hits with $c \lt 120{,}000$.

A naive search would try every $a$ for every $c$, which is far too much work. The implementations succeed because they exploit the exact multiplicative structure of $\operatorname{rad}(abc)$, and they scan candidate values in increasing order of radical so that most candidates are rejected very early.

Mathematical Approach

The key point is that once $c$ is fixed, the whole problem is controlled by the radicals of $a$, $b=c-a$, and $c$. That turns the search from a generic triple loop into a heavily pruned one-parameter scan.

Radicals and Pairwise Coprimality

Write

$$\operatorname{rad}(n)=\prod_{p\mid n} p.$$

Because $c=a+b$, the condition $\gcd(a,b)=1$ automatically forces the other two gcds to be 1 as well:

$$\gcd(a,c)=\gcd(a,a+b)=\gcd(a,b)=1,\qquad \gcd(b,c)=\gcd(b,a+b)=\gcd(a,b)=1.$$

So in every abc-hit the three numbers are pairwise coprime, and therefore

$$\operatorname{rad}(abc)=\operatorname{rad}(a)\operatorname{rad}(b)\operatorname{rad}(c).$$

This identity is what makes the pruning in the code mathematically sound. The test is not approximate: once coprimality is known, the radical of the product is exactly the product of the radicals.

Reducing the Search to One Free Variable

For a fixed $c$, the equation $a+b=c$ determines $b$ from $a$:

$$b=c-a.$$

The inequality $a \lt b$ is equivalent to $a \lt c/2$, so only values

$$1 \le a \lt \frac{c}{2}$$

need to be checked. For each such $a$, we test the two conditions

$$\gcd(a,c-a)=1,\qquad \operatorname{rad}(a)\operatorname{rad}(c-a)\operatorname{rad}(c) \lt c.$$

If $\operatorname{rad}(c)\ge c$, then the radical inequality is already impossible because $\operatorname{rad}(a)\ge 1$ and $\operatorname{rad}(b)\ge 1$. That is why many values of $c$ terminate immediately.

The Radical-Order Pruning Bound

Let $r_c=\operatorname{rad}(c)$. Any abc-hit must satisfy

$$\operatorname{rad}(a)\operatorname{rad}(b)r_c \lt c.$$

In particular, $\operatorname{rad}(a)$ must be small. The implementations sort all positive integers below the limit by $\operatorname{rad}(n)$ and, for each fixed $c$, scan that list only until the radical gets too large.

The cutoff used in the code is the integer quotient

$$q=\left\lfloor \frac{c}{r_c}\right\rfloor.$$

Once $\operatorname{rad}(a)\ge q$, the scan can stop. If $\operatorname{rad}(a) \gt q$, then already $\operatorname{rad}(a)r_c \gt c$. If $\operatorname{rad}(a)=q$, then either $qr_c=c$, which fails the strict inequality immediately, or $qr_c \lt c \lt (q+1)r_c$. In that second case $b \gt 1$, so $\operatorname{rad}(b)\ge 2$, and therefore

$$\operatorname{rad}(a)\operatorname{rad}(b)r_c \ge 2qr_c \ge c.$$

So the radical-sorted list really does allow an early break, not merely a heuristic skip. The extra check $a \lt c/2$ is still needed, because the list is sorted by radical, not by numeric value.

Worked Example: $c=9$

Take $c=9$. Then $\operatorname{rad}(9)=3$, so the radical cutoff is

$$q=\left\lfloor \frac{9}{3}\right\rfloor=3.$$

Only candidates with $\operatorname{rad}(a) \lt 3$ can matter, and $a \lt 9/2$ leaves $a=1,2,4$.

For $a=1$, we get $b=8$, and

$$\operatorname{rad}(1)\operatorname{rad}(8)\operatorname{rad}(9)=1\cdot 2\cdot 3=6 \lt 9,$$

with $\gcd(1,8)=1$. So $(1,8,9)$ is an abc-hit.

For $a=2$, $b=7$, and the radical product is $2\cdot 7\cdot 3=42$, which fails. For $a=4$, $b=5$, and the product is $2\cdot 5\cdot 3=30$, which also fails. This is exactly the pattern used at scale: a cheap radical cutoff finds a short candidate list, then the full radical product and gcd test decide the survivors.

How the Code Works

The C++, Python, and Java implementations all follow the same two-stage algorithm: first build radicals efficiently, then enumerate $c$ and scan only the low-radical values of $a$.

Precomputing the Radicals

The implementation starts with an array filled with 1. It then runs a prime sieve. Whenever a prime $p$ is discovered, every multiple of $p$ has its radical multiplied by $p$ exactly once. After this pass, each entry stores the product of the distinct prime factors of its index.

At the same time, the program prepares a collection of all integers $1,2,\dots,119{,}999$ arranged in increasing order of $\operatorname{rad}(n)$, with the numeric value used as a tie-breaker. That ordering is the backbone of the later pruning.

Scanning Each $c$ Efficiently

For each $c$ from 3 up to $119{,}999$, the implementation computes $r_c=\operatorname{rad}(c)$ and the quotient $q=\lfloor c/r_c\rfloor$. It then walks through the pre-sorted candidate list and stops as soon as the candidate radical reaches $q$.

Each candidate still has to pass three exact filters: it must satisfy $a \lt c/2$, it defines $b=c-a$, and then the code checks both

$$\operatorname{rad}(a)\operatorname{rad}(b)\operatorname{rad}(c) \lt c$$

and

$$\gcd(a,b)=1.$$

Whenever both tests succeed, the current $c$ is added to the answer. One implementation also includes a small checkpoint: for the reduced bound $c \lt 1000$, the partial sum is $12523$.

Complexity Analysis

Let $N=120{,}000$. The sieve that computes all radicals costs $O(N\log\log N)$, and sorting the integers by radical costs $O(N\log N)$.

The dominant work is the search over candidates. In the worst case, the radical cutoff may still leave many entries, so the main loop is $O(N^2)$ candidate visits, with an additional $O(\log N)$ factor for the gcd test if one writes that cost explicitly. In practice it runs much faster than a full quadratic scan because most values of $c$ have a small quotient $c/\operatorname{rad}(c)$, so the radical-ordered loop breaks early. Memory usage is $O(N)$.

Footnotes and References

Problem page: Project Euler 127 - abc-hits
Radical of an integer: Wikipedia - Radical of an integer
Coprime integers: Wikipedia - Coprime integers
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
ABC conjecture background: Wikipedia - ABC conjecture

Problem 127 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <utility>
#include <vector>

namespace {

using i64 = std::int64_t;

struct Options {
    int limit = 120000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 3;
}

i64 solve(const int limit) {
    std::vector<int> rad(static_cast<std::size_t>(limit + 1), 1);
    std::vector<bool> is_prime(static_cast<std::size_t>(limit + 1), true);
    is_prime[0] = false;
    is_prime[1] = false;

    for (int p = 2; p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int m = p; m <= limit; m += p) {
            rad[static_cast<std::size_t>(m)] *= p;
            if (m > p) {
                is_prime[static_cast<std::size_t>(m)] = false;
            }
        }
    }

    std::vector<std::pair<int, int>> by_rad;
    by_rad.reserve(static_cast<std::size_t>(limit));
    for (int n = 1; n < limit; ++n) {
        by_rad.push_back({rad[static_cast<std::size_t>(n)], n});
    }
    std::sort(by_rad.begin(), by_rad.end());

    i64 sum_c = 0;

    for (int c = 3; c < limit; ++c) {
        const int rc = rad[static_cast<std::size_t>(c)];
        const int threshold = c / rc;

        for (const auto& [ra, a] : by_rad) {
            if (ra >= threshold) {
                break;
            }
            if (a >= c / 2) {
                continue;
            }

            const int b = c - a;
            const i64 rad_product = static_cast<i64>(ra) *
                                    static_cast<i64>(rad[static_cast<std::size_t>(b)]) *
                                    static_cast<i64>(rc);
            if (rad_product >= c) {
                continue;
            }
            if (std::gcd(a, b) != 1) {
                continue;
            }
            sum_c += c;
        }
    }

    return sum_c;
}

bool run_checkpoints() {
    if (solve(1000) != 12523) {
        std::cerr << "Checkpoint failed for limit=1000" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

# Problem 127: abc-hits
# Sum c for abc-hits: a+b=c, gcd(a,b)=1, rad(a*b*c) < c.

from math import gcd

def solve():
    limit = 120000
    rad = [1] * limit
    is_prime = [True] * limit
    for p in range(2, limit):
        if not is_prime[p]: continue
        for m in range(p, limit, p):
            rad[m] *= p
            if m > p: is_prime[m] = False
    
    # Sort values by rad for efficient pruning
    by_rad = sorted(range(1, limit), key=lambda n: (rad[n], n))
    
    total = 0
    for c in range(3, limit):
        rc = rad[c]
        threshold = c // rc
        for a in by_rad:
            if rad[a] >= threshold: break
            if a >= c // 2 + 1: continue
            b = c - a
            if rad[a] * rad[b] * rc >= c: continue
            if gcd(a, b) != 1: continue
            total += c
    
    print(total)

solve()

Java

import java.util.*;

public class Euler127 {
    public static void main(String[] args) {
        int limit = 120000;
        int[] rad = new int[limit];
        Arrays.fill(rad, 1);
        boolean[] isPrime = new boolean[limit];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int p = 2; p < limit; p++) {
            if (!isPrime[p])
                continue;
            for (int m = p; m < limit; m += p) {
                rad[m] *= p;
                if (m > p)
                    isPrime[m] = false;
            }
        }
        int[][] byRad = new int[limit - 1][2];
        for (int i = 0; i < limit - 1; i++) {
            byRad[i][0] = rad[i + 1];
            byRad[i][1] = i + 1;
        }
        Arrays.sort(byRad, (a, b) -> a[0] != b[0] ? a[0] - b[0] : a[1] - b[1]);
        long total = 0;
        for (int c = 3; c < limit; c++) {
            int rc = rad[c];
            int th = c / rc;
            for (int[] ar : byRad) {
                if (ar[0] >= th)
                    break;
                int a = ar[1];
                if (a >= (c + 1) / 2)
                    continue;
                int b = c - a;
                if ((long) ar[0] * rad[b] * rc >= c)
                    continue;
                if (gcd(a, b) != 1)
                    continue;
                total += c;
            }
        }
        System.out.println(total);
    }

    static int gcd(int a, int b) {
        while (b != 0) {
            int t = b;
            b = a % b;
            a = t;
        }
        return a;
    }
}