Problem 134: Prime Pair Connection

Project Euler Problem 134 Solution

EulerSolve provides an optimized solution for Project Euler Problem 134, Prime Pair Connection, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each pair of consecutive primes $p_1, p_2$ with $5 \le p_1 \le 10^6$, let $S(p_1,p_2)$ be the smallest positive integer that is divisible by $p_2$ and whose decimal expansion ends with the digits of $p_1$. The task is to sum $S(p_1,p_2)$ over all such consecutive prime pairs. A brute-force scan through multiples of $p_2$ would be wasteful. The suffix condition is really a congruence modulo a power of ten, so every prime pair can be solved by one modular inverse and one multiplication. Mathematical Approach Let $d$ be the number of decimal digits of $p_1$, and set $M = 10^d$. Saying that a number ends with the digits of $p_1$ is equivalent to saying that it is congruent to $p_1$ modulo $M$. Write the unknown number as a multiple of the next prime If $n$ is divisible by $p_2$, then $n = p_2 k$ for some positive integer $k$. The suffix requirement becomes $$p_2 k \equiv p_1 \pmod{M}.$$ This is the central congruence used by the implementations. It solves directly for the multiplier $k$, and then reconstructs $n$ as $p_2 k$. Why a modular inverse exists Because $p_1 \ge 5$ and $p_2$ is the next prime, we always have $p_2 \ge 7$. Therefore $p_2$ is not divisible by 2 or 5, while $M = 10^d$ has no prime factors other than 2 and 5....

Detailed mathematical approach

Problem Summary

For two consecutive primes $p_1\lt p_2$, let $S$ be the smallest positive integer that is divisible by $p_2$ and whose decimal representation ends with the digits of $p_1$. For $p_1=19$ and $p_2=23$ the statement gives $S=1219$, and it notes that such a number exists for every pair of consecutive primes except $p_1=3$, $p_2=5$. We must find the sum of $S$ over all pairs of consecutive primes with $5\le p_1\le1\,000\,000$. ^[1]

The whole solution in one sentence: ending with the digits of $p_1$ means leaving the remainder $p_1$ after division by a power of ten $M$, and since $p_2$ has no factor in common with $10$, the smallest multiple of $p_2$ with this remainder is computed in one step with a modular inverse, a number $u$ for which $p_2u$ leaves the remainder $1$ after division by $M$.

Mathematical Approach

1. Small pairs by hand

For $(p_1,p_2)=(5,7)$ the multiples $7,14,21,28,35$ of $7$ first end in $5$ at $S=35$, and for $(7,11)$ the first multiple of $11$ ending in $7$ is $S=77$. For $(11,13)$ it is quicker to test the numbers that end in $11$, namely $11,111,211,\dots$: their remainders after division by $13$ are $11,7,3,12,8,4,0$, so $S=611=13\cdot47$. In the same way $(13,17)$, $(17,19)$, $(19,23)$ and $(23,29)$ give $S=1513=17\cdot89$, $817=19\cdot43$, $1219=23\cdot53$ and $2523=29\cdot87$.

Such searches are slow for large primes: going through the multiples of $p_2$ takes up to $999981$ steps for one pair and about $3.5\cdot10^{10}$ steps for all pairs, so we look for a formula.

2. An ending is a remainder

Let $d$ be the number of decimal digits of $p_1$ and $M=10^d$, so that $10^{d-1}\le p_1\lt M$; for $p_1=19$, $M=100$. Division by $M$ leaves the last $d$ digits of a number as the remainder, as in $1219=12\cdot100+19$. As $p_1$ has exactly $d$ digits, a positive integer $n$ ends with the digits of $p_1$ exactly when this remainder is $p_1$; with $1000$ instead of $100$, the remainder $19$ would mean the ending $019$. In the notation of modular arithmetic, $n\equiv p_1\pmod M$: the difference $n-p_1$ is a multiple of $M$. ^[2] The candidates are therefore $p_1$, $p_1+M$, $p_1+2M$, and so on, and $S$ is the smallest positive solution of

