Problem 200: Find the 200th Prime-proof Sqube Containing the Contiguous Sub-string "200"

Project Euler Problem 200 Solution

EulerSolve provides an optimized solution for Project Euler Problem 200, Find the 200th Prime-proof Sqube Containing the Contiguous Sub-string "200", with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A sqube is a number of the form $$n=p^2q^3,\qquad p\neq q,\qquad p,q\text{ prime}.$$ The task is to find the 200th sqube whose decimal expansion contains the contiguous substring "200" and which is prime-proof , meaning that every valid one-digit change produces a composite number. The arithmetic form $p^2q^3$ is the main structural restriction, and the decimal pattern and prime-proof property are two additional filters applied to that sparse family. There is no closed formula for the 200th answer. The practical route is therefore to enumerate all squbes below a bound, sort them, and test each candidate with exact criteria that mirror the statement of the problem. Mathematical Approach The solution is best understood as a search over a highly structured set, not as a scan over all integers. The number-theoretic part identifies which values can even be candidates, and the decimal and primality parts then certify which of those candidates survive....

Detailed mathematical approach

Problem Summary

A sqube is a number of the form $p^2q^3$, where $p$ and $q$ are distinct primes, such as $200=5^2\cdot2^3$ or $120072949=23^2\cdot61^3$; the first five squbes are $72, 108, 200, 392$ and $500$. A number is prime-proof if no change of a single decimal digit turns it into a prime, and $200$ is the smallest such number. The prime-proof squbes whose decimal representation contains the three consecutive digits $200$ begin with $200$ and $1992008$, and we must find the $200$th of them. ^[1] Changing a digit means replacing it by a different digit; like the published programs, we never replace the first digit by $0$, and Section 4 shows that the other reading gives the same answer.

The whole solution in one sentence: we list all squbes up to a bound in increasing order, keep those whose digits contain $200$ and for which every one-digit change is composite, decided by a fast primality test, and enlarge the bound until $200$ squbes are kept.

Mathematical Approach

Squbes are rare, so we build them from pairs of primes instead of scanning all integers, and test the two digit conditions only on them.

1. The smallest squbes

The seven squbes below $1000$ are $72=3^2\cdot2^3$, $108=2^2\cdot3^3$, $200=5^2\cdot2^3$, $392=7^2\cdot2^3$, $500=2^2\cdot5^3$, $675=5^2\cdot3^3$ and $968=11^2\cdot2^3$. The first two use the same primes in exchanged roles, so it matters which prime is squared and which is cubed. The squbes form sequence A143610 of the On-Line Encyclopedia of Integer Sequences. ^[2]

Claim. Every sqube $x$ comes from exactly one ordered pair $(p,q)$ of distinct primes with $x=p^2q^3$.

Proof. By the fundamental theorem of arithmetic, the factorisation of $x$ into primes is unique up to the order of the factors, ^[3] and $x=p^2q^3$ is that factorisation. So $p$ is the only prime that occurs in it exactly twice, and $q$ the only one that occurs exactly three times. $\square$

Running through all ordered pairs of distinct primes therefore produces every sqube exactly once. The numbers $p^3q^2$ are the same squbes with the letters exchanged, so a second loop over them only repeats the list.

2. Listing all squbes up to a bound

Fix a bound $L$. The smallest cube of a prime is $2^3=8$ and the smallest square is $2^2=4$, so every sqube $p^2q^3\le L$ satisfies

$$p\le\sqrt{L/8},\qquad q\le\sqrt[3]{L/4}.$$

For $L\ge32$ the second bound is at most the first, because $(L/4)^2\le(L/8)^3$ is equivalent to $L\ge32$, so one list of the primes up to $\sqrt{L/8}$ serves for $p$ and for $q$. The sieve of Eratosthenes produces it: going through $2, 3, 4,\dots$, it takes each number $p$ that is still unmarked, which is a prime, and marks $p^2, p^2+p,\dots$ as composite. It may start at $p^2$, since a smaller multiple $kp$ with $2\le k\lt p$ was already marked by a prime factor of $k$, and stop once $p^2$ passes the end of the range, since a composite has a prime factor at most its square root (Section 8). ^[4] For fixed $p$ the products grow with $q$, so the loop over $q$ stops at the first product above $L$, and the loop over $p$ stops once $8p^2\gt L$. The products come out of order, since $p=2$ gives $108, 500, 1372,\dots$ and $p=3$ gives $72, 1125, 3087,\dots$, so the list is sorted afterwards.

