Problem 234: Semidivisible Numbers

Project Euler Problem 234 Solution

EulerSolve provides an optimized solution for Project Euler Problem 234, Semidivisible Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each integer $n$, let $p$ and $q$ be the consecutive primes satisfying $p \lt \sqrt n \lt q$. The number $n$ is called semidivisible when exactly one of those two primes divides it. The task is to compute $$S(L)=\sum_{\substack{n\le L\\ n\text{ semidivisible}}} n$$ for the given limit $L=999966663333$. A direct search over all $n\le L$ would be hopeless. The successful idea is to stop thinking about individual integers and instead group them by the pair of consecutive primes surrounding $\sqrt n$. Once that pair is fixed, the semidivisible condition becomes a clean inclusion-exclusion calculation on an interval. Mathematical Approach The three implementations all exploit the same structure: partition the line into prime-square bands, compute one closed-form contribution for each band, and add those contributions. Prime-square bands identify the relevant primes Suppose $p$ and $q$ are consecutive primes. Then the condition $$p \lt \sqrt n \lt q$$ is equivalent to $$p^2 \lt n \lt q^2.$$ So every semidivisible number belongs to exactly one interval between two consecutive prime squares. Inside that whole band, the lower surrounding prime is always $p$ and the upper surrounding prime is always $q$....

Detailed mathematical approach

Problem Summary

For an integer $n\ge4$, the lower prime square root $\operatorname{lps}(n)$ is the largest prime $\le\sqrt n$, and the upper prime square root $\operatorname{ups}(n)$ is the smallest prime $\ge\sqrt n$; for example $\operatorname{lps}(4)=2=\operatorname{ups}(4)$, $\operatorname{lps}(1000)=31$ and $\operatorname{ups}(1000)=37$. The number $n$ is semidivisible if one of $\operatorname{lps}(n)$ and $\operatorname{ups}(n)$ divides $n$, but not both. The semidivisible numbers up to $15$ are $8$, $10$ and $12$, with sum $30$; $15$ itself is not semidivisible, because it is a multiple of both $\operatorname{lps}(15)=3$ and $\operatorname{ups}(15)=5$; and the $92$ semidivisible numbers up to $1000$ have the sum $34825$. We must find the sum $S(L)$ of all semidivisible numbers $n\le L$ for $L=999\,966\,663\,333$. ^[1]

The whole solution in one sentence: between the squares of two consecutive primes $p\lt q$ the two prime square roots are always $p$ and $q$, so the semidivisible numbers there are the multiples of exactly one of $p$ and $q$, whose sum follows from three sums of arithmetic progressions, and adding these sums over all primes $p$ with $p^2\le L$ gives $S(L)$.

Mathematical Approach

Testing every $n\le L$ separately would take about $10^{12}$ tests. Instead we find out where the two prime square roots change, turn "one of them divides $n$, but not both" into sums of multiples, and evaluate those sums in closed form.

1. Small cases

Working out the two roots for $4\le n\le49$ shows that they stay the same between two consecutive prime squares:

for $n=4$ both roots are $2$, and for $5\le n\le8$ they are $2$ and $3$;
for $n=9$ both are $3$, and for $10\le n\le24$ they are $3$ and $5$;
for $n=25$ both are $5$, for $26\le n\le48$ they are $5$ and $7$, and for $n=49$ both are $7$.

For instance, $\sqrt{16}=4$ is not prime, so $\operatorname{lps}(16)=3$ and $\operatorname{ups}(16)=5$. The semidivisible numbers up to $49$ are $8$; $10, 12, 18, 20, 21, 24$; and $28, 30, 40, 42, 45$, the first terms of the sequence A157939 of the OEIS. ^[2] The numbers up to $49$ divisible by both roots are the prime squares $4, 9, 25, 49$, whose two roots coincide, and the products $6=2\cdot3$, $15=3\cdot5$ and $35=5\cdot7$ of their two roots. Sections 2 and 4 show that this is no accident.

2. Bands between consecutive prime squares

