Problem 127: abc-hits
View on Project EulerProject 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
The radical \(\operatorname{rad}(n)\) of a positive integer \(n\) is the product of its distinct prime factors; for example \(504=2^3\cdot3^2\cdot7\), so \(\operatorname{rad}(504)=2\cdot3\cdot7=42\). A triple of positive integers \((a,b,c)\) is an abc-hit if \(\gcd(a,b)=\gcd(a,c)=\gcd(b,c)=1\), \(a\lt b\), \(a+b=c\) and \(\operatorname{rad}(abc)\lt c\). The statement's example is \((5,27,32)\): the three numbers are pairwise coprime, \(5\lt27\), \(5+27=32\), and \(\operatorname{rad}(5\cdot27\cdot32)=\operatorname{rad}(4320)=30\lt32\). There are \(31\) abc-hits with \(c\lt1000\), and their values of \(c\) add up to \(12523\). We must find the sum of \(c\) over all abc-hits with \(c\lt120000\); a value of \(c\) that occurs in several hits is counted once for each hit. [1] The name refers to the abc conjecture of number theory. [3]
The whole solution in one sentence: in a hit the three numbers share no prime factor, so the radical of \(abc\) is the product of the three radicals, which forces \(\operatorname{rad}(a)\) to be smaller than \(c/\operatorname{rad}(c)\), and once a sieve has computed all radicals and all numbers have been sorted by radical, each \(c\) only has to examine the short beginning of that sorted list.
Mathematical Approach
We look at small hits first, then reduce the four conditions to two simple tests per candidate \(a\), find a cut-off that rejects almost all candidates at once, and arrange the search so that the cut-off is applied by simply stopping.
1. Radicals and the first hits
A radical forgets the exponents of the prime factorisation: \(\operatorname{rad}(8)=2\), \(\operatorname{rad}(9)=3\), \(\operatorname{rad}(12)=6\) and \(\operatorname{rad}(72)=6\), while \(\operatorname{rad}(1)=1\), because \(1\) has no prime factors and an empty product is \(1\). Hence \(\operatorname{rad}(n)\) divides \(n\), and \(\operatorname{rad}(n)=n\) exactly when \(n\) is squarefree, that is, divisible by no square other than \(1\). [2] [4]
The smallest candidate, \((1,2,3)\), is not a hit, since \(\operatorname{rad}(6)=6\gt3\). The first hit is \((1,8,9)\), with \(\operatorname{rad}(72)=6\lt9\). Testing every triple with \(c\lt100\) gives exactly six hits, \((1,8,9)\), \((5,27,32)\), \((1,48,49)\), \((1,63,64)\), \((1,80,81)\) and \((32,49,81)\), and \(9+32+49+64+81+81=316\). Each of them contains high prime powers, such as \(8=2^3\), \(27=3^3\), \(49=7^2\) and \(81=3^4\): a hit needs numbers whose prime factors are few and repeated, so that the radical is much smaller than the number itself. The value \(c=81\) occurs twice; in the same way the \(31\) hits below \(1000\) have only \(23\) different values of \(c\), and the checkpoint \(12523\) counts every hit. [6]
2. Hits never run out
For every \(n\ge1\) the triple \((1,\,64^n-1,\,64^n)\) is a hit. The number \(b=64^n-1\) is divisible by \(64-1=63=9\cdot7\), because \(x-1\) divides \(x^n-1\). Write \(b=9m\); the primes of \(9m\) are those of \(3m\), so \(\operatorname{rad}(b)\le3m=b/3\). Consecutive numbers are coprime, so the three numbers are pairwise coprime, and \(\operatorname{rad}(abc)=1\cdot\operatorname{rad}(b)\cdot2\le2b/3\lt c\). For \(n=1\) this is the hit \((1,63,64)\) with radical \(42\), and for \(n=2\) it is \((1,4095,4096)\) with radical \(2730\). [3] So there are infinitely many hits, and the bound on \(c\) is essential.
The abc conjecture measures a triple by its quality \(\log c/\log\operatorname{rad}(abc)\), which exceeds \(1\) exactly for hits, for instance about \(1.226\) for \((1,8,9)\) and \(1.019\) for \((5,27,32)\), and states that for every \(\varepsilon\gt0\) only finitely many triples have a quality above \(1+\varepsilon\). The highest quality known, about \(1.6299\), belongs to \(2+3^{10}\cdot109=23^5\), found by Eric Reyssat. [3] Hits are sparse: there are \(1\), \(6\), \(31\) and \(120\) of them with \(c\) below \(10\), \(100\), \(1000\) and \(10000\). [7]
3. One gcd is enough, and the radical splits
If \(a+b=c\) and \(\gcd(a,b)=1\), then also \(\gcd(a,c)=\gcd(b,c)=1\): a number that divides \(a\) and \(c\) divides \(c-a=b\), hence divides \(\gcd(a,b)=1\), and in the same way for \(b\). This is the rule \(\gcd(a+mb,b)=\gcd(a,b)\) with \(m=1\), [5] so the first condition reduces to \(\gcd(a,b)=1\).
When \(a\), \(b\) and \(c\) are pairwise coprime, no prime divides two of them, so the distinct primes of \(abc\) are those of \(a\), \(b\) and \(c\) without repetition, and
$$\operatorname{rad}(abc)=\operatorname{rad}(a)\,\operatorname{rad}(b)\,\operatorname{rad}(c).$$
Without coprimality the product of the radicals can be larger. For \((8,24,32)\) it is \(2\cdot6\cdot2=24\), while \(8\cdot24\cdot32=6144=2^{11}\cdot3\) has radical \(6\); this triple even satisfies \(\operatorname{rad}(abc)\lt c\), but it is only \((1,3,4)\) multiplied by \(8\) and not a hit. For a fixed \(c\), the number \(b=c-a\) is determined by \(a\), and \(a\lt b\) means \(a\lt c/2\). So the hits with a given \(c\) correspond exactly to the integers \(a\) with \(1\le a\lt c/2\), \(\gcd(a,c-a)=1\) and
$$\operatorname{rad}(a)\,\operatorname{rad}(c-a)\,\operatorname{rad}(c)\lt c.$$
Testing all these \(a\) for every \(c\lt120000\) would mean about \(3.6\cdot10^9\) pairs, almost all of them hopeless.
4. A cut-off for \(\operatorname{rad}(a)\)
Because \(\operatorname{rad}(c)\) divides \(c\), the quotient \(q(c)=c/\operatorname{rad}(c)\) is a whole number; for example \(q(32)=16\), \(q(81)=27\) and \(q(49)=7\). Claim: every hit satisfies \(\operatorname{rad}(a)\lt q(c)\). Indeed \(\operatorname{rad}(c-a)\ge1\), so \(\operatorname{rad}(a)\operatorname{rad}(c)\)\({}\le\operatorname{rad}(a)\operatorname{rad}(c-a)\operatorname{rad}(c)\)\({}\lt c\), and dividing by \(\operatorname{rad}(c)\) gives the claim.
If \(c\) is squarefree, then \(q(c)=1\) and no \(a\) is left; more directly, \(\operatorname{rad}(abc)\ge\operatorname{rad}(c)=c\). Over three fifths of all integers are squarefree, the density being \(6/\pi^2\approx0.608\), [4] and among the \(119997\) values \(3\le c\lt120000\) there are \(72946\) squarefree ones. If \(q(c)=2\), only \(a=1\) has a radical below \(2\). Larger cut-offs also leave few candidates, because most numbers have large radicals: only \(111\) numbers below \(120000\) have a radical below \(q(49)=7\), namely \(1\), the powers of \(2\), \(3\) and \(5\), and the numbers \(2^i3^j\) with \(i,j\ge1\).
5. Searching in order of radical
The cut-off concerns \(\operatorname{rad}(a)\), not \(a\), so the search uses one list of all \(n\) with \(1\le n\lt120000\), sorted by radical and, for equal radicals, by \(n\). It begins with \(1\); then \(2,4,8,\dots,65536\), the \(16\) numbers with radical \(2\); then \(3,9,\dots,59049\) with radical \(3\); then \(5,25,\dots,78125\) with radical \(5\); then \(6,12,18,24,36,\dots\) with radical \(6\). The list is built once and serves every \(c\).
For a given \(c\) we walk along the list from the start and stop at the first entry whose radical is at least \(q(c)\): every later entry has a radical at least as large, so by Section 4 it cannot belong to a hit. Entries with \(a\ge c/2\) are skipped, since the list is not sorted by size. For every other entry we form \(b=c-a\) and test the product condition and the gcd condition of Section 3. Each number occurs exactly once in the list, so every hit with this \(c\) is found exactly once, and by Section 3 a candidate that passes both tests is a hit.
6. Worked example: \(c=32\)
Here \(\operatorname{rad}(32)=2\) and \(q(32)=16\). The walk passes the \(244\) entries with a radical below \(16\) and stops at the next entry, \(17\); no number has radical \(16\), because a radical is always squarefree. Of the \(244\) entries, \(229\) are at least \(16\) and are skipped, for example \(16\), \(32\), \(27\) and \(25\). The \(15\) numbers \(a=1,\dots,15\) are tested, in the list order \(1\); \(2,4,8\); \(3,9\); \(5\); \(6,12\); \(7\); \(10\); \(11\); \(13\); \(14\); \(15\). The product \(\operatorname{rad}(a)\operatorname{rad}(32-a)\cdot2\) is below \(32\) only for \(a=8\), where it is \(2\cdot6\cdot2=24\), and for \(a=5\), where it is \(5\cdot3\cdot2=30\). The gcd test removes \(a=8\), since \(\gcd(8,24)=8\), so \((5,27,32)\) is the only hit with \(c=32\).
For \(c=81\) we have \(q(81)=27\). Four candidates pass the product test, \(a=1\), \(32\), \(9\) and \(27\) in list order, with products \(30\), \(42\), \(54\) and \(54\). The gcd test removes \(9\) and \(27\), since \(\gcd(9,72)=9\) and \(\gcd(27,54)=27\), and leaves the two hits \((1,80,81)\) and \((32,49,81)\).
7. All radicals from one sieve
The radicals come from the sieve of Eratosthenes with one extra multiplication. [8] Start with the value \(1\) for every \(n\) and go through \(p=2,3,4,\dots\) in order. A \(p\) that has not been marked is prime, since no smaller prime divides it; then every multiple \(m\) of \(p\) has its value multiplied by \(p\), and the multiples above \(p\) are marked as composite. Each \(m\) is reached once for each of its distinct prime factors, whatever the exponents, so in the end its value is \(\operatorname{rad}(m)\): \(12\) is multiplied by \(2\) and by \(3\) and ends with \(6\), and so does \(72=2^3\cdot3^2\). The work is \(\sum_{p\le N}\lfloor N/p\rfloor\) multiplications, which grows like \(N\ln\ln N\) because the sum of \(1/p\) over the primes up to \(N\) grows like \(\ln\ln N\). [8] For \(N=120000\) these are \(321652\) multiplications in the C++ program, whose sieve also covers the index \(N\) itself, and \(321649\) in Python and Java.
How the Code Works
The C++ program [9] keeps its settings in Options (limit \(120000\), checkpoints on). parse_arguments accepts --skip-checkpoints and --limit= followed by digits, which parse_int_after_prefix converts without checking for int overflow; the program exits with status \(1\) on any other argument or on a limit below \(3\). Unless --skip-checkpoints is given, main first calls run_checkpoints, which checks that solve(1000) returns the statement's \(12523\); if it does not, main returns status \(2\). Then it prints solve(limit).
solve runs the sieve of Section 7 on a vector rad of int values and a vector is_prime, both indexed up to limit. It fills by_rad with the pairs \((\operatorname{rad}(n),n)\) for every \(n\) from \(1\) to limit - 1 and sorts it with std::sort, which compares pairs by radical first and then by \(n\): this is the list of Section 5. For each \(c\) from \(3\) to limit - 1 it computes rc, the radical of \(c\), and threshold = c / rc, the exact quotient \(q(c)\), and walks through by_rad: break when the radical reaches threshold, continue when a >= c / 2, then the product of the three radicals and finally std::gcd(a, b). A candidate that passes both tests adds \(c\) to sum_c. The cheap product test comes first and rejects almost everything: at the full limit only \(23665\) of the \(9796595\) tested candidates reach the gcd. The product is formed in 64-bit integers (i64), and this is needed: it is below \(c\cdot\operatorname{rad}(b)\lt c^2\), but the largest product formed at the full limit is \(14282923830\), more than six times \(2^{31}\).
The Python and Java versions follow the same steps with the limit fixed at \(120000\) and without the checkpoint. Python sorts the numbers themselves with the key (rad[n], n) and uses math.gcd and unbounded integers. Java stores the pairs as rows of int[][] byRad, sorts them with a comparator on radical and number, forms the product after a cast to long, computes the gcd with its own Euclidean loop and accumulates in a long.
The three versions write \(a\lt c/2\) in three ways. Java skips \(a\ge\lceil c/2\rceil\), which is exactly the complement of \(a\lt c/2\). Python skips \(a\ge\lfloor c/2\rfloor+1\) and therefore also examines \(a=c/2\) for even \(c\); this candidate is never accepted, because if it passes the product test, the gcd test rejects it, as \(\gcd(c/2,c/2)=c/2\ge2\). C++ skips \(a\ge\lfloor c/2\rfloor\) and therefore never tests \(a=(c-1)/2\), \(b=a+1\), for odd \(c\). This loses nothing here: none of the triples \((k,k+1,2k+1)\) with \(2k+1\lt2\cdot10^8\), which include the \(59999\) omitted below the limit, is a hit. That check is finite, but hits of this form are not expected at all: such a hit would give the triple \((1,\,4k(k+1),\,(2k+1)^2)\), and since \(k(k+1)\) is even and \(k\), \(k+1\), \(2k+1\) are pairwise coprime, its radical is \(\operatorname{rad}(k)\operatorname{rad}(k+1)\operatorname{rad}(2k+1)\), which is below \(2k+1\), the square root of its \(c\). Its quality would exceed \(2\), while the highest known quality is about \(1.6299\) (Section 2).
Complexity and Verification
Let \(N=120000\). The sieve takes \(O(N\log\log N)\) time and the sort \(O(N\log N)\). The walk for \(c\) examines every entry with a radical below \(q(c)\) and the one entry that stops it. In the worst case a walk covers the whole list, so the simple bounds are \(O(N^2)\) list visits and \(O(N^2\log N)\) time, since one gcd takes \(O(\log N)\) steps of Euclid's algorithm (for \(x\ge y\), the remainder of \(x\) modulo \(y\) is below \(x/2\)). The cut-off keeps the walks short: for \(N=120000\) they examine \(15080384\) entries in total, \(72946\) walks stop at the very first entry, \(5163792\) entries are skipped as too large and \(9796595\) candidates reach the product test, while a direct search over all \(a\lt c/2\) would test \(3599880001\) pairs, about \(240\) times as many. Between \(N=60000\) and \(N=2000000\), doubling \(N\) multiplies the number of examined entries by about \(2.5\), a growth like \(N^{1.33}\) in this range; this is a measurement, not a proved bound. A sharper cut-off, \(\operatorname{rad}(a)\lt q(c)/2\) from \(\operatorname{rad}(b)\ge2\), would keep only \(8195016\) of these \(14960387\) entries; the programs use the simpler one. Memory is \(O(N)\), about \(1.4\) MB for the radicals and the sorted pairs in C++. On one desktop computer the C++ and Java programs took a fraction of a second and the Python program about one second.
The C++ checkpoint reproduces the statement's \(12523\). For this article an independent brute force tested every pair \(a\lt b\) with \(a+b=c\) literally, with all three gcds and the radical of \(abc\) as the product of the distinct primes dividing \(a\), \(b\) or \(c\), without the product rule of Section 3 and without any cut-off. It finds exactly the same hits as the fast method for every \(c\lt30000\), \(234\) hits whose values of \(c\) add up to \(2540437\). The C++ program agrees with the brute force for every limit from \(3\) to \(3000\); copies of the Python and Java programs with the limit changed give \(12523\), \(441085\) and \(2540437\) for the limits \(1000\), \(10000\) and \(30000\); the fast method finds \(1268\) hits below \(10^6\), as listed in OEIS A216370; [7] and all three programs print the same value for the full limit. These checks are not part of the published solutions.
Footnotes and References
The references below supply the background facts used above. The reduction to two tests per candidate, the cut-off \(\operatorname{rad}(a)\lt c/\operatorname{rad}(c)\), the operation counts and the remark on hits with \(b=a+1\) are derived explicitly in this article.
- Project Euler 127 — abc-hits. The official statement defines rad(n) as the product of the distinct prime factors of n (rad(504) = 42) and calls (a, b, c) an abc-hit if gcd(a, b) = gcd(a, c) = gcd(b, c) = 1, a < b, a + b = c and rad(abc) < c. It gives the example (5, 27, 32) with rad(4320) = 30, states that the 31 abc-hits with c < 1000 have Σc = 12523, and asks for Σc over all abc-hits with c < 120000.
- Wikipedia — Radical of an integer. The radical of n is the product of the distinct primes dividing n, for example rad(504) = 42; it is the largest squarefree divisor of n, and it equals n when n is squarefree.
- Wikipedia — abc conjecture. For every ε > 0 only finitely many triples of coprime positive integers with a + b = c satisfy c > rad(abc)^(1+ε); the quality of a triple is log c / log rad(abc). The triples a = 1, b = 2^(6n) − 1, c = 2^(6n) satisfy c > rad(abc) because 9 divides b. Among the highest-quality triples listed there, the highest quality, about 1.6299, was found by Eric Reyssat for a = 2, b = 3^10 · 109 = 6436341, c = 23^5 = 6436343.
- Wikipedia — Square-free integer. An integer is squarefree if no square number other than 1 divides it, and it is squarefree if and only if it equals its radical. The squarefree numbers have natural density 6/π² ≈ 0.6079, so over three fifths of the integers are squarefree.
- Wikipedia — Greatest common divisor. Two numbers are coprime if their greatest common divisor is 1, and gcd(a + m·b, b) = gcd(a, b) for every integer m; with m = 1 this is the step from gcd(a, b) to gcd(c, b) used in Section 3.
- OEIS A130510 — ABC conjecture: values of c in the list of abc-hits. The entry calls a triple of positive integers with a + b = c, gcd(a, b) = 1 and c > rad(abc) an abc-hit. Its first 31 terms, 9, 32, 49, 64, 81, 81, 125, 128, …, 961, 968, are the values of c of the 31 hits below 1000, repeated values included.
- OEIS A216370 — Number of ABC triples with quality q > 1 and c < 10^n. The sequence begins 1, 6, 31, 120, 418, 1268 for n = 1, 2, …, 6; its comment defines these triples as the abc-hits.
- Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime as composite, starting with 2, and the next unmarked number is the next prime. It needs O(n log log n) operations because the prime harmonic series grows like log log n.
- C++ — Euler127.cpp. The linked immutable C++ revision contains
parse_arguments,solve,run_checkpointsand the checkpoint solve(1000) = 12523. It is the implementation record for this article.
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;
}
}