Problem 111: Primes with Runs

Project Euler Problem 111 Solution

EulerSolve provides an optimized solution for Project Euler Problem 111, Primes with Runs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each digit $d \in \{0,1,\dots,9\}$, let $M(10,d)$ be the largest number of times that $d$ can appear in a 10-digit prime, and let $S(10,d)$ be the sum of all 10-digit primes that attain this maximum repetition count. The goal is to compute $$\sum_{d=0}^{9} S(10,d).$$ The important point is that the optimization is digit-specific. For each fixed digit $d$, we are not looking for all primes with many repeated digits in general; we are looking for primes with as many copies of that one chosen digit as possible. Mathematical Approach The implementations solve the problem by turning "maximum repetition" into the equivalent question "how few positions must be changed before a prime appears?" That reformulation matches the combinatorial structure of the search exactly. Minimal replacement count instead of maximal repetition Fix a digit $d$. Let $r_d$ be the smallest integer $r \ge 1$ for which there exists a 10-digit prime with exactly $r$ digits different from $d$....

Detailed mathematical approach

Problem Summary

A prime with $n$ digits is a prime from $10^{n-1}$ to $10^n-1$, so its first digit is not $0$. For a digit $d$ from $0$ to $9$, $M(n,d)$ is the largest number of times that $d$ occurs among the digits of an $n$-digit prime, $N(n,d)$ is the number of $n$-digit primes in which $d$ occurs exactly $M(n,d)$ times, and $S(n,d)$ is the sum of these primes. The copies of $d$ need not stand next to each other, despite the word "runs" in the title: $1151$ counts as a prime with three ones. The statement works out four digits completely. No four-digit prime has four equal digits, since $1111=11\cdot101$, $2222$ is even, and so on. Nine four-digit primes contain three ones, namely $1117$, $1151$, $1171$, $1181$, $1511$, $1811$, $2111$, $4111$ and $8111$, so $M(4,1)=3$, $N(4,1)=9$ and $S(4,1)=22275$; for the digit $0$ only $M(4,0)=2$ is possible, with $N(4,0)=13$ primes. For $d=0,1,\dots,9$ the whole table reads

$M(4,d)$: $2, 3, 3, 3, 3, 3, 3, 3, 3, 3$
$N(4,d)$: $13, 9, 1, 12, 2, 1, 1, 9, 1, 7$
$S(4,d)$: $67061, 22275, 2221, 46214, 8888,$
$5557, 6661, 57863, 8887, 48073$

and the ten values $S(4,d)$ add up to $273700$. We must find the sum of $S(10,d)$ over all ten digits $d$. ^[1]

The whole solution in one sentence: for each digit $d$ we start from the string $dd\dots d$ of ten copies of $d$, change as few of its places as possible, first one and then two, test every ten-digit number made in this way for primality, and stop at the first number of changed places that produces a prime, because the primes found there are exactly the primes with the most copies of $d$. The sections below explain why this is right, why one or two changes are always enough for ten digits, and how a ten-digit number is proved prime with a few hundred multiplications.

Mathematical Approach

We first redo a few rows of the four-digit table by hand, then turn the question from "most copies" into "fewest changes", count the candidates, remove most of them with simple divisibility rules, estimate how many primes to expect, run the search for ten digits, explain the primality test that the programs use, and finally look at other lengths.

1. Four digits by hand

Take $d=2$. The number $2222$ is even, so no four-digit prime has four twos. A number with three twos has exactly one other digit, and that digit must stand in the last place: otherwise the number ends in $2$ and is even. The last digit cannot be $0$, $4$, $6$ or $8$ either, and not $5$, because then the number is divisible by $5$; so it is $1$, $3$, $7$ or $9$, and only four numbers remain: $2221$, $2223$, $2227$ and $2229$. Now $2223=3\cdot741$, $2229=3\cdot743$ and $2227=17\cdot131$, while $2221$ is prime: none of the fifteen primes up to $47$ divides it, and Section 10 shows why primes up to $\sqrt{2221}\approx47.1$ are enough. So $M(4,2)=3$, $N(4,2)=1$ and $S(4,2)=2221$, exactly as in the table.

Take $d=0$. A four-digit number cannot start with $0$, and a prime with more than one digit cannot end in $0$, because then it is divisible by $10$. So only the two middle digits can be zeros, and the candidates are the $36$ numbers $a00b$ with a first digit $a$ from $1$ to $9$ and a last digit $b$ from $1,3,7,9$. One of them, $1009$, is prime, so $M(4,0)=2$; Section 7 finds all $13$.

Take $d=1$. Four ones give $1111=11\cdot101$, which is not prime. Changing one digit, for example the last one into $7$, gives the prime $1117$, so $M(4,1)=3$. To find $N(4,1)$ and $S(4,1)$ we must look at every number with exactly three ones; Section 7 does this and finds the nine primes of the statement.

So the method is: change as few places as possible, and let the first and the last digit discard most candidates before any primality test.

2. From the most copies to the fewest changes

Fix $n$ and $d$. For a string $x$ of $n$ digits, in which for the moment a leading $0$ is allowed, count the places where the digit of $x$ is not $d$, and call this number $k$; then $d$ occurs exactly $n-k$ times in $x$. We say that $x$ lies on level $k$. Every such string, and in particular every $n$-digit number, lies on exactly one of the levels $0,1,\dots,n$. For $d=1$ and $n=4$, the primes $1117$ and $1151$ lie on level $1$, the prime $1009$ on level $3$ and the prime $2003$ on level $4$.