Let $p\lt q$ be consecutive primes, so that no prime lies strictly between them. We call the integers $n$ with $p^2\lt n\lt q^2$ the band of $p$; the bands of $2$, $3$ and $5$ run from $5$ to $8$, from $10$ to $24$ and from $26$ to $48$.

Claim. Every $n$ in the band of $p$ has $\operatorname{lps}(n)=p$ and $\operatorname{ups}(n)=q$.

Proof. Taking square roots keeps the order of positive numbers, so $p\lt\sqrt n\lt q$. A prime $r$ with $p\lt r\le\sqrt n$, or with $\sqrt n\le r\lt q$, would lie strictly between $p$ and $q$, which is impossible. Hence the largest prime $\le\sqrt n$ is $p$, and the smallest prime $\ge\sqrt n$ is $q$. $\square$

For $n=p^2$ both roots equal $p$, which divides $n$, so a prime square is never semidivisible, which is why the list of the statement lacks $4$ and $9$. Every $n\ge4$ is a prime square or lies in exactly one band: if $p$ is the largest prime with $p^2\le n$ and $q$ is the next prime, then $p^2\le n\lt q^2$.

3. One divisor but not both

Fix a band and a number $n$ in it. Put $a=1$ if $p\mid n$ and $a=0$ otherwise, and define $b$ in the same way for $q$. The four cases $(a,b)=(0,0)$, $(1,0)$, $(0,1)$, $(1,1)$ give $a+b-2ab=0$, $1$, $1$, $0$, so $a+b-2ab$ is $1$ exactly when $n$ is semidivisible and $0$ otherwise: it is the exclusive or of the two conditions, the indicator of a symmetric difference of two sets. ^[3] The product $ab$ is $1$ exactly when $pq\mid n$: if $n=pm$ and also $q\mid n$, then $q$ divides the product $pm$ but not the prime $p$, so $q\mid m$ by Euclid's lemma, ^[4] and $pq\mid n$. Multiplying by $n$ and adding over the numbers $A\le n\le B$ of the band gives its contribution

$$\begin{aligned}\mathcal B_p&=M_p(A,B)+M_q(A,B)\\&\quad-2M_{pq}(A,B),\end{aligned}$$

where $M_d(A,B)$ is the sum of the multiples of $d$ in the interval $[A,B]$. We take $A=p^2+1$ and $B=\min(L,q^2-1)$, so that only the last band can be cut short at $L$. By Section 2,

$$S(L)=\sum_{p^2\le L}\mathcal B_p,$$

where $p$ runs over the primes and $q$ is always the prime after $p$.

4. Only one multiple of $pq$ in a band

Bertrand's postulate says that for every integer $m\gt1$ some prime lies strictly between $m$ and $2m$. ^[5] With $m=p$ it gives $q\lt2p$, so $2pq\gt q^2$, while $pq\gt p^2$. Hence $pq$ is the only multiple of $pq$ in the band of $p$: $M_{pq}(A,B)$ is $pq$ for a complete band, and $pq$ or $0$ for a band cut at $L$. So the numbers $n\ge4$ divisible by both prime square roots are exactly the prime squares and the products $pq$ of consecutive primes, as seen in Section 1.

5. Sums of multiples

The multiples of $d$ in $[A,B]$ are the numbers $dk$ with $k_0\le k\le k_1$, where $\lceil x\rceil$ and $\lfloor x\rfloor$ round $x$ up and down to an integer and

$$k_0=\left\lceil\frac{A}{d}\right\rceil,\qquad k_1=\left\lfloor\frac{B}{d}\right\rfloor.$$

If $k_0\gt k_1$ there are none and $M_d(A,B)=0$. Otherwise they form an arithmetic progression with $c=k_1-k_0+1$ terms, from $f=dk_0$ to $\ell=dk_1$. Writing the sum once forwards and once backwards and adding term by term gives $c$ pairs, each with the sum $f+\ell$, ^[6] so

$$M_d(A,B)=\frac{c\,(f+\ell)}{2}.$$

The product $c\,(f+\ell)$ is twice a sum of integers, so it is always even. For $d=5$, $A=26$ and $B=48$ we get $k_0=6$, $k_1=9$, $c=4$, $f=30$ and $\ell=45$, so $M_5(26,48)=4\cdot75/2=150$, which is indeed $30+35+40+45$.