$$\begin{aligned}n&\equiv p_1\pmod M,\\n&\equiv0\pmod{p_2}.\end{aligned}$$

3. Why only the pair $(3,5)$ fails

The pair $(2,3)$ works, with $S=12$, but $(3,5)$ fails because $5$ divides $10$: every number that ends in $3$ leaves the remainder $3$ after division by $5$. For $p_1\ge5$ there is no such obstacle: $p_2\ge7$ is a prime other than $2$ and $5$ and $M=2^d5^d$, so the greatest common divisor $\gcd(p_2,M)$ is $1$.

Now take the $p_2$ numbers $p_1+tM$ with $t=0,1,\dots,p_2-1$. If two of them, with $t\ne t'$, left the same remainder after division by $p_2$, the prime $p_2$ would divide $(t-t')M$, and so, by Euclid's lemma, ^[3] it would divide $M$ or $t-t'$. Neither is possible, as $\gcd(p_2,M)=1$ and $0\lt|t-t'|\lt p_2$. So the remainders are all different, each of $0,1,\dots,p_2-1$ occurs exactly once, and exactly one of the numbers is a multiple of $p_2$; for $(11,13)$ the remainders are $11,7,3,12,8,4,0,9,5,1,10,6,2$. Hence $S$ exists, and

$$S\le p_1+(p_2-1)M\lt p_2M.$$

This is the case of two moduli of the Chinese remainder theorem: the two conditions of Section 2 have exactly one solution with $0\le n\lt Mp_2$. ^[4]

4. Dividing by $p_2$ modulo $M$

The positive multiples of $p_2$ are the numbers $n=p_2k$ with $k\ge1$, and the ending condition becomes

$$p_2k\equiv p_1\pmod M.$$

Modulo $M$, dividing by $p_2$ means multiplying by an inverse, an integer $u$ with $p_2u\equiv1\pmod M$; for example $23\cdot87=2001$, so $87$ is an inverse of $23$ modulo $100$. Since $\gcd(p_2,M)=1$, an inverse exists and is unique modulo $M$. ^[5] If $p_2k\equiv p_1$, multiplying by $u$ gives $k\equiv up_2k\equiv up_1\pmod M$; conversely, every $k\equiv up_1$ satisfies $p_2k\equiv p_2up_1\equiv p_1$. So the multipliers that work are exactly the integers congruent modulo $M$ to

$$k_0=p_1u\bmod M.$$

5. The smallest multiplier gives $S$

The remainder $k_0$ lies between $0$ and $M-1$, and it is not $0$, since otherwise $p_1\equiv p_2k_0\equiv0\pmod M$, which is impossible for $0\lt p_1\lt M$. So the positive multipliers that work are $k_0$, $k_0+M$, $k_0+2M$, and so on, and since $p_2k$ grows with $k$, the smallest valid number is

$$S=p_2k_0=p_2\cdot\bigl(p_1u\bmod M\bigr).$$

Only $p_2\bmod M$ matters for $u$: it differs from $p_2$ by a multiple of $M$, so it has the same inverses.

6. The extended Euclidean algorithm

Euclid's algorithm divides with remainder again and again, each time dividing the previous divisor by the remainder:

$$\begin{aligned}100&=4\cdot23+8,\\23&=2\cdot8+7,\\8&=1\cdot7+1,\\7&=7\cdot1+0.\end{aligned}$$

A division $a=qb+r$ keeps the common divisors, since a number divides $a$ and $b$ exactly when it divides $b$ and $r$; so the last nonzero remainder, $1$, is $\gcd(100,23)$. Going back up, with $7=23-2\cdot8$ and then $8=100-4\cdot23$, writes it as a combination of the two numbers:

$$\begin{aligned}1&=8-7=3\cdot8-23\\&=3\cdot(100-4\cdot23)-23\\&=3\cdot100-13\cdot23.\end{aligned}$$