Levels turn the maximum into a minimum. Let $k^*$ be the smallest level that contains at least one prime. Then

$$M(n,d)=n-k^*,$$

and the primes in which $d$ occurs $M(n,d)$ times are exactly the primes on level $k^*$; so $N(n,d)$ is their number and $S(n,d)$ is their sum. The proof is short. A prime on level $k$ contains $n-k$ copies of $d$. The levels below $k^*$ contain no prime, so no prime has more than $n-k^*$ copies, and the primes on level $k^*$ have exactly that many. Conversely, a prime with $n-k^*$ copies of $d$ has $k^*$ other digits and therefore lies on level $k^*$. Since every $n$-digit prime lies on some level, $k^*$ exists as soon as there is an $n$-digit prime at all; for ten digits there are more than $400$ million (Section 8).

This is the algorithm: examine level $0$, then level $1$, then level $2$, and so on, and stop at the first level that contains a prime. The levels above $k^*$ are never looked at, and that is what makes the search fast; as Section 9 shows, for ten digits every digit stops at level $1$ or level $2$.

The requested total adds the ten sums $S(n,d)$, so a prime that is maximal for several digits counts several times. For two digits this is common, since $23$ is maximal for $0$, $2$ and $3$ (Section 11), but for four and ten digits it never happens, because the other digits of a maximal prime are too few to reach the maximum of a second digit.

3. Level 0: one repeated digit

Level $0$ holds a single string, the repdigit $dd\dots d$. With the repunit $R_n=11\dots1=(10^n-1)/9$, made of $n$ ones, it equals $d\cdot R_n$. For $d=0$ it is not an $n$-digit number at all. For $d\ge2$ and $n\ge2$ it is divisible by $d$ and larger than $d$, so it is not prime; the statement makes this point with $2222$, which is even divisible by $22$. Only $d=1$ needs a thought: is $R_n$ itself prime?

If $n=ab$ with $a,b\ge2$, then $R_n$ is divisible by $R_a$: cut the $ab$ ones into $b$ blocks of $a$ ones, so that

$$\begin{aligned}R_{ab}&=R_a\,\bigl(1+10^{a}+10^{2a}\\&\quad{}+\dots+10^{(b-1)a}\bigr).\end{aligned}$$

For example $R_4=1111=11\cdot101$ and $R_6=111111=111\cdot1001$. So $R_n$ can be prime only when $n$ is prime. For $n=10=2\cdot5$, both $R_2=11$ and $R_5=11111$ divide $R_{10}$: $R_{10}=11\cdot101010101$${}=11111\cdot100001$, and in full $R_{10}=11\cdot41\cdot271\cdot9091$. ^[2] Level $0$ therefore contains no prime for ten digits and any $d$, and the search may start at level $1$. A prime $n$ is not enough, though: $R_3=3\cdot37$ and $R_5=41\cdot271$. The repunit $R_n$ is prime only for $n=2,19,23,317,1031,\dots$, ^[2] and Section 11 shows what this means for other lengths.

4. How many numbers lie on a level

A number on level $k$ is built in two steps: choose the $k$ places that get a digit other than $d$, then choose that digit, one of nine, for each chosen place. The number of ways to choose $k$ places out of $n$ is the binomial coefficient

$$\binom nk=\frac{n!}{k!\,(n-k)!},$$

where $n!=1\cdot2\cdots n$ and $0!=1$. ^[3] For instance, two places out of four can be chosen in $\binom42=6$ ways: the first and second, first and third, first and fourth, second and third, second and fourth, or third and fourth. So level $k$ consists of

$$\binom nk\cdot9^k$$

strings of digits, some of which start with $0$ and are not $n$-digit numbers. For $n=4$ this gives $4\cdot9=36$ strings on level $1$ and $6\cdot81=486$ on level $2$. For $n=10$ it gives $90$ on level $1$, $45\cdot81=3645$ on level $2$, and already $120\cdot729=87480$ on level $3$.

Each string is produced exactly once, because both choices can be read off from the string itself: the chosen places are the places whose digit is not $d$, and the chosen digits are the digits standing there. So a search that runs through all choices meets every candidate exactly once, and no prime is ever counted twice.

5. The first and the last digit

Two facts about the decimal system remove most candidates without a single division. The first digit of an $n$-digit number is not $0$. The last digit decides divisibility by $2$ and by $5$: a number $x$ with last digit $b$ is $x=10q+b$, and since $10$ is divisible by $2$ and by $5$, $x$ is divisible by $2$ exactly when $b$ is, and by $5$ exactly when $b$ is. ^[4] A prime with at least two digits is neither $2$ nor $5$, so its last digit is $1$, $3$, $7$ or $9$. This has three consequences.

$d=0$: both the first and the last place must hold a digit other than $0$. So level $1$ contains no prime and $M(n,0)\le n-2$. On level $2$ the two changed places are forced, and the candidates are the $36$ numbers $a\cdot10^{n-1}+b$ with $a$ from $1$ to $9$ and $b$ from $1,3,7,9$.
$d=2,4,5,6,8$: the last place must be changed. On level $1$ it is the only changed place, and just four numbers remain: $dd\dots d1$, $dd\dots d3$, $dd\dots d7$ and $dd\dots d9$.
$d=1,3,7,9$: every place may be changed. For ten digits, level $1$ has $90$ strings; one of them starts with $0$, and six end in $0$, $2$, $4$, $5$, $6$ or $8$, which leaves $83$ candidates.

6. Digit sums and the factor 3