There are $7$ squbes below $10^3$, $197$ below $10^6$, $3959$ below $10^9$ and $85026$ below $10^{12}$. For a fixed cube prime $q$, the prime $p$ runs up to $\sqrt{L/q^3}$, so by the prime number theorem, which gives about $y/\ln y$ primes up to $y$, ^[5] this family has roughly $\sqrt{L/q^3}/\ln\sqrt{L/q^3}$ members. Small cube primes dominate: below $10^{12}$ the $30253$ squbes $8p^2$ are more than a third of all. Summed over $q$, the count grows roughly like $\sqrt L/\ln L$, so for large $L$ a doubling multiplies it by about $\sqrt2\,\ln L/\ln(2L)$, a little less than $\sqrt2$.

3. Finding the digits 200

Write $x=\sum_i a_i10^i$ with decimal digits $a_i$, where $a_0$ is the units digit. Dividing by $10^k$ and rounding down removes the last $k$ digits, and the remainder modulo $1000$ keeps the next three:

$$\begin{aligned}&\Bigl\lfloor\frac{x}{10^k}\Bigr\rfloor\bmod1000\\&\quad=100a_{k+2}+10a_{k+1}+a_k.\end{aligned}$$

This equals $200$ exactly when $a_{k+2}=2$ and $a_{k+1}=a_k=0$. So $x$ contains the block $200$ if and only if this remainder is $200$ for some $k\ge0$, and only the $k$ with $\lfloor x/10^k\rfloor\ge200$ need to be tried. For $1992008$ the remainders for $k=0,1,2,3,4$ are $8, 200, 920, 992, 199$, so the block is found at $k=1$; the example $120072949$ contains it at $k=5$.

4. Prime-proof numbers

A number with $d$ digits has $9d-1$ one-digit changes: each of the $d-1$ lower digits can be replaced by $9$ other digits, and the leading digit by $8$, since it may not become $0$. For $200$ the $26$ changes are

$100, 300, 400,\dots,900$, all multiples of $100$;
$210, 220,\dots,290$, all even;
$201,\dots,209$: the even ones, $205=5\cdot41$, $201=3\cdot67$, $203=7\cdot29$, $207=3^2\cdot23$ and $209=11\cdot19$.

None of them is prime, so $200$ is prime-proof. Checking all smaller numbers shows that each has a prime among its changes, for example $100$ becomes the prime $101$, so $200$ indeed comes first. The prime-proof numbers form sequence A118118, ^[6] but there a leading digit may become $0$, so $202$ and $205$, whose changes $002=2$ and $005=5$ are primes, are missing. For every sqube below $10^{12}$ that contains $200$, however, both readings agree.

5. The last digit decides

Claim. Let $x\ge10$ end in $0, 2, 4, 5, 6$ or $8$, and put $b=10\lfloor x/10\rfloor$. Then $x$ is prime-proof if and only if none of the four numbers $b+1$, $b+3$, $b+7$, $b+9$ is prime.

Proof. A change of any digit except the last keeps the last digit, so the new number is divisible by $2$ or by $5$. It has as many digits as $x$, because the leading digit never becomes $0$, so it is at least $10$ and therefore composite. A change of the last digit gives $b+c$ with a digit $c$ different from the last digit of $x$. For $c\in\{0,2,4,5,6,8\}$ this number is again divisible by $2$ or $5$ and at least $10$, hence composite. The digits $c=1,3,7,9$ all differ from the last digit of $x$ and give the four numbers of the claim, which are therefore the only changes that can be prime. $\square$