6. Worked example: the band of $5$

The band of $5$ runs from $A=26$ to $B=48$. Section 5 gives $M_5(26,48)=150$, $M_7(26,48)=28+35+42=105$ and $M_{35}(26,48)=35$, so

$$\mathcal B_5=150+105-2\cdot35=185.$$

The semidivisible numbers of this band are $28, 30, 40, 42, 45$ from Section 1, and indeed $28+30+40+42+45=185$. The term $-2M_{35}$ removes $35$, which both $M_5$ and $M_7$ contain.

7. Worked example: the two checkpoints

For $L=15$ only the bands of $2$ and $3$ are needed. The band of $2$ runs from $5$ to $8$, with $M_2=6+8=14$, $M_3=6$ and $M_6=6$, so $\mathcal B_2=14+6-12=8$. The band of $3$ is cut to the numbers from $10$ to $15$, with $M_3=12+15=27$, $M_5=10+15=25$ and $M_{15}=15$, so $\mathcal B_3=27+25-30=22$. Hence $S(15)=8+22=30$.

For $L=1000$ we have $31^2=961\le1000\lt37^2$, so the bands of the eleven primes from $2$ to $31$ are needed. The ten complete bands, of $p=2, 3, 5, 7, 11,$ $13, 17, 19, 23, 29$, contribute $8$, $105$, $185$, $1216$, $737$, $3211$, $1649$, $5757$, $14451$ and $4553$, coming from $1$, $6$, $5$, $14$, $5$, $14$, $5$, $13$, $21$ and $5$ semidivisible numbers. The band of $31$ is cut to the numbers from $962$ to $1000$. It contains the multiples $992=32\cdot31$, $962=26\cdot37$ and $999=27\cdot37$, while $31\cdot37=1147$ exceeds $1000$; so $M_{31}=992$, $M_{37}=1961$, $M_{1147}=0$ and $\mathcal B_{31}=2953$. Altogether $S(1000)=34825$ from $92$ numbers, the second checkpoint.

8. How many numbers a band holds

Let $g=q-p$. The multiples of $p$ in the band are the numbers $kp$ with $p+1\le k\le\lfloor q^2/p\rfloor$, because $p\nmid q^2$; since $q^2=p^2+2gp+g^2$, there are $2g+\lfloor g^2/p\rfloor$ of them. The multiples of $q$ are the numbers $kq$ with $\lfloor p^2/q\rfloor+1\le k\le q-1$, because $q\nmid p^2$; since $p^2=q^2-2gq+g^2$, there are $2g-1-\lfloor g^2/q\rfloor$ of them. Both lists contain $pq$, which is not semidivisible, so a complete band holds

$$4g-3+\left\lfloor\frac{g^2}{p}\right\rfloor-\left\lfloor\frac{g^2}{q}\right\rfloor$$

semidivisible numbers, simply $4g-3$ whenever $g^2\lt p$. Twin primes, with $g=2$ and $p\ge5$, give $5$ numbers, as for $p=5, 11, 17, 29$ in Section 7. Ignoring the floor terms and the cut last band, the gaps add up to roughly $\sqrt L$, which suggests about $4\sqrt L-3\pi(\sqrt L)$ semidivisible numbers up to $L$, where $\pi(x)$ counts the primes up to $x$. ^[7] The exact counts are $36\,310$ up to $10^8$ and $3\,764\,437$ up to $L=999\,966\,663\,333$, so even listing them one by one would be fast.

9. The bands for the given limit

For $L=999\,966\,663\,333$ we have $\lfloor\sqrt L\rfloor=999\,983$, because $999\,983^2=999\,966\,000\,289$${}\le L\lt999\,984^2$, and $999\,983$ is the largest prime below $10^6$. The primes with $p^2\le L$ are therefore the $\pi(10^6)=78\,498$ primes below $10^6$. ^[7] All their bands are complete except the last: the prime after $999\,983$ is $1\,000\,003$, so the band of $999\,983$ is cut at $L$ after only $L-999\,983^2=663\,044$ numbers. This piece holds no semidivisible number: the next multiple of $999\,983$ is $999\,983\cdot999\,984$${}=999\,967\,000\,272\gt L$, and the nearest multiples of $1\,000\,003$ are $999\,965\,999\,889\lt999\,983^2$ and $999\,966\,999\,892\gt L$.