The next filter needs remainders. We write $a\equiv b\pmod m$, read "$a$ is congruent to $b$ modulo $m$", when $a-b$ is divisible by $m$, that is, when $a$ and $b$ leave the same remainder after division by $m$. For example $10\equiv1\pmod3$, because $10-1=9$. Congruences can be added and multiplied like equations, so $10^i\equiv1\pmod 3$ for every $i$, and a number $x$ with the digits $c_{n-1},\dots,c_1,c_0$ satisfies

$$x=\sum_{i=0}^{n-1}c_i\,10^i\equiv\sum_{i=0}^{n-1}c_i\pmod 3.$$

So every number leaves the same remainder modulo $3$ as its digit sum, and it is divisible by $3$ exactly when its digit sum is. ^[4] For example, $2223$ has the digit sum $9$, and $2223=3\cdot741$.

On level $1$ this filter is especially simple. For ten digits, the nine copies of $d$ contribute $9d$ to the digit sum, a multiple of $3$; for four digits, the three copies contribute $3d$. In both cases a candidate is divisible by $3$ exactly when its one changed digit is $0$, $3$, $6$ or $9$. For the four numbers $dd\dots db$ of Section 5 this removes $b=3$ and $b=9$ and leaves two candidates. Of the $83$ ten-digit candidates of $d=1$ it leaves $46$: each of the first nine places keeps the five replacements $2,4,5,7,8$, and the last place only $7$, so $9\cdot5+1=46$; for $d=7$ the same count, with $1,2,4,5,8$ and $1$, gives $46$ again. For $d=3$ and $d=9$ it leaves $9\cdot6+2=56$, because $d$ itself is one of the multiples of $3$: each of the first nine places keeps all six digits $1,2,4,5,7,8$, and the last place both $1$ and $7$.

A second remainder fact shows that level-$1$ candidates of ten digits avoid some small factors altogether. Such a candidate, with the digit $r$ in the place of value $10^i$, is

$$x=d\,R_{10}+(r-d)\,10^{i}.$$

If $p$ is one of the prime factors $11$, $41$, $271$, $9091$ of $R_{10}$, then $x\equiv(r-d)\,10^i\pmod p$, and this is not $0$: the difference $r-d$ lies between $-9$ and $9$ and is not $0$, so $p$ divides neither it nor $10^i$. So no such candidate is divisible by $11$, $41$, $271$ or $9091$. For $p=11$ this is the alternating-sum rule: ^[4] the alternating digit sum of a ten-digit repdigit is $0$, and changing one digit moves it by $\pm(r-d)$. The programs do not compute digit sums; their primality test rejects multiples of $3$ at once. The rules only explain the counts that the search meets.

7. The four-digit table, completely

Now every row of the table can be checked. For $d=1$, level $1$ has $36$ strings. The string $0111$ is not a four-digit number, and $1110$, $1112$, $1114$, $1115$, $1116$, $1118$ have a forbidden last digit, so $29$ candidates remain. Their digit sum is $3+r$ for the changed digit $r$, so those with $r=0,3,6,9$ are divisible by $3$, and $16$ candidates survive:

$2111, 4111, 5111, 7111, 8111,$
$1211, 1411, 1511, 1711, 1811,$
$1121, 1141, 1151, 1171, 1181,$
$1117.$

Seven of them are composite: $5111=19\cdot269$, $7111=13\cdot547$, $1211=7\cdot173$, $1411=17\cdot83$, $1711=29\cdot59$, $1121=19\cdot59$ and $1141=7\cdot163$. The other nine are the primes of the statement, and their sum is $22275$.

The digits $3$, $7$ and $9$ work the same way: the filters leave $20$, $16$ and $20$ candidates, of which $12$, $9$ and $7$ are prime. For $d=2,4,5,6,8$, the two candidates $ddd1$ and $ddd7$ remain. Here $2221$ is prime and $2227=17\cdot131$; $4441$ and $4447$ are both prime, which gives $N(4,4)=2$ and $S(4,4)=8888$; $5551=7\cdot13\cdot61$ and $5557$ is prime; $6661$ is prime and $6667=59\cdot113$; $8881=83\cdot107$ and $8887$ is prime. For $d=0$, level $1$ contains no prime by Section 5, and on level $2$ the filters leave the $36$ numbers $a00b$. Their digit sum is $a+b$, and $24$ of them are not divisible by $3$. Of these, $11$ are composite, for example $1001=7\cdot11\cdot13$ and $2009=7^2\cdot41$, and $13$ are prime:

$1009, 2003, 3001, 4001, 4003, 4007,$
$5003, 5009, 6007, 7001, 8009,$
$9001, 9007,$

with the sum $67061$. All ten rows agree with the table, and the ten sums add up to $273700$. Only $56$ of the $1061$ four-digit primes take part.

8. How many primes should we expect?

Why do one or two changes suffice, and why do some digits need two? The prime number theorem says that the number of primes up to $x$ is about $x/\ln x$, so near $x$ roughly one number in $\ln x$ is prime. ^[5] For ten-digit numbers $\ln10^{10}\approx23.0$. The true count agrees: there are $455\,052\,511$ primes below $10^{10}$ and $50\,847\,534$ below $10^9$, ^[5] so $404\,204\,977$ of the $9\cdot10^9$ ten-digit numbers are prime, about one in $22.3$.