Modulo $100$ this says $23\cdot(-13)\equiv1$, so $-13$, or $-13+100=87$, is an inverse of $23$. In general the extended Euclidean algorithm finds integers $x,y$ with $ax+My=\gcd(a,M)$, Bézout's identity; when the gcd is $1$, $x$ reduced modulo $M$ is the inverse of $a$. ^[6] By Lamé's theorem there are at most five times as many division steps as the smaller number has digits, ^[3] at most $30$ here; no pair of the problem needs more than $23$.

7. Worked example

For $(19,23)$: $M=100$, $u=87$ by Section 6, $p_1u=19\cdot87=1653$, so $k_0=53$ and $S=23\cdot53=1219$, the value of the statement.

If $p_2$ has more digits than $p_1$, the reduction modulo $M$ of Section 5 becomes visible. For $(997,1009)$: $1009\bmod1000=9$ and $9\cdot889=8001$, so $u=889$; then $997\cdot889=886333$ gives $k_0=333$ and $S=1009\cdot333=335997$. Six pairs of the problem are of this kind, with $p_1=7$, $97$, $997$, $9973$, $99991$ and $999983$.

8. A second formula as a check

Swapping the roles of the two conditions, write $n=p_1+tM$ and require $tM\equiv-p_1\pmod{p_2}$. With an inverse $v$ of $M$ modulo $p_2$, the smallest $t\ge0$ is $t_0=(-p_1v)\bmod p_2$, and $S=p_1+t_0M$. For $(19,23)$, $100\bmod23=8$ and $8\cdot3=24\equiv1$, so $v=3$, and $-19\cdot3=-57\equiv12\pmod{23}$ gives $S=1219$ again. Both formulas give the unique solution below $p_2M$ of Section 3, so they always agree; the second serves as a check below.

9. Primes, pairs and the size of the sum

The primes come from the sieve of Eratosthenes. Start with the integers from $2$ to a bound $L$; for each $q$ not yet crossed out, in increasing order while $q^2\le L$, cross out $q^2,q^2+q,q^2+2q,\dots$. Smaller multiples of $q$ have a smaller prime factor and are already crossed out, and every composite number up to $L$ has a prime factor at most $\sqrt L$, so exactly the primes remain. ^[7] The last pair is $(999983,1000003)$, so $L$ must be at least $1000003$; the programs use $L=1\,200\,000$. The prime-counting function gives $\pi(10^6)=78498$ primes up to one million, ^[8] so without $2$ and $3$ there are $78496$ pairs. Every $S$ is below $p_2M\le1000003\cdot10^6$, so the sum is below $8\cdot10^{16}$, far below $2^{63}\approx9.2\cdot10^{18}$: signed 64-bit integers suffice.

How the Code Works

The C++ program ^[9] reads its options in parse_arguments: --limit= followed by digits, read into an int without an overflow check, with default $1\,000\,000$ and at least $5$, and --skip-checkpoints. Unless that flag is given, run_checkpoints checks that pair_value(19, 23) is $1219$ and otherwise ends the program with exit code $2$. solve sieves with primes_up_to(limit + 200000), a std::vector<bool> version of Section 9, walks through the neighbouring entries primes[i] and primes[i + 1], skips $p_1\lt5$, stops at the first $p_1$ that is at least limit and adds pair_value(p1, p2) to a 64-bit sum. So it sums over the $p_1$ below limit: for the composite $10^6$ exactly the range of the statement, while a prime limit such as --limit=7 would leave out the pair $(7,11)$.

pair_value is Sections 2 to 5: it builds $M=10^d$ by multiplying by $10$ once per digit of $p_1$, takes $u$ from mod_inverse(p2 % mod, mod), forms $k_0$ as (p1 * inv) % mod and returns $p_2k_0$. mod_inverse calls the recursive extended_gcd(a, b, x, y), which returns $a$ with $x=1$, $y=0$ when $b=0$, and otherwise solves $(b,a\bmod b)$ with coefficients $x_1,y_1$ and sets $x=y_1$, $y=x_1-\lfloor a/b\rfloor y_1$, the back-substitution of Section 6. For $23$ and $100$ it returns $x=-13$, and because x % mod keeps the sign in C++, mod_inverse adds $M$ and returns $87$. It returns $-1$ when the gcd is not $1$, which Section 3 rules out. The modular arithmetic and the sum use std::int64_t; the product $p_1u$ is below $M^2\le10^{12}$.