10. Primes up to $\sqrt L$: the linear sieve

The bands need every prime up to $\sqrt L$ and the prime after the last of them. The programs find them with a linear sieve, a variant of the sieve of Eratosthenes that crosses out every composite number exactly once. ^[8] It runs through $i=2,3,\dots,N$. If $i$ has not been crossed out, it is a new prime. Then $i$ crosses out $i\cdot p$ for the primes $p=2,3,5,\dots$ found so far, up to and including the first $p$ that divides $i$, as long as $i\cdot p\le N$. For example, $i=12$ crosses out only $24$, because $2$ divides $12$; $36$ is left for $i=18$, and $i=15$ crosses out $30$ and $45$.

Claim. Every composite $m\le N$ is crossed out exactly once, as $i\cdot p$ with $p$ its smallest prime factor.

Proof. Write $m=i\cdot p$ with $p$ the smallest prime factor of $m$. All prime factors of $i$ are at least $p$, so $p\le i$ is already in the list when $i$ is handled; no smaller prime divides $i$, and $i\cdot p'\lt m\le N$ for every smaller prime $p'$, so the loop for $i$ reaches $p$ and crosses out $m$. Conversely, the loop for any $i'$ uses only primes $p'$ up to the smallest prime factor of $i'$, so $p'$ is the smallest prime factor of $i'p'$; hence a crossing-out $m=i'p'$ has $p'=p$ and $i'=i$. As $m$ is crossed out at the step $i\lt m$, only primes join the list. $\square$

The programs sieve up to $N=\lfloor\sqrt L\rfloor+110=1\,000\,093$ and find $78\,503$ primes. The prime after the last $p$ with $p^2\le L$ can lie beyond such a margin: below $10^6$ the largest gap between consecutive primes is $114$, from $492\,113$ to $492\,227$, ^[9] and for a limit just above $492\,113^2$ the sieve stops at $492\,223$. So the programs append the first prime above $N$, found by trial division, that is, by dividing it by the primes up to its square root. For the given limit it is $1\,000\,099$, and it is not used, as the loop stops at $1\,000\,003$.

How the Code Works

The C++ program ^[10] starts in main. parse_arguments accepts --limit=L, whose default is $999\,966\,663\,333$, and --skip-checkpoints; unknown arguments and limits below $4$ end the program with exit code $1$, and the digit parser does not check for 64-bit overflow. run_checkpoints requires solve_u64(15) to be $30$ and solve_u64(1000) to be $34825$, and otherwise ends the program with exit code $2$. Then the result of solve_i128 is printed by i128_to_string.

solve_i128 sets sqrt_limit to the integer part of a long double square root of the limit plus $10$ and calls primes_up_to(sqrt_limit + 100), the linear sieve of Section 10 with a byte array is_composite; the margin also absorbs any rounding error of the root. A while loop then tests candidate $=$ sqrt_limit + 101 and the following numbers by trial division and appends the first prime; for the given limit it tests $1\,000\,094$ to $1\,000\,099$.

The main loop takes p = primes[i] and q = primes[i + 1] and stops once $p^2$ exceeds the limit. It sets left $=p^2+1$ and right $=\min(L,q^2-1)$, the $A$ and $B$ of Section 3; left > right happens only when the limit is itself a prime square, such as $L=49$, and the empty band is skipped. sum_multiples_in_range(d, left, right) is Section 5: it forms the first multiple ((left + d - 1) / d) * d, the last multiple (right / d) * d, the count $c=(\ell-f)/d+1$ and the result $(f+\ell)\,c/2$, or $0$ when $f\gt\ell$. The loop adds sum_p + sum_q - 2 * sum_pq, the $\mathcal B_p$ of Section 3, to total; it computes $M_{pq}$ with the general function, although by Section 4 it is $pq$ or $0$.