Our candidates are not ordinary numbers: they end in $1$, $3$, $7$ or $9$, and their digit sum is not divisible by $3$. Only $4$ of every $10$ numbers end in such a digit, and only $2$ of every $3$ are not divisible by $3$. Every prime above $5$ passes both conditions, so a number that passes them is $\tfrac{10}4\cdot\tfrac32=3.75$ times as likely to be prime as a number taken at random: about $3.75/23\approx0.16$, one chance in six. This is a heuristic, not a theorem, but it describes the search well.

$d=1,3,7,9$: level $1$ has $46$ or $56$ candidates, and the sum of $3.75/\ln x$ over them is $8.2$, $9.6$, $7.6$ and $9.2$. The search finds $11$, $7$, $9$ and $8$ primes.
$d=2,4,5,6,8$: level $1$ has only two candidates, so about $0.33$ primes are expected, and both candidates are composite with probability about $(1-0.16)^2\approx0.7$. For these five digits level $1$ is a lottery.
Level $2$ has $230$ candidates for $d=2$ and for $d=8$ (Section 9), and about $40$ and $38$ primes are expected. For $d=0$ it has only $24$ candidates, and about $4$ primes are expected.

9. The search for ten digits

Level $0$ contains no prime for any digit (Section 3). On level $1$, the digit $0$ has no candidates (Section 5). The digits $1$, $3$, $7$ and $9$ find $11$, $7$, $9$ and $8$ primes, so $M(10,d)=9$ for them. Each of the digits $2$, $4$, $5$, $6$, $8$ keeps the two candidates $dd\dots d1$ and $dd\dots d7$, nine copies of $d$ followed by $1$ or $7$. For $d=4$, $5$ and $6$ exactly one of the two is prime, so $M(10,d)=9$ and $N(10,d)=1$. For $d=2$ and $d=8$ both are composite:

$$\begin{aligned}2222222221&=89\cdot24968789,\\2222222227&=419\cdot5303633,\\8888888881&=31\cdot286738351,\\8888888887&=67453\cdot131779.\end{aligned}$$

So the digits $0$, $2$ and $8$ go on to level $2$. For $d=0$ the candidates are the $36$ numbers $a\cdot10^9+b$ of Section 5; $24$ of them are not divisible by $3$, and $8$ are prime. For $d=2$ and $d=8$, the last place and one other place change. The last place takes one of $1,3,7,9$, the first place one of the eight non-zero digits other than $d$, and each of the eight middle places one of nine digits, so there are $4\cdot(8+8\cdot9)=320$ candidates. The eight unchanged digits have the digit sum $16$ or $64$, both $\equiv1\pmod3$, so the test of Section 6 removes a candidate exactly when its other changed digit $r$ and its last digit $b$ satisfy $r+b\equiv2\pmod3$. A middle place then keeps $6$, $7$, $6$, $7$ digits for $b=1,3,7,9$, $26$ in all, and the first place, where $0$ is not allowed, keeps $22$, so $8\cdot26+22=230$ candidates remain; $39$ of them are prime for $d=2$ and $32$ for $d=8$. A $0$ may stand in any of the eight middle places: $5$ and $8$ of these primes contain one, and a search that never put in a $0$ would lose them. For $d=0,1,\dots,9$ the search therefore ends with

$M(10,d)$: $8, 9, 8, 9, 9, 9, 9, 9, 8, 9$
$N(10,d)$: $8, 11, 39, 7, 1, 1, 1, 9, 32, 8$

These are $117$ primes in all, and the requested number is their sum, which we leave to the program. After the place-value filters of Section 5, only $4\cdot83+5\cdot4+36+2\cdot320$${}=1028$ numbers have to be tested for primality.

10. Proving that a ten-digit number is prime

The simplest test is trial division. If $x=ab$ with $1\lt a\le b$, then $a^2\le ab=x$, so a composite $x$ has a prime factor no larger than $\sqrt x$. For a ten-digit $x$ we have $\sqrt x\lt10^5$, and there are $9592$ primes below $10^5$, ^[5] so at most that many divisions decide whether $x$ is prime. That would be fast enough for the $1028$ candidates of the C++ program, but that program uses a test that needs only a few hundred multiplications even for a prime; the Python and Java programs call library tests instead, as the code section explains.

It starts from Fermat's little theorem: if $p$ is prime and $a$ is not divisible by $p$, then $a^{p-1}\equiv1\pmod p$. ^[6] For example $2^6=64=9\cdot7+1\equiv1\pmod7$. The example also shows why: doubling $1,2,3,4,5,6$ modulo $7$ gives $2,4,6,1,3,5$, the same six numbers in another order, so $2^6\cdot6!\equiv6!\pmod7$, and $6!$, which $7$ does not divide, can be cancelled. For any prime $p$, multiplying by $a$ only reorders the remainders $1,\dots,p-1$: neither $ai$ nor $ai-aj=a(i-j)$ with $i\ne j$ is divisible by $p$, because a prime that divides a product divides one of its factors. So if $a^{x-1}\not\equiv1\pmod x$ for some $a$ not divisible by $x$, then $x$ is certainly composite. The converse is false: $341=11\cdot31$ is composite, and yet $2^{340}\equiv1\pmod{341}$. ^[6]

The Miller–Rabin test adds one observation. Modulo a prime $p$, the only numbers whose square is $1$ are $1$ and $-1$: if $y^2\equiv1$, then $p$ divides $y^2-1=(y-1)(y+1)$ and therefore $y-1$ or $y+1$. For an odd $x$ and a base $a$ not divisible by $x$, write $x-1=2^s\,t$ with $t$ odd and look at the list