So such an $x$ is prime-proof exactly when its decade $b,\dots,b+9$ contains no prime. The decade from $200$ to $209$ is the first one without a prime, and the first prime-proof numbers $200, 202, 204, 205, 206, 208$ all lie in it. For $1992008$ the four numbers are $1992001=11\cdot47\cdot3853$, $1992003=3\cdot13^2\cdot3929$, $1992007=23\cdot257\cdot337$ and $1992009=3\cdot17\cdot139\cdot281$, so $1992008$ is prime-proof. A sqube is divisible by $2$ or $5$ exactly when $2$ or $5$ is one of $p$ and $q$, and for these squbes four primality tests settle the question.

6. Squbes ending in 1, 3, 7 or 9

If $x$ ends in $1, 3, 7$ or $9$, changes of the other digits keep that last digit, and Section 5 no longer rules them out: each of the $9d-10$ changes of the other digits and of three changes of the last digit may be prime, $101$ numbers when $d=12$. With so many candidates a prime almost always appears. The smallest prime-proof number of this kind is $212159=53\cdot4003$, ^[6] and of the $44819$ squbes below $10^{12}$ that end in $1, 3, 7$ or $9$ only two are prime-proof, $21\,683\,329\,727=457^2\cdot47^3$ and $283\,998\,957\,017=7603^2\cdot17^3$, neither containing $200$. By contrast, $30022$ of the $40207$ squbes below $10^{12}$ that are divisible by $2$ or $5$ are prime-proof. The example $120072949$ is typical: it contains $200$, but its change $120072947$ is prime, and so are six other changes.

Of the $933$ squbes below $10^{12}$ that contain $200$, $478$ are divisible by $2$ or $5$, and $332$ of these are prime-proof, while none of the other $455$ is. So more than $200$ qualifying squbes lie below $10^{12}$.

7. Why the bound can grow

Claim. If the sorted list of qualifying squbes up to $L$ has at least $k$ entries, its $k$-th entry is the $k$-th qualifying sqube overall.

Proof. Whether a sqube qualifies depends only on the number, not on $L$. By Section 2 the list contains every qualifying sqube up to $L$, and every other one is larger than $L$, hence larger than all entries. $\square$

A program may therefore start with any bound and double it until $200$ squbes qualify, or take a bound known to suffice, such as $10^{12}$ by Section 6.

8. Trial division

If $n=ab$ with $1\lt a\le b$, then $a^2\le n$, so a composite $n$ has a prime factor at most $\sqrt n$. Every prime above $3$ leaves the remainder $1$ or $5$ on division by $6$, because $6k$, $6k+2$, $6k+3$ and $6k+4$ are divisible by $2$ or $3$. ^[7] After $2$ and $3$ it is therefore enough to try the pairs $i$, $i+2$ for $i=5, 11, 17,\dots$ while $i^2\le n$: about $\sqrt n/6\approx1.7\cdot10^5$ pairs for a prime $n$ near $10^{12}$, and far fewer for most composite $n$, which have a small prime factor.

9. The Miller–Rabin test

For an odd $n\gt2$, write $n-1=2^sd$ with $d$ odd. If $n$ is prime and does not divide $a$, Fermat's little theorem gives $a^{n-1}\equiv1\pmod n$, and in the chain $a^d, a^{2d},\dots,a^{2^sd}=a^{n-1}$ each number is the square of the one before. Modulo a prime the only square roots of $1$ are $\pm1$, because $n\mid(y-1)(y+1)$ forces $n\mid y-1$ or $n\mid y+1$; so the first $1$ in the chain is either $a^d$ itself or comes right after a $-1$. Hence $a^d\equiv1$, or $a^{2^rd}\equiv-1$ for some $r$ with $0\le r\lt s$; if neither holds, $a$ proves that $n$ is composite. ^[8]