Every value fits easily into a signed 64-bit integer, whose maximum is $2^{63}-1\approx9.2\cdot10^{18}$: the largest $(f+\ell)\,c$ is about $3.3\cdot10^{14}$, and the total is less than $L$ times the $3\,764\,437$ semidivisible numbers of Section 8, which is below $3.8\cdot10^{18}$. So Java's 64-bit long suffices; the C++ program still sums in __int128_t, the signed 128-bit type of GCC and Clang, and Python's integers are unbounded. ^[11] The Python and Java programs have the same three functions, named primesUpTo, sumMultiplesInRange and solve in Java, but no options or checkpoints; Python's solve takes the limit as an optional argument, Java's solve() fixes it, and both take the square root in double precision, which the margin makes harmless. Java halves $f+\ell$ when it is even and $c$ otherwise; this is exact, because $f+\ell=2f+(c-1)d$ is odd only when $c$ is even.

Complexity and Verification

The linear sieve takes $O(N)$ time with $N\approx\sqrt L$, since each composite is crossed out once and each $i$ adds at most one further test, and the main loop makes $\pi(\sqrt L)$ passes, $78\,498$ here, each with three closed-form sums: $O(\sqrt L)$ steps in all. The search for the extra prime tests the $h$ numbers from $N+1$ up to the next prime, with at most $\pi(\sqrt{N+h})$ divisions each, and $h$ is bounded by a gap between consecutive primes: here $h=6$, with $177$ divisions, and for every $L\le10^{12}$ there are at most $114$ candidates and $1148$ divisions. The memory is $O(\sqrt L)$: $N+1=1\,000\,094$ bytes of flags and the list of $78\,504$ primes, 4 bytes each in C++ and integer objects in Python and Java. On the test machine (clang -O2, OpenJDK 17, Python 3.11) the C++, Java and Python programs take about $0.004$, $0.05$ and $0.15$ seconds.

Only the C++ program checks itself, with the two checkpoints. For this article, a brute force that finds $\operatorname{lps}(n)$ and $\operatorname{ups}(n)$ for each $n$ from a table of primes, without bands, was compared with the published C++ program run with --limit: they agree for every $L$ from $4$ to $3000$, for $L=10^4,10^5,\dots,10^9$ and at $997^2$, $9973^2$ and their neighbours, and so does Python's solve(L). Run over all $n\le999\,966\,663\,333$, in about half an hour of compiled code, the brute force gives the number that all three programs print. At $L=492\,113^2+5\cdot492\,113$, where the band of $492\,113$ holds $10$ semidivisible numbers and needs the appended prime $492\,227$, the C++ program and Python's solve(L) agree with a separate walk through the multiples in each band. The count formula of Section 8 holds in all $78\,497$ complete bands. None of these checks is part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the band decomposition, the band formula, the band count and all numbers in the examples are derived or computed in this article.

Project Euler 234 — Semidivisible Numbers. The official statement defines lps(n), the largest prime ≤ √n, and ups(n), the smallest prime ≥ √n, for n ≥ 4, calls n semidivisible if one of them divides n but not both, gives lps(4) = 2 = ups(4), lps(1000) = 31, ups(1000) = 37 and the sums 30 of the semidivisible numbers up to 15 and 34825 of the 92 up to 1000, and asks for the sum of all semidivisible numbers not exceeding 999966663333.
OEIS A157939 — Numbers k divisible by prevprime(sqrt(k)) or nextprime(sqrt(k)) but not both. The sequence begins 8, 10, 12, 18, 20, 21, 24, 28, 30, 40, 42, 45, the semidivisible numbers of Section 1.
Wikipedia — Symmetric difference. The indicator function of the symmetric difference of two sets is the exclusive or, the sum modulo 2, of their indicator functions; for the values 0 and 1 this is a + b − 2ab.
Wikipedia — Euclid's lemma. If a prime divides the product of two integers, it divides at least one of them.
Wikipedia — Bertrand's postulate. For every n > 1 there is a prime p with n < p < 2n, which Chebyshev proved in 1852; equivalently, consecutive primes satisfy pₙ₊₁ < 2pₙ.
Wikipedia — Arithmetic progression. The sum of n terms is n/2 times the sum of the first and the last term; it is derived by writing the sum forwards and backwards and adding.
Wikipedia — Prime-counting function. π(x) is the number of primes up to x, and the table of values gives π(10⁶) = 78,498.
cp-algorithms — Linear Sieve. For each i the sieve marks i · p for the primes p up to the smallest prime factor of i, so every composite number is reached exactly once, through its smallest prime factor, in O(n) time.
Wikipedia — Prime gap. The table of maximal prime gaps lists the gap 114 after the prime 492113 as the largest one below 10⁶.
C++ — Euler234.cpp. The linked immutable revision contains the functions and the two checkpoints described above. It is the implementation record for this article.
GCC manual — 128-bit Integers documents the __int128 extension, Java — Primitive Data Types gives the range −2⁶³ to 2⁶³ − 1 of long, and Python — numeric types states that integers have unlimited precision.