$$a^{t},\ a^{2t},\ a^{4t},\ \dots,\ a^{2^{s}t}\pmod x,$$

in which every entry is the square of the one before and the last entry is $a^{x-1}$. If $x$ is prime, the last entry is $1$ by Fermat's theorem, and the entry just before the first $1$, if there is one, must be $-1$. So a prime always passes for the base $a$: either $a^t\equiv1$, or some entry before the last is $\equiv-1\pmod x$. For the prime $13$ and $a=2$, for instance, $12=2^2\cdot3$, and the list $2^3=8$, $8^2=64\equiv12\equiv-1$, $1$ reaches $1$ right after $-1$. A composite number that passes is called a strong pseudoprime to the base $a$. ^[7]

Look at $341$ again. Here $340=2^2\cdot85$, and the list for $a=2$ is $2^{85}\equiv32$ and $32^2=1024=3\cdot341+1\equiv1$. The first $1$ follows $32$, not $340\equiv-1$, so $341$ fails and is composite. The failure even exposes a factor: $341$ divides $32^2-1=31\cdot33$ without dividing $31$ or $33$, and indeed $\gcd(31,341)=31$ and $\gcd(33,341)=11$. ^[7] Some composite numbers do pass for one base. For $2047=23\cdot89$ we have $2046=2\cdot1023$, and since $2^{11}=2048\equiv1$, also $2^{1023}=(2^{11})^{93}\equiv1\pmod{2047}$; this is the smallest strong pseudoprime to the base $2$. ^[8] The base $3$ exposes it: $3^{1023}\equiv1565$, which is neither $1$ nor $-1$, and with $s=1$ there is no further entry.

Several bases together leave no gap below a known bound: no composite number below $341\,550\,071\,728\,321$ passes for all seven bases $2,3,5,7,11,13,17$. ^[7] Ten-digit numbers lie far below this bound, so for them the test with these seven bases is not a matter of probability but a proof. The bound is sharp: $341\,550\,071\,728\,321$${}=10670053\cdot32010157$ passes for all seven bases.

The powers are computed by repeated squaring. ^[9] In binary, $85=1010101_2=64+16+4+1$. Squaring $2$ six times modulo $341$ gives $2^1,2^2,2^4,\dots,2^{64}$${}\equiv2,4,16,256,64,4,16$, and then $2^{85}\equiv16\cdot64\cdot16\cdot2\equiv32$: nine multiplications instead of $84$. For a ten-digit $x$, the number $x-1=2^st$ is below $2^{34}$, so $t$ has at most $34-s$ binary digits. The loop of the C++ program computes $a^t$ with one squaring per binary digit of $t$ and one product per digit $1$, two multiplications more than the hand calculation, which skips the last squaring and the first product by $1$. The list needs at most $s-1$ more squarings, so one base costs at most $2(34-s)+s-1\le66$ multiplications. The products themselves need care: two remainders below $10^{10}$ can have a product close to $10^{20}$, more than the largest 64-bit number, $2^{64}-1\approx1.8\cdot10^{19}$.

11. Other lengths

A full sieve over all numbers with $n$ digits, which does not use the method above, gives the values of $M(n,d)$ for $d=0,1,\dots,9$:

$n=3$: $1, 2, 2, 2, 2, 2, 2, 2, 2, 2$
$n=4$: $2, 3, 3, 3, 3, 3, 3, 3, 3, 3$
$n=5$: $3, 4, 4, 4, 4, 3, 3, 4, 4, 4$
$n=6$: $4, 5, 4, 5, 5, 5, 5, 5, 5, 5$
$n=7$: $5, 6, 5, 6, 5, 5, 5, 6, 5, 6$
$n=8$: $6, 7, 7, 7, 6, 7, 7, 7, 7, 7$
$n=9$: $7, 8, 8, 8, 8, 7, 8, 8, 8, 8$
$n=10$: $8, 9, 8, 9, 9, 9, 9, 9, 8, 9$

The pattern is the one that Sections 5 and 8 predict. For these lengths the digit $0$ always reaches $n-2$, the most that Section 5 allows, and the digits $1$, $3$, $7$, $9$ always reach $n-1$. The heuristic of Section 8 explains the second fact: their level $1$ offers several candidates for every place, so the number of candidates grows like $n$, while the chance of each one falls like $1/n$, because $\ln10^n=n\ln10$; for every length from $3$ to $300$ the expected number of primes lies between about $6.9$ and $11.4$. The digits $2$, $4$, $5$, $6$, $8$ keep only four candidates on level $1$ for every length, and they drop below $n-1$ whenever all four are composite: for the lengths $3$ to $10$ in $12$ of the $40$ cases, for the lengths $11$ to $19$ in $28$ of $45$, and each time to $n-2$.

For longer numbers the pattern of the digit $0$ breaks. The number of its candidates $a\cdot10^{n-1}+b$ stays at $36$ for every length, but each of them is prime with a chance that falls like $1/n$, and only about $39/n$ primes are expected. The first length at which all of them are composite is $n=15$, so $M(15,0)=12$, and the search has to go on to level $3$; the next such length is $n=24$. In the other direction, level $0$ can contain a prime: for $n=2$, $19$ and $23$ the repunit $R_n$ is prime (Section 3), so $M(n,1)=n$ and $N(n,1)=1$. Two digits also show a degenerate case: no two-digit prime contains a $0$, so $M(2,0)=0$, and all $21$ two-digit primes count for $d=0$. The level search of Section 2 handles all of this correctly, provided that it starts at level $0$.

How the Code Works