The Python program does the same with integers of unlimited size, taking the primes from SymPy's sieve.primerange(5, limit + 200000), the primes $p$ with $5\le p\lt1\,200\,000$, and the inverse from the built-in pow(p2 % mod, -1, mod) of Python 3.8 or later. The Java program sieves a boolean[] of length limit + 200001, keeps the primes in an ArrayList<Integer>, inverts with BigInteger.modInverse and sums in a long. ^[10] Both have the limit built in and no checkpoint.

Complexity and Verification

The sieve up to $L=1\,200\,000$ takes $O(L\log\log L)$ time and $O(L)$ memory for its flags and the list of primes. ^[7] Each of the $78496$ pairs needs at most $6$ multiplications by $10$ and one modular inverse. In C++ that is an extended Euclidean algorithm with at most $23$ division steps, so the recursion is at most $25$ calls deep: the first call only swaps its arguments and the last one has $b=0$; Python and Java use library routines. Counting an arithmetic operation on numbers of this size as one step, and as $\pi(L)$ grows like $L/\log L$, ^[8] the pairs cost $O(\pi(L)\log L)=O(L)$ and the whole run $O(L\log\log L)$, a few million steps. On a desktop computer each program finished in well under a second.

Besides the C++ checkpoint, the following was checked for this article, outside the published solutions. All three programs print the same sum. The brute force of Section 1 over the multiples of $p_2$ agrees with pair_value for all $78496$ pairs, and so does the formula of Section 8 with $v=M^{p_2-2}\bmod p_2$, computed without the Euclidean algorithm; it is an inverse of $M$ because Euler's theorem gives $M^{p_2-1}\equiv1\pmod{p_2}$ for the prime $p_2$. ^[5] Decimal-string tests agree for all $1227$ pairs with $p_1\le10^4$. For testing other programs: the pairs with $5\le p_1\le100$ give the sum $69155$, and those with $5\le p_1\le1000$ give $36941222$, which --limit=100 and --limit=1000 reproduce.

Footnotes and References

The references supply the background facts; the derivations are given in this article.

Project Euler 134 — Prime Pair Connection. The official statement gives S = 1219 for p₁ = 19 and p₂ = 23, notes that only the pair p₁ = 3, p₂ = 5 has no such number, and asks for the sum of S over all pairs of consecutive primes with 5 ≤ p₁ ≤ 1000000.
Wikipedia — Modular arithmetic. Integers are congruent modulo m when their difference is a multiple of m, equivalently when they leave the same remainder after division by m.
Wikipedia — Euclidean algorithm. It states Bézout's identity, Euclid's lemma (a prime that divides a product divides one of its factors) and Lamé's theorem (at most five times as many division steps as the smaller integer has decimal digits).
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N, a system of congruences has exactly one solution x with 0 ≤ x < N.
Wikipedia — Modular multiplicative inverse. An inverse of a modulo m exists if and only if gcd(a, m) = 1 and is unique modulo m; it is found with the extended Euclidean algorithm, and for a prime m Euler's theorem gives the inverse a^(m − 2).
Wikipedia — Extended Euclidean algorithm. It computes integers x and y with ax + by = gcd(a, b); for coprime a and b, x is the inverse of a modulo b.
Wikipedia — Sieve of Eratosthenes. Multiples of each prime p are crossed out from p² on, the sieve can stop once p² exceeds the bound n, and it takes O(n log log n) time and O(n) memory.
Wikipedia — Prime-counting function. Its table gives π(10⁶) = 78498, and the prime number theorem says that π(x) / (x / log x) tends to 1.
C++ — Euler134.cpp. The linked immutable revision, with pair_value, extended_gcd and the checkpoint for (19, 23), is the implementation record for this article.
Python — pow: since Python 3.8, pow(base, −1, mod) returns an inverse of base modulo mod. SymPy — Ntheory functions: Sieve.primerange(a, b) generates the primes in [a, b). Java — BigInteger.modInverse returns this⁻¹ mod m and throws an ArithmeticException if none exists.

Problem 134 source code

C++

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

namespace {

using i64 = std::int64_t;

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

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

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

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

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

i64 extended_gcd(i64 a, i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 x1 = 0;
    i64 y1 = 0;
    const i64 g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

i64 mod_inverse(i64 a, i64 mod) {
    i64 x = 0;
    i64 y = 0;
    const i64 g = extended_gcd(a, mod, x, y);
    if (g != 1) {
        return -1;
    }
    x %= mod;
    if (x < 0) {
        x += mod;
    }
    return x;
}

i64 pair_value(const int p1, const int p2) {
    i64 mod = 1;
    int x = p1;
    while (x > 0) {
        mod *= 10;
        x /= 10;
    }

    const i64 inv = mod_inverse(p2 % mod, mod);
    const i64 k = (static_cast<i64>(p1) * inv) % mod;
    return static_cast<i64>(p2) * k;
}

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

    for (int i = 2; static_cast<i64>(i) * i <= limit; ++i) {
        if (!sieve[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            sieve[static_cast<std::size_t>(j)] = false;
        }
    }

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

i64 solve(const int limit) {
    const std::vector<int> primes = primes_up_to(limit + 200000);

    i64 sum = 0;
    for (std::size_t i = 0; i + 1 < primes.size(); ++i) {
        const int p1 = primes[i];
        const int p2 = primes[i + 1];
        if (p1 < 5) {
            continue;
        }
        if (p1 >= limit) {
            break;
        }
        sum += pair_value(p1, p2);
    }

    return sum;
}

bool run_checkpoints() {
    if (pair_value(19, 23) != 1219) {
        std::cerr << "Checkpoint failed for pair (19,23)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

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

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

Python

# Problem 134: Prime pair connection
# For consecutive primes p1, p2 (5 <= p1 <= 10^6), find S = smallest number
# ending in p1 that is divisible by p2. Sum all S values.

from sympy import sieve

def solve():
    limit = 1000000
    primes = list(sieve.primerange(5, limit + 200000))
    total = 0
    for i in range(len(primes) - 1):
        p1 = primes[i]
        p2 = primes[i + 1]
        if p1 >= limit: break
        mod = 1
        x = p1
        while x > 0: mod *= 10; x //= 10
        inv = pow(p2 % mod, -1, mod)
        k = (p1 * inv) % mod
        total += p2 * k
    print(total)

solve()

Java

import java.math.BigInteger;
import java.util.*;

public class Euler134 {
    public static void main(String[] args) {
        int limit = 1000000;
        boolean[] sieve = new boolean[limit + 200001];
        Arrays.fill(sieve, true);
        sieve[0] = sieve[1] = false;
        for (int i = 2; (long) i * i < sieve.length; i++)
            if (sieve[i])
                for (int j = i * i; j < sieve.length; j += i)
                    sieve[j] = false;
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i < sieve.length; i++)
            if (sieve[i])
                primes.add(i);

        long total = 0;
        for (int i = 0; i + 1 < primes.size(); i++) {
            int p1 = primes.get(i), p2 = primes.get(i + 1);
            if (p1 < 5)
                continue;
            if (p1 >= limit)
                break;
            long mod = 1;
            int x = p1;
            while (x > 0) {
                mod *= 10;
                x /= 10;
            }
            long inv = BigInteger.valueOf(p2 % mod).modInverse(BigInteger.valueOf(mod)).longValue();
            long k = (p1 * inv) % mod;
            total += (long) p2 * k;
        }
        System.out.println(total);
    }
}