For $n=2047=23\cdot89$ we have $n-1=2\cdot1023$. Since $2^{11}=2048\equiv1$, also $2^{1023}=(2^{11})^{93}\equiv1$, so the base $2$ is fooled; but $3^{1023}\equiv1565$ is neither $1$ nor $-1\equiv2046$, so the base $3$ exposes $n$. A few fixed bases are enough: the seven bases $2$, $325$, $9375$, $28178$, $450775$, $9780504$, $1795265022$ found by Jim Sinclair decide every $n\lt2^{64}$ correctly, ^[9] and so do the twelve primes up to $37$. ^[8] Each base costs one modular power, which binary exponentiation computes by reading the exponent bit by bit: it squares the base at every step and multiplies it into the result where the bit is $1$, about $\log_2n$ squarings and at most as many multiplications. ^[10]

How the Code Works

In the C++ program, ^[11] parse_arguments reads the optional --target=, default $200$, --initial-limit=, default $10^9$, and --skip-checkpoints; an unknown option, an empty or non-digit value, a target below $1$ or a bound below $100$ gives exit code $1$, but the digit parsers do not check for overflow, so an oversized number is not reliably rejected. Unless the checkpoints are skipped, run_checkpoints requires first_squbes(5) to return $72, 108, 200, 392, 500$ and solve(2, 2000000) to return $1992008$, the facts of the statement, or exits with code $2$. Then main prints the result of solve.

solve runs Section 7 in a loop: generate_squbes(limit) lists the squbes, those that pass contains_200 and is_prime_proof are kept, and if fewer than the target are found, limit is doubled and everything is computed again. generate_squbes follows Section 2: primes_up_to sieves up to $\lfloor\sqrt{L/8}\rfloor+5$ in a std::vector<std::uint8_t>, the products are formed as unsigned 128-bit integers and stored as std::uint64_t, and the list is sorted; the std::unique after sorting is only a safeguard, as Section 1 rules out repeats. first_squbes doubles a bound from $1000$ in the same way, and contains_200 is the loop of Section 3.

is_prime_proof changes digits arithmetically: for pow10 $=1, 10, 100,\dots$ it subtracts the digit at that place times pow10, adds every other digit times pow10, skipping $0$ at the leading place, and returns false at the first prime. Because it starts with the units digit, for a sqube divisible by $2$ or $5$ it tries first the only changes that can be prime by Section 5. is_prime divides by the twelve primes up to $37$ and then runs Section 9 with Sinclair's bases, skipping a base that is a multiple of $n$. ^[9] pow_mod is the binary exponentiation of Section 9, and mul_mod multiplies two residues as an unsigned __int128 before reducing, ^[12] so nothing overflows for $n\lt2^{64}$.

The Python and Java programs fix the bound at $10^{12}$, enough by Section 6, and have no options or checkpoints. They put $p^2q^3$ and also $p^3q^2$ into a set, a Python set sorted afterwards or a Java TreeSet<Long>; by Section 1 the second loop only repeats squbes, which the set ignores. Both change digits in the decimal string from the leading digit on, convert back with int or Long.parseLong, and test primality by the trial division of Section 8. Python integers have arbitrary precision, and in Java every product the loops form stays below $5\cdot10^{12}$, far inside the range of the 64-bit long.

Complexity and Verification

Count each arithmetic operation on 64-bit numbers as one step, and let $M$ squbes lie below $L$, $S$ of them containing $200$. The sieve costs $O(\sqrt L\log\log L)$, ^[4] the generation loops $O(M+\sqrt L)$, and the ordering $O(M\log M)$, by a final sort in C++ and Python and by the insertions into the TreeSet in Java. Each of the $S$ candidates needs at most $9d-1$ primality tests, $d$ being its number of digits, and stops at the first prime. A C++ test costs at most twelve divisions and seven modular powers of about $40$ squarings, so all tests cost $O(Sd\log L)$; trial division needs up to $\sqrt n/6$ loop steps for a prime $n$, which bounds Python and Java by $O(Sd\sqrt L)$, far above the actual work. The memory is $O(\sqrt L+M)$. In the default run each doubling multiplied the number of squbes by about $1.36$, so all C++ rounds together generated about $3.5$ times as many as the last round.

Instrumented copies of the programs count the work: in the main search, without the checkpoints, the C++ program calls is_prime $72756$ times over all its rounds, and only $3924$ calls get past the trial division to Miller–Rabin; Python and Java call their primality test $29634$ times before the $200$th match, with $21\,076\,981$ trial-division steps in total. On a laptop with an Apple M4 Max chip, using GCC 16 with -O2, OpenJDK 17 and CPython 3.11, a run takes about $0.01$ s in C++, $0.1$ s in Java and $1.2$ s in Python, with peak memory of about $10$, $70$ and $28$ MB.

For this article all three programs were run: the C++ checkpoints pass, and all three print the same number. The C++ is_prime agrees with a sieve for every $n\le10^8$ and with SymPy's isprime on all $95790$ one-digit changes of the $933$ squbes below $10^{12}$ that contain $200$; contains_200 agrees with a string search for every $x\le10^7$; the criterion of Section 5 agrees with the full test on all $478$ of these squbes that are divisible by $2$ or $5$. An independent Python program using SymPy, once with the full test and once with that criterion, finds the same $200$th sqube. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the criterion of Section 5 and all counts are derived or computed in this article.

Project Euler 200 — Prime-proof Squbes. The official statement defines squbes and prime-proof numbers, names 200 and 1992008 as the first two prime-proof squbes containing 200, and asks for the 200th prime-proof sqube containing the contiguous sub-string 200.
OEIS A143610. The numbers p²·q³ with distinct primes p and q: 72, 108, 200, 392, 500, 675, 968, … The entry refers to this problem.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a prime or a product of primes, unique up to the order of the factors.
Wikipedia — Sieve of Eratosthenes. Marking the multiples of p may start at p² and stop once p² exceeds n; the sieve needs O(n log log n) operations.
Wikipedia — Prime number theorem. The number π(x) of primes up to x satisfies π(x) ∼ x / ln x.
OEIS A118118. Composite numbers that stay composite when any single decimal digit is changed, called prime-proof after Project Euler: 200, 204, 206, 208, 320, …, without 202 and 205. A comment gives 212159 as the first odd term not ending in 5.
Wikipedia — Primality test. Trial division needs only divisors up to √n, and since every prime above 3 is 1 or 5 modulo 6, it suffices to test 2, 3 and the numbers 6k + 1 and 6k + 5.
Wikipedia — Miller–Rabin primality test. The strong probable prime test with n − 1 = 2ˢd, proved from the fact that ±1 are the only square roots of 1 modulo a prime; the bases 2, 3, 5, …, 37 suffice for every n < 2⁶⁴.
Deterministic variants of the Miller–Rabin primality test. The page lists Jim Sinclair's seven bases 2, 325, 9375, 28178, 450775, 9780504, 1795265022 for every n < 2⁶⁴ and warns that a base which is a multiple of n needs special handling.
Wikipedia — Modular exponentiation. The right-to-left binary method squares the base for every bit of the exponent e and multiplies it into the result where the bit is 1, reducing modulo m, in O(log e) steps.
C++ — Euler200.cpp. The linked immutable revision contains is_prime, is_prime_proof, generate_squbes, solve and run_checkpoints. It is the implementation record for this article.
GCC manual — 128-bit Integers. GCC provides unsigned __int128 on targets with an integer mode wide enough for 128 bits.

Problem 200 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    int target = 200;
    u64 initial_limit = 1000000000ULL;
    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_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_int_after_prefix(arg, "--target=", options.target) ||
            parse_u64_after_prefix(arg, "--initial-limit=", options.initial_limit)) {
            continue;
        }

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

    return options.target >= 1 && options.initial_limit >= 100;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    u64 result = 1 % mod;
    base %= mod;

    while (exp > 0) {
        if (exp & 1ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1;
    }
    return result;
}

bool is_prime(const u64 n) {
    if (n < 2) {
        return false;
    }
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0) {
            return false;
        }
    }

    u64 d = n - 1;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1;
        ++s;
    }

    constexpr u64 bases[] = {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};
    for (u64 a : bases) {
        if (a % n == 0) {
            continue;
        }
        u64 x = pow_mod(a, d, n);
        if (x == 1 || x == n - 1) {
            continue;
        }

        bool witnessed = true;
        for (int r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1) {
                witnessed = false;
                break;
            }
        }
        if (witnessed) {
            return false;
        }
    }

    return true;
}

std::vector<int> primes_up_to(const int n) {
    std::vector<std::uint8_t> is_prime_small(static_cast<std::size_t>(n + 1), 1);
    is_prime_small[0] = 0;
    is_prime_small[1] = 0;

    for (int i = 2; i * i <= n; ++i) {
        if (!is_prime_small[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            is_prime_small[static_cast<std::size_t>(j)] = 0;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime_small[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

bool contains_200(u64 x) {
    while (x >= 200) {
        if (x % 1000ULL == 200ULL) {
            return true;
        }
        x /= 10ULL;
    }
    return false;
}

bool is_prime_proof(const u64 x) {
    u64 pow10 = 1;
    while (pow10 <= x) {
        const int digit = static_cast<int>((x / pow10) % 10ULL);
        const u64 base = x - static_cast<u64>(digit) * pow10;

        for (int repl = 0; repl <= 9; ++repl) {
            if (repl == digit) {
                continue;
            }
            if (pow10 * 10ULL > x && repl == 0) {
                continue;
            }

            const u64 candidate = base + static_cast<u64>(repl) * pow10;
            if (is_prime(candidate)) {
                return false;
            }
        }

        pow10 *= 10ULL;
    }

    return true;
}

std::vector<u64> generate_squbes(const u64 limit) {
    const u64 max_p = static_cast<u64>(std::sqrt(static_cast<long double>(limit / 8ULL))) + 5ULL;
    const std::vector<int> primes = primes_up_to(static_cast<int>(max_p));

    std::vector<u64> squbes;
    for (int p : primes) {
        const u128 p2 = static_cast<u128>(p) * static_cast<u128>(p);
        if (p2 * 8U > limit) {
            break;
        }

        for (int q : primes) {
            if (q == p) {
                continue;
            }

            const u128 q3 = static_cast<u128>(q) * q * q;
            const u128 value = p2 * q3;
            if (value > limit) {
                break;
            }

            squbes.push_back(static_cast<u64>(value));
        }
    }

    std::sort(squbes.begin(), squbes.end());
    squbes.erase(std::unique(squbes.begin(), squbes.end()), squbes.end());
    return squbes;
}

std::vector<u64> first_squbes(const std::size_t k) {
    u64 limit = 1000;
    while (true) {
        std::vector<u64> squbes = generate_squbes(limit);
        if (squbes.size() >= k) {
            squbes.resize(k);
            return squbes;
        }
        limit *= 2;
    }
}

u64 solve(const int target, u64 limit) {
    while (true) {
        const std::vector<u64> squbes = generate_squbes(limit);

        std::vector<u64> matches;
        matches.reserve(squbes.size() / 4);

        for (u64 x : squbes) {
            if (!contains_200(x)) {
                continue;
            }
            if (is_prime_proof(x)) {
                matches.push_back(x);
            }
        }

        if (static_cast<int>(matches.size()) >= target) {
            return matches[static_cast<std::size_t>(target - 1)];
        }

        limit *= 2ULL;
    }
}

bool run_checkpoints() {
    const std::vector<u64> first = first_squbes(5);
    const std::vector<u64> expected = {72ULL, 108ULL, 200ULL, 392ULL, 500ULL};
    if (first != expected) {
        std::cerr << "Checkpoint failed for first squbes" << '\n';
        return false;
    }

    if (solve(2, 2000000ULL) != 1992008ULL) {
        std::cerr << "Checkpoint failed for second prime-proof sqube containing 200" << '\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.target, options.initial_limit) << '\n';
    return 0;
}

Python

# Problem 200: Find the 200th prime-proof sqube containing "200".
# A sqube is p^2*q^3 or p^3*q^2 for distinct primes p,q.
# Prime-proof: replacing any digit doesn't produce a prime.

from math import isqrt

def solve():
    def sieve_primes(n):
        s = bytearray(b'\x01') * (n+1)
        s[0] = s[1] = 0
        for i in range(2, isqrt(n)+1):
            if s[i]:
                for j in range(i*i, n+1, i): s[j] = 0
        return [i for i in range(n+1) if s[i]]

    def is_prime(n):
        if n < 2: return False
        if n < 4: return True
        if n % 2 == 0 or n % 3 == 0: return False
        i = 5
        while i*i <= n:
            if n % i == 0 or n % (i+2) == 0: return False
            i += 6
        return True

    def contains_200(x):
        while x >= 200:
            if x % 1000 == 200: return True
            x //= 10
        return False

    def is_prime_proof(x):
        s = str(x)
        for i in range(len(s)):
            orig = int(s[i])
            for d in range(10):
                if d == orig: continue
                if i == 0 and d == 0: continue
                candidate = int(s[:i] + str(d) + s[i+1:])
                if is_prime(candidate): return False
        return True

    limit = 10**12
    primes = sieve_primes(isqrt(limit) + 10)
    squbes = set()
    for p in primes:
        p2 = p*p
        if p2 * 8 > limit: break
        for q in primes:
            if q == p: continue
            val = p2 * q * q * q
            if val > limit: break
            squbes.add(val)
    for p in primes:
        p3 = p*p*p
        if p3 * 4 > limit: break
        for q in primes:
            if q == p: continue
            val = p3 * q * q
            if val > limit: break
            squbes.add(val)

    squbes = sorted(squbes)
    count = 0
    for x in squbes:
        if not contains_200(x): continue
        if is_prime_proof(x):
            count += 1
            if count == 200:
                print(x)
                return

solve()

Java

import java.util.*;

public class Euler200 {
    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        if (n < 4)
            return true;
        if (n % 2 == 0 || n % 3 == 0)
            return false;
        for (long i = 5; i * i <= n; i += 6)
            if (n % i == 0 || n % (i + 2) == 0)
                return false;
        return true;
    }

    static boolean contains200(long x) {
        while (x >= 200) {
            if (x % 1000 == 200)
                return true;
            x /= 10;
        }
        return false;
    }

    static boolean isPrimeProof(long x) {
        String s = Long.toString(x);
        for (int i = 0; i < s.length(); i++) {
            int orig = s.charAt(i) - '0';
            for (int d = 0; d <= 9; d++) {
                if (d == orig)
                    continue;
                if (i == 0 && d == 0)
                    continue;
                long c = Long.parseLong(s.substring(0, i) + d + s.substring(i + 1));
                if (isPrime(c))
                    return false;
            }
        }
        return true;
    }

    public static void main(String[] args) {
        long limit = (long) 1e12;
        int maxP = (int) Math.sqrt(limit / 8.0) + 10;
        boolean[] sieve = new boolean[maxP + 1];
        Arrays.fill(sieve, true);
        sieve[0] = sieve[1] = false;
        for (int i = 2; i * i <= maxP; i++)
            if (sieve[i])
                for (int j = i * i; j <= maxP; j += i)
                    sieve[j] = false;
        int[] primes = new int[maxP];
        int pc = 0;
        for (int i = 2; i <= maxP; i++)
            if (sieve[i])
                primes[pc++] = i;
        TreeSet<Long> squbes = new TreeSet<>();
        for (int i = 0; i < pc; i++) {
            long p = primes[i];
            if (p * p * 8 > limit)
                break;
            for (int j = 0; j < pc; j++) {
                if (i == j)
                    continue;
                long q = primes[j];
                long v = p * p * q * q * q;
                if (v > limit)
                    break;
                squbes.add(v);
            }
        }
        for (int i = 0; i < pc; i++) {
            long p = primes[i];
            if (p * p * p * 4 > limit)
                break;
            for (int j = 0; j < pc; j++) {
                if (i == j)
                    continue;
                long q = primes[j];
                long v = p * p * p * q * q;
                if (v > limit)
                    break;
                squbes.add(v);
            }
        }
        int count = 0;
        for (long x : squbes) {
            if (!contains200(x))
                continue;
            if (isPrimeProof(x)) {
                count++;
                if (count == 200) {
                    System.out.println(x);
                    return;
                }
            }
        }
    }
}