Problem 234 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using u64 = std::uint64_t;
using i128 = __int128_t;

struct Options {
    u64 limit = 999966663333ULL;
    bool run_checkpoints = true;
};

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

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(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_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 4;
}

std::vector<int> primes_up_to(const int n) {
    std::vector<std::uint8_t> is_composite(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;

    for (int i = 2; i <= n; ++i) {
        if (!is_composite[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
        for (int p : primes) {
            const long long v = 1LL * i * p;
            if (v > n) {
                break;
            }
            is_composite[static_cast<std::size_t>(v)] = 1;
            if (i % p == 0) {
                break;
            }
        }
    }

    return primes;
}

i128 sum_multiples_in_range(const u64 d, const u64 left, const u64 right) {
    if (left > right) {
        return 0;
    }
    const u64 first = ((left + d - 1) / d) * d;
    const u64 last = (right / d) * d;
    if (first > last) {
        return 0;
    }

    const u64 count = (last - first) / d + 1;
    return (static_cast<i128>(first) + static_cast<i128>(last)) * static_cast<i128>(count) / 2;
}

std::string i128_to_string(i128 value) {
    if (value == 0) {
        return "0";
    }

    bool neg = value < 0;
    if (neg) {
        value = -value;
    }

    std::string out;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    if (neg) {
        out.push_back('-');
    }
    std::reverse(out.begin(), out.end());
    return out;
}

i128 solve_i128(const u64 limit) {
    const int sqrt_limit = static_cast<int>(std::sqrt(static_cast<long double>(limit))) + 10;
    std::vector<int> primes = primes_up_to(sqrt_limit + 100);

    // Ensure one prime above sqrt(limit) exists.
    int candidate = sqrt_limit + 101;
    while (true) {
        bool is_prime = true;
        for (int p : primes) {
            if (1LL * p * p > candidate) {
                break;
            }
            if (candidate % p == 0) {
                is_prime = false;
                break;
            }
        }
        if (is_prime) {
            primes.push_back(candidate);
            break;
        }
        ++candidate;
    }

    i128 total = 0;
    for (std::size_t i = 0; i + 1 < primes.size(); ++i) {
        const u64 p = static_cast<u64>(primes[i]);
        const u64 q = static_cast<u64>(primes[i + 1]);

        const u64 p2 = p * p;
        if (p2 > limit) {
            break;
        }

        const u64 q2 = q * q;
        const u64 left = p2 + 1;
        const u64 right = std::min(limit, q2 - 1);
        if (left > right) {
            continue;
        }

        const i128 sum_p = sum_multiples_in_range(p, left, right);
        const i128 sum_q = sum_multiples_in_range(q, left, right);
        const i128 sum_pq = sum_multiples_in_range(p * q, left, right);
        total += sum_p + sum_q - 2 * sum_pq;
    }

    return total;
}

u64 solve_u64(const u64 limit) {
    return static_cast<u64>(solve_i128(limit));
}

bool run_checkpoints() {
    if (solve_u64(15ULL) != 30ULL) {
        std::cerr << "Checkpoint failed for limit 15" << '\n';
        return false;
    }
    if (solve_u64(1000ULL) != 34825ULL) {
        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 << i128_to_string(solve_i128(options.limit)) << '\n';
    return 0;
}

Python

import math

def primes_up_to(n):
    is_composite = bytearray(n + 1)
    primes = []
    
    for i in range(2, n + 1):
        if not is_composite[i]:
            primes.append(i)
        for p in primes:
            v = i * p
            if v > n:
                break
            is_composite[v] = 1
            if i % p == 0:
                break
    return primes

def sum_multiples_in_range(d, left, right):
    if left > right:
        return 0
        
    first = ((left + d - 1) // d) * d
    last = (right // d) * d
    
    if first > last:
        return 0
        
    count = (last - first) // d + 1
    return (first + last) * count // 2

def solve(limit=999966663333):
    sqrt_limit = int(math.sqrt(limit)) + 10
    primes = primes_up_to(sqrt_limit + 100)
    
    # Ensure one prime above sqrt(limit) exists
    candidate = sqrt_limit + 101
    while True:
        is_prime = True
        for p in primes:
            if p * p > candidate:
                break
            if candidate % p == 0:
                is_prime = False
                break
        if is_prime:
            primes.append(candidate)
            break
        candidate += 1
        
    total = 0
    for i in range(len(primes) - 1):
        p = primes[i]
        q = primes[i + 1]
        
        p2 = p * p
        if p2 > limit:
            break
            
        q2 = q * q
        left = p2 + 1
        right = min(limit, q2 - 1)
        
        if left > right:
            continue
            
        sum_p = sum_multiples_in_range(p, left, right)
        sum_q = sum_multiples_in_range(q, left, right)
        sum_pq = sum_multiples_in_range(p * q, left, right)
        
        total += sum_p + sum_q - 2 * sum_pq
        
    return str(total)

if __name__ == '__main__':
    print(solve())

Java

import java.util.ArrayList;
import java.util.List;

public class Euler234 {
    static List<Integer> primesUpTo(int n) {
        byte[] isComposite = new byte[n + 1];
        List<Integer> primes = new ArrayList<>();

        for (int i = 2; i <= n; ++i) {
            if (isComposite[i] == 0) {
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) i * p;
                if (v > n)
                    break;
                isComposite[(int) v] = 1;
                if (i % p == 0)
                    break;
            }
        }
        return primes;
    }

    static long sumMultiplesInRange(long d, long left, long right) {
        if (left > right)
            return 0;

        long first = ((left + d - 1) / d) * d;
        long last = (right / d) * d;

        if (first > last)
            return 0;

        long count = (last - first) / d + 1;

        // Use BigInteger arithmetic to avoid overflow or cast carefully to avoid sign
        // extension issues?
        // Since limits are up to 10^12, first+last is 2*10^12, and count can be large.
        // Product is < 10^18, so long is safe for intermediate values in this problem.
        long sum = first + last;

        // Account for odd parity to divide by 2 safely
        if (sum % 2 == 0) {
            return (sum / 2) * count;
        } else {
            return sum * (count / 2);
        }
    }

    public static String solve() {
        long limit = 999966663333L;
        int sqrtLimit = (int) Math.sqrt(limit) + 10;
        List<Integer> primes = primesUpTo(sqrtLimit + 100);

        int candidate = sqrtLimit + 101;
        while (true) {
            boolean isPrime = true;
            for (int p : primes) {
                if ((long) p * p > candidate)
                    break;
                if (candidate % p == 0) {
                    isPrime = false;
                    break;
                }
            }
            if (isPrime) {
                primes.add(candidate);
                break;
            }
            candidate++;
        }

        long total = 0;
        for (int i = 0; i < primes.size() - 1; ++i) {
            long p = primes.get(i);
            long q = primes.get(i + 1);

            long p2 = p * p;
            if (p2 > limit)
                break;

            long q2 = q * q;
            long left = p2 + 1;
            long right = Math.min(limit, q2 - 1);

            if (left > right)
                continue;

            long sumP = sumMultiplesInRange(p, left, right);
            long sumQ = sumMultiplesInRange(q, left, right);
            long sumPQ = sumMultiplesInRange(p * q, left, right);

            total += sumP + sumQ - 2 * sumPQ;
        }

        return String.valueOf(total);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}