The C++ program ^[10] reads two optional arguments in parse_arguments: --digits=n sets the length, $10$ by default, and --skip-checkpoints skips the self-check. Unless it is skipped, run_checkpoints first computes the four-digit total with the same code and stops with an error unless it equals $273700$. Then solve adds sum_for_digit(n, d) over the ten digits and prints the result.

sum_for_digit is the level search of Section 2; its loop variable replacements is the level $k$, starting at $1$ by Section 3. For each level, choose_positions lists all sets of $k$ places in increasing order by recursion; it tries a place only while enough places remain to its right, so every set of Section 4 is produced once. A vector digits holds $n$ copies of $d$, and for each set of places assign_replacements puts every digit other than $d$ into each chosen place, again by recursion. It never restores a place, which is safe because every chosen place is overwritten before the next number is formed; only after the whole set are the chosen places reset to $d$. At the end of each branch the function applies the two filters of Section 5, a leading $0$ and a last digit that is even or $5$, builds the value digit by digit as $10v+c$, and calls is_prime. When a whole level has been searched and has yielded at least one prime, the primes are sorted, duplicates are removed, which is only a precaution because Section 4 shows that there are none, and their sum is returned; the higher levels are never generated.

is_prime is the test of Section 10. It returns false below $2$ and divides by the twelve primes from $2$ to $37$, which already settles the small examples $13$, $341$ and $2047$ of Section 10. Then it writes $x-1=2^s\,t$, keeping $t$ in a variable d that has nothing to do with the repeated digit, and runs the strong test for the bases $2,3,5,7,11,13,17$: pow_mod computes $a^t\bmod x$ by repeated squaring, and a loop squares at most $s-1$ more times, looking for $x-1$. Every product is taken by mul_mod in the 128-bit type unsigned __int128, because, as Section 10 noted, the product of two ten-digit remainders can exceed the 64-bit range.

The Python program follows the same plan with library tools. itertools.combinations lists the sets of places in increasing order, ^[11] an inner recursive function fill assigns the digits and restores each place on the way back, and sympy.isprime tests primality; its answer is definitive for numbers below $2^{64}$. ^[12] The program skips only a leading zero, so it tests $8667$ numbers instead of $1028$, and it adds sum(set(primes)) at the first level that contains a prime. Python integers have no size limit, so no overflow care is needed. The Java program enumerates the sets of places with an array pos that advances like an odometer to the next set in increasing order, and the digits with an array repl that counts through all $10^k$ digit tuples and skips those containing $d$. Its isPrime accepts $2$ and $3$, rejects other multiples of $2$ and $3$, and passes the rest to BigInteger.isProbablePrime(30), whose answer true is correct with probability above $1-2^{-30}$. ^[13] Only $847$ numbers reach that call, the results agree with the exact tests of the other two programs, and every sum fits easily into a long.

The option --digits carries three assumptions of the default length to other lengths. The search starts at level $1$, which is safe exactly when the repunit $R_n$ is composite (Section 3); the seven bases of Section 10 are a proof only below $341\,550\,071\,728\,321$, that is, for at most $14$ digits; and the candidates and the total are unsigned 64-bit numbers. For $n=2$ the option prints $3107$ instead of the correct $2865$, because it misses the prime $11$ on level $0$ and adds the seven two-digit primes with a single $1$ instead. For $3\le n\le17$ it prints the correct totals, but from $n=18$ on the total exceeds $2^{64}$ and wraps around, for $n=19$ the program would also miss the repunit prime $R_{19}$, and from $n=20$ on even a candidate can exceed $2^{64}$.

Complexity and Verification

For one digit, level $k$ costs $\binom nk\cdot9^k$ branches in the C++ and Python programs, and the search stops at the first level with a prime. For ten digits this means level $1$ for seven digits and levels $1$ and $2$ for the digits $0$, $2$ and $8$: $10\cdot90+3\cdot3645=11835$ branches in total, of which the C++ filters leave $1028$ for is_prime. The Java program runs through all $10^k$ digit tuples for each set of places, $14500$ in total. One call of is_prime costs at most twelve trial divisions and seven rounds of at most $66$ multiplications of 128-bit numbers; measured, a prime takes about $2$ to $3$ microseconds and an average odd number about $0.3$, because most composite numbers fail at the first trial division or the first base. The C++ search itself takes less than a millisecond, and the whole run, start-up included, less than a hundredth of a second; the Python program takes about $0.15$ seconds, most of it for loading SymPy, and the Java program about $0.05$ seconds, mostly start-up. The memory is tiny: C++ and Python keep the sets of places of one level, at most $45$ for ten digits, and the primes of one level, while Java keeps only its two arrays and a running sum. For other lengths the search visits $\sum_{k\le k^*}\binom nk 9^k$ branches, each with $O(n)$ digit operations and at most one primality test; while $k^*\le2$, as Sections 8 and 11 lead us to expect for moderate $n$, the number of branches grows like $n^2$. A full sieve of all ten-digit numbers, by comparison, has to look at $9\cdot10^9$ numbers and took about $16$ seconds in C++.

By default the C++ program checks itself with the four-digit total $273700$ of the statement before it computes anything else; the Python and Java programs have no such check. All three programs print the same number. For this article we checked more, with code that is not part of the published solutions:

A segmented sieve of Eratosthenes over all $n$-digit numbers, for every $n$ from $2$ to $10$, counted the copies of every digit in every prime. It reproduces the whole four-digit table of the statement, the lists of Section 11 and, for ten digits, the values $M(10,d)$ and $N(10,d)$ of Section 9; its ten-digit total is the same as the output of all three programs.
The C++ program with --digits=n gives the totals of the sieve for $3\le n\le9$ and those of an exact level search with Python integers for $11\le n\le17$; for $n=2$, $18$ and $19$ it prints the wrong values explained above.
The routine is_prime agrees with a sieve on every number below $10^8$, and with trial division on every odd ten-digit number that differs from a repdigit in at most two places; it also accepts the composite number $341\,550\,071\,728\,321$, so its bound is real.
A separate level search with a primality test that is definitive below $2^{64}$ confirmed the candidate counts and factorisations of Sections 7 to 9, the absence of the factors $11$, $41$, $271$, $9091$ on level $1$, the expected counts of Section 8, and, for every $n\le19$, the statements of Section 11, including $M(15,0)=12$; a factorisation of each of the $36$ numbers $a\cdot10^{14}+b$ confirms that none of them is prime. A separate check of the lengths $20$ to $24$, with a primality proof for one prime of each length up to $23$ and a compositeness witness for each of the $36$ numbers of length $24$, confirms that $24$ is the next such length.

Footnotes and References

The references below supply the background facts used above. The level reformulation, the candidate counts, the divisibility arguments, the estimates and the ten-digit search are derived explicitly in this article.

Project Euler 111 — Primes with Runs. The official statement defines M(n, d) as the maximum number of repeated digits d in an n-digit prime, N(n, d) as the number of such primes and S(n, d) as their sum, gives the complete table for 4-digit primes, including M(4, 1) = 3, N(4, 1) = 9, S(4, 1) = 22275 and M(4, 0) = 2, N(4, 0) = 13, states that the sum of all S(4, d) is 273700, and asks for the sum of all S(10, d).
Wikipedia — Repunit. The decimal repunit is R_n = (10^n − 1)/9; if a divides n, then R_n is divisible by R_a, so a repunit with a composite number of digits is composite. The table of factorisations gives R_10 = 11 · 41 · 271 · 9091, and the decimal repunit primes occur for n = 2, 19, 23, 317, 1031, …
Wikipedia — Binomial coefficient. The binomial coefficient C(n, k) is the number of ways to choose an unordered subset of k elements from a set of n elements, and C(n, k) = n!/(k!(n − k)!).
Wikipedia — Divisibility rule. A number is divisible by 2 when its last digit is even and by 5 when its last digit is 0 or 5; it is divisible by 3 or 9 when its digit sum is, because 10 leaves the remainder 1; and it is divisible by 11 when its alternating digit sum is, because 10 leaves the remainder −1 modulo 11.
Wikipedia — Prime number theorem. The number of primes up to x is asymptotic to x/log x, so the probability that a random integer not greater than N is prime is close to 1/log N. Its table lists π(10^5) = 9,592, π(10^9) = 50,847,534 and π(10^10) = 455,052,511.
Wikipedia — Fermat's little theorem. If p is prime and a is not divisible by p, then a^(p−1) ≡ 1 (mod p). The converse fails: 341 = 11 × 31 is a pseudoprime to the base 2, the first one found, by Sarrus in 1820.
Wikipedia — Miller–Rabin primality test. With n − 1 = 2^s · d and d odd, n is a strong probable prime to the base a if a^d ≡ 1 or a^(2^r·d) ≡ −1 (mod n) for some 0 ≤ r < s; every prime passes because 1 and −1 are the only square roots of 1 modulo a prime. The article works through n = 341 with the base 2, a pseudoprime but not a strong pseudoprime, where a greatest common divisor reveals the factor 31, and states that for n < 341,550,071,728,321 the bases 2, 3, 5, 7, 11, 13 and 17 are enough.
Wikipedia — Strong pseudoprime. The first strong pseudoprimes to the base 2 are 2047, 3277, 4033, 4681, 8321, …, so 2047 is the smallest one.
Wikipedia — Modular exponentiation. The right-to-left binary method writes the exponent in binary, squares the base once per binary digit and multiplies the result by the current square for every digit 1; its running time is proportional to the logarithm of the exponent.
C++ — Euler111.cpp. The linked immutable C++ revision contains parse_arguments, mul_mod, pow_mod, is_prime with trial division by the primes up to 37 and the Miller–Rabin bases 2 to 17, choose_positions, assign_replacements, sum_for_digit, solve and run_checkpoints. It is the implementation record for this article.
Python — itertools.combinations. The combinations of r elements are emitted in lexicographic order according to the order of the input, elements are treated as unique based on their position, and the number of combinations is n!/(r!(n − r)!).
SymPy — Number theory (isprime). The documentation of isprime states that for n < 2^64 the answer is definitive, using deterministic Miller–Rabin tests for such numbers, while larger n get a strong BPSW test with no known counterexamples.
Java — BigInteger. isProbablePrime(certainty) returns false only if the number is definitely composite; if it returns true, the probability that the number is prime exceeds 1 − 1/2^certainty.

Problem 111 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    int digits = 10;
    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, "--digits=", options.digits)) {
            continue;
        }

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

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

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

bool is_prime(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 == 0ULL) {
            return false;
        }
    }

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

    for (u64 a : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL}) {
        if (a % n == 0ULL) {
            continue;
        }
        u64 x = pow_mod(a, d, n);
        if (x == 1ULL || x == n - 1) {
            continue;
        }

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

    return true;
}

void choose_positions(const int n,
                      const int k,
                      const int start,
                      std::vector<int>& current,
                      std::vector<std::vector<int>>& all) {
    if (static_cast<int>(current.size()) == k) {
        all.push_back(current);
        return;
    }

    for (int i = start; i <= n - (k - static_cast<int>(current.size())); ++i) {
        current.push_back(i);
        choose_positions(n, k, i + 1, current, all);
        current.pop_back();
    }
}

void assign_replacements(const int idx,
                         const std::vector<int>& positions,
                         const int repeated_digit,
                         std::vector<int>& digits,
                         std::vector<u64>& primes) {
    if (idx == static_cast<int>(positions.size())) {
        if (digits[0] == 0) {
            return;
        }
        const int last = digits.back();
        if ((last % 2) == 0 || last == 5) {
            return;
        }

        u64 value = 0;
        for (int d : digits) {
            value = value * 10ULL + static_cast<u64>(d);
        }
        if (is_prime(value)) {
            primes.push_back(value);
        }
        return;
    }

    for (int d = 0; d <= 9; ++d) {
        if (d == repeated_digit) {
            continue;
        }
        digits[static_cast<std::size_t>(positions[static_cast<std::size_t>(idx)])] = d;
        assign_replacements(idx + 1, positions, repeated_digit, digits, primes);
    }
}

u64 sum_for_digit(const int n, const int repeated_digit) {
    for (int replacements = 1; replacements <= n; ++replacements) {
        std::vector<std::vector<int>> position_sets;
        std::vector<int> current;
        choose_positions(n, replacements, 0, current, position_sets);

        std::vector<u64> primes;
        std::vector<int> digits(static_cast<std::size_t>(n), repeated_digit);

        for (const std::vector<int>& positions : position_sets) {
            assign_replacements(0, positions, repeated_digit, digits, primes);
            for (int pos : positions) {
                digits[static_cast<std::size_t>(pos)] = repeated_digit;
            }
        }

        if (!primes.empty()) {
            std::sort(primes.begin(), primes.end());
            primes.erase(std::unique(primes.begin(), primes.end()), primes.end());
            u64 sum = 0;
            for (u64 p : primes) {
                sum += p;
            }
            return sum;
        }
    }

    return 0;
}

u64 solve(const int digits) {
    u64 total = 0;
    for (int d = 0; d <= 9; ++d) {
        total += sum_for_digit(digits, d);
    }
    return total;
}

bool run_checkpoints() {
    if (solve(4) != 273700ULL) {
        std::cerr << "Checkpoint failed for n=4" << '\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.digits) << '\n';
    return 0;
}

Python

# Problem 111: Primes with runs
# Sum of all 10-digit primes with maximum repeated digits.

from sympy import isprime
from itertools import combinations

def solve():
    n = 10
    total = 0
    for d in range(10):
        for replacements in range(1, n + 1):
            primes = []
            positions = list(combinations(range(n), replacements))
            for pos in positions:
                digits = [d] * n
                def fill(idx, digs):
                    if idx == len(pos):
                        if digs[0] == 0: return
                        val = int(''.join(map(str, digs)))
                        if isprime(val):
                            primes.append(val)
                        return
                    for r in range(10):
                        if r == d: continue
                        digs[pos[idx]] = r
                        fill(idx + 1, digs)
                        digs[pos[idx]] = d
                fill(0, digits)
            if primes:
                total += sum(set(primes))
                break
    print(total)

solve()

Java

import java.math.BigInteger;

public class Euler111 {
    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        if (n < 4)
            return true;
        if (n % 2 == 0 || n % 3 == 0)
            return false;
        BigInteger bi = BigInteger.valueOf(n);
        return bi.isProbablePrime(30);
    }

    public static void main(String[] args) {
        int N = 10;
        long total = 0;
        for (int d = 0; d <= 9; d++) {
            boolean found = false;
            for (int rep = 1; rep <= N && !found; rep++) {
                // rep = number of positions replaced (non-d digits)
                long sum = 0;
                boolean any = false;
                int[] pos = new int[rep];
                for (int i = 0; i < rep; i++)
                    pos[i] = i;
                // iterate over all combinations of 'rep' positions out of N
                outer: while (true) {
                    int[] digits = new int[N];
                    java.util.Arrays.fill(digits, d);
                    // try all replacement values
                    int[] repl = new int[rep];
                    inner: while (true) {
                        boolean valid = true;
                        for (int i = 0; i < rep; i++) {
                            if (repl[i] == d) {
                                valid = false;
                                break;
                            }
                            digits[pos[i]] = repl[i];
                        }
                        if (valid && digits[0] != 0) {
                            long val = 0;
                            for (int i = 0; i < N; i++)
                                val = val * 10 + digits[i];
                            if (isPrime(val)) {
                                sum += val;
                                any = true;
                            }
                        }
                        // restore
                        for (int i = 0; i < rep; i++)
                            digits[pos[i]] = d;
                        // next replacement combo
                        int carry = rep - 1;
                        while (carry >= 0) {
                            repl[carry]++;
                            if (repl[carry] < 10)
                                break;
                            repl[carry] = 0;
                            carry--;
                        }
                        if (carry < 0)
                            break;
                    }
                    // next combination
                    int k = rep - 1;
                    while (k >= 0 && pos[k] == N - rep + k)
                        k--;
                    if (k < 0)
                        break;
                    pos[k]++;
                    for (int i = k + 1; i < rep; i++)
                        pos[i] = pos[i - 1] + 1;
                }
                if (any) {
                    total += sum;
                    found = true;
                }
            }
        }
        System.out.println(total);
    }
}