Problem 789: Minimal Pairing Modulo $p$

Project Euler Problem 789 Solution

EulerSolve provides an optimized solution for Project Euler Problem 789, Minimal Pairing Modulo $p$, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For an odd prime $p$, pair the numbers $1,2,\dots,p-1$ into $\frac{p-1}{2}$ disjoint pairs. Each pair $\{a,b\}$ contributes the least positive residue of $ab \bmod p$, and the goal is to minimize the total sum over all pairs. The implementations solve the large target prime by converting this pairing problem into a modular optimization problem inside $(\mathbb{Z}/p\mathbb{Z})^\times$, then reconstructing the corresponding witness integer. Mathematical Approach Let $n=\frac{p-1}{2}$. If a pairing produces pair residues $r_1,\dots,r_n\in\{1,\dots,p-1\}$, then the total pairing cost is $$S(p)=\min \sum_{i=1}^{n} r_i.$$ Step 1: Every valid pairing gives a product constraint If the pairs are $\{a_i,b_i\}$, then $r_i\equiv a_i b_i \pmod p$. Because the pairing uses each nonzero residue exactly once, we get $$\prod_{i=1}^{n} r_i \equiv \prod_{k=1}^{p-1} k \equiv -1 \pmod p,$$ where the last congruence is Wilson's theorem. So every admissible pairing yields a multiset of positive residues whose product is $-1$ modulo $p$. Step 2: Split the total into baseline plus excess Every pair contributes at least $1$, so $$\sum_{i=1}^{n} r_i = n + \sum_{i=1}^{n}(r_i-1).$$ The constant term $n$ is unavoidable. The real optimization is the excess $\sum (r_i-1)$....

Detailed mathematical approach

Problem Summary

For an odd prime $p$, pair the numbers $1,2,\dots,p-1$ into $\frac{p-1}{2}$ disjoint pairs. Each pair $\{a,b\}$ contributes the least positive residue of $ab \bmod p$, and the goal is to minimize the total sum over all pairs. The implementations solve the large target prime by converting this pairing problem into a modular optimization problem inside $(\mathbb{Z}/p\mathbb{Z})^\times$, then reconstructing the corresponding witness integer.

Mathematical Approach

Let $n=\frac{p-1}{2}$. If a pairing produces pair residues $r_1,\dots,r_n\in\{1,\dots,p-1\}$, then the total pairing cost is

$$S(p)=\min \sum_{i=1}^{n} r_i.$$

Step 1: Every valid pairing gives a product constraint

If the pairs are $\{a_i,b_i\}$, then $r_i\equiv a_i b_i \pmod p$. Because the pairing uses each nonzero residue exactly once, we get

$$\prod_{i=1}^{n} r_i \equiv \prod_{k=1}^{p-1} k \equiv -1 \pmod p,$$

where the last congruence is Wilson's theorem. So every admissible pairing yields a multiset of positive residues whose product is $-1$ modulo $p$.

Step 2: Split the total into baseline plus excess

Every pair contributes at least $1$, so

$$\sum_{i=1}^{n} r_i = n + \sum_{i=1}^{n}(r_i-1).$$

The constant term $n$ is unavoidable. The real optimization is the excess $\sum (r_i-1)$.

If a residue $c>1$ is composite, say $c=uv$, then

$$c-1=(u-1)+(v-1)+(u-1)(v-1)\ge (u-1)+(v-1).$$

This is why prime factors are the natural atomic pieces: splitting a composite residue into smaller factors never increases the excess. The implementations therefore search for a witness of the form

$$W=\prod_q q^{e_q}\equiv -1 \pmod p,$$

with $q$ prime and $e_q\ge 0$, minimizing

$$E(p)=\sum_q e_q(q-1).$$

For the large target prime, this excess is tiny compared with $n$, so the factorized search matches the pairing objective cleanly. The small-prime checks in the C++ implementation confirm the relation

$$S(p)=\frac{p-1}{2}+E(p).$$

Step 3: Compress the search into residue-cost states

Fix a current excess budget $B$. Any prime with $q-1>B$ cannot appear, so only primes $q\le B+1$ matter. An exponent vector determines two values:

$$r=\prod_q q^{e_q}\bmod p,\qquad c=\sum_q e_q(q-1).$$

For each reachable residue $r$, only the smallest cost $c$ matters. Many exponent vectors can therefore be compressed into one best state per residue class.

Step 4: Use meet-in-the-middle

Split the allowed primes into two groups. Enumerate all reachable residue-cost states on the left and on the right, always staying within the same budget $B$. If the left side produces $r_A$, then the right side must produce

$$r_B\equiv -r_A^{-1}\pmod p$$

so that $r_A r_B\equiv -1\pmod p$. Because $p$ is prime, every nonzero residue has an inverse. The cheapest compatible left-right combination gives the best excess under budget $B$.

The implementations further speed up the merge step by computing many inverses in one batch, using prefix products and one Fermat-style inverse instead of inverting every left residue separately.

Step 5: Reconstruct the exact witness

After the smallest feasible excess has been found, a second exact search recovers the exponents $e_q$ that attain it. Their product

$$W=\prod_q q^{e_q}$$

is then built as an arbitrary-precision integer. This witness is what the implementations finally print for the large prime from the problem.

Worked Example: $p=7$

Here $n=\frac{7-1}{2}=3$, and the target congruence is $W\equiv -1\equiv 6 \pmod 7$.

The cheapest prime-factor witness is

$$W=2\cdot 3=6,$$

with excess

$$E(7)=(2-1)+(3-1)=3.$$

So the predicted minimum pairing cost is

$$S(7)=3+3=6.$$

An explicit pairing is

$$\{1,3\},\qquad \{2,4\},\qquad \{5,6\},$$

whose pair residues are $3$, $1$, and $2$. Their sum is $3+1+2=6$, and their product is $3\cdot 1\cdot 2=6\equiv -1\pmod 7$, exactly as the theory predicts. The $p=5$ checkpoint works the same way: $W=4=2^2$, the excess is $2$, and the optimal pairing sum is $4$.

How the Code Works

The C++, Python, and Java implementations start with a small excess budget and enlarge it gradually until a feasible witness appears. For each budget they generate all primes up to $B+1$, split them into two groups, and recursively enumerate every exponent pattern whose weighted cost stays within the budget. Each side stores only the cheapest cost for each residue modulo $p$.

After the two maps are built, the implementation searches for matching residues whose product is $-1$ modulo $p$. When the first successful budget is found, a second recursive pass reconstructs the exact exponents and multiplies the corresponding primes to obtain the final big-integer witness $W$. The C++ implementation also checks small primes by exhaustive pairing search, confirming that the computed excess matches the true minimum pairing total after adding the baseline $\frac{p-1}{2}$.

Complexity Analysis

Let $B^*$ be the minimal excess cost. For a given budget $B$, sieving primes up to $B+1$ costs $O(B\log\log B)$. The dominant work is enumerating exponent vectors satisfying

$$\sum_q e_q(q-1)\le B.$$

If $N_A(B)$ and $N_B(B)$ are the numbers of states explored on the two sides, then one round costs

$$O\bigl(B\log\log B + N_A(B) + N_B(B)\bigr)$$

time. Memory usage is

$$O\bigl(M_A(B)+M_B(B)\bigr),$$

where $M_A(B)$ and $M_B(B)$ are the numbers of stored residues in the two hash maps. In practice the budget stays small, so the running time depends much more on the optimal excess than on the huge size of $p$.

Footnotes and References

Problem page: https://projecteuler.net/problem=789
Wilson's theorem: Wikipedia — Wilson's theorem
Multiplicative group modulo $n$: Wikipedia — Multiplicative group of integers modulo n
Fermat's little theorem: Wikipedia — Fermat's little theorem
Meet-in-the-middle: Wikipedia — Meet-in-the-middle
Matching in graph theory: Wikipedia — Matching

Problem 789 source code

C++

#include <algorithm>
#include <boost/multiprecision/cpp_int.hpp>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <unordered_map>
#include <utility>
#include <vector>

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = long long;
using boost::multiprecision::cpp_int;

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

static u32 mod_pow(u32 a, u64 e, u32 mod) {
    u32 r = 1;
    while (e > 0) {
        if (e & 1ULL) {
            r = mul_mod(r, a, mod);
        }
        a = mul_mod(a, a, mod);
        e >>= 1ULL;
    }
    return r;
}

static std::vector<int> primes_upto(int n) {
    std::vector<bool> sieve(n + 1, true);
    sieve[0] = sieve[1] = false;
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (!sieve[i]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            sieve[j] = false;
        }
    }

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

static void enumerate_min_costs(
    const std::vector<int>& primes,
    int idx,
    int cost,
    u32 residue,
    int limit,
    u32 mod,
    std::unordered_map<u32, int>& best
) {
    if (idx == static_cast<int>(primes.size())) {
        auto it = best.find(residue);
        if (it == best.end() || cost < it->second) {
            best[residue] = cost;
        }
        return;
    }

    const int q = primes[idx];
    const int w = q - 1;
    const int max_e = (limit - cost) / w;

    u32 r = residue;
    for (int e = 0; e <= max_e; ++e) {
        enumerate_min_costs(primes, idx + 1, cost + e * w, r, limit, mod, best);
        r = mul_mod(r, static_cast<u32>(q), mod);
    }
}

static bool find_combo_exact(
    const std::vector<int>& primes,
    int idx,
    int rem_cost,
    u32 residue,
    u32 target,
    u32 mod,
    std::vector<int>& exps
) {
    if (idx == static_cast<int>(primes.size())) {
        return rem_cost == 0 && residue == target;
    }

    const int q = primes[idx];
    const int w = q - 1;
    const int max_e = rem_cost / w;

    u32 r = residue;
    for (int e = 0; e <= max_e; ++e) {
        exps[idx] = e;
        if (find_combo_exact(primes, idx + 1, rem_cost - e * w, r, target, mod, exps)) {
            return true;
        }
        r = mul_mod(r, static_cast<u32>(q), mod);
    }

    return false;
}

struct SolveResult {
    int extra_cost;
    cpp_int product;
};

static SolveResult solve_prime(u32 p) {
    int limit = 80;

    while (true) {
        const int pmax = std::min<int>(static_cast<int>(p) - 1, limit + 1);
        auto primes = primes_upto(pmax);
        const int split = std::min<int>(4, static_cast<int>(primes.size()));

        std::vector<int> pa(primes.begin(), primes.begin() + split);
        std::vector<int> pb(primes.begin() + split, primes.end());

        std::unordered_map<u32, int> bestA;
        std::unordered_map<u32, int> bestB;

        if (p > 1'000'000U) {
            bestA.reserve(5'000'000);
            bestB.reserve(2'000'000);
        }

        enumerate_min_costs(pa, 0, 0, 1, limit, p, bestA);
        enumerate_min_costs(pb, 0, 0, 1, limit, p, bestB);

        const u32 target = p - 1;

        std::vector<u32> keys;
        std::vector<int> costsA;
        keys.reserve(bestA.size());
        costsA.reserve(bestA.size());
        for (const auto& [res, c] : bestA) {
            keys.push_back(res);
            costsA.push_back(c);
        }

        std::vector<u32> prefix(keys.size() + 1, 1);
        for (std::size_t i = 0; i < keys.size(); ++i) {
            prefix[i + 1] = mul_mod(prefix[i], keys[i], p);
        }
        u32 inv_total = mod_pow(prefix.back(), static_cast<u64>(p) - 2ULL, p);

        std::vector<u32> invs(keys.size(), 1);
        for (std::size_t i = keys.size(); i-- > 0;) {
            invs[i] = mul_mod(inv_total, prefix[i], p);
            inv_total = mul_mod(inv_total, keys[i], p);
        }

        int best_total = std::numeric_limits<int>::max();
        int best_ca = -1;
        int best_cb = -1;
        u32 best_ra = 0;
        u32 best_rb = 0;

        for (std::size_t i = 0; i < keys.size(); ++i) {
            const u32 ra = keys[i];
            const int ca = costsA[i];
            const u32 need = mul_mod(target, invs[i], p);
            const auto it = bestB.find(need);
            if (it == bestB.end()) {
                continue;
            }

            const int cb = it->second;
            const int total = ca + cb;
            if (total < best_total) {
                best_total = total;
                best_ca = ca;
                best_cb = cb;
                best_ra = ra;
                best_rb = need;
            }
        }

        if (best_total != std::numeric_limits<int>::max() && best_total <= limit) {
            std::vector<int> expsA(pa.size(), 0);
            std::vector<int> expsB(pb.size(), 0);

            const bool okA = find_combo_exact(pa, 0, best_ca, 1, best_ra, p, expsA);
            const bool okB = find_combo_exact(pb, 0, best_cb, 1, best_rb, p, expsB);
            if (!okA || !okB) {
                throw std::runtime_error("reconstruction failed");
            }

            cpp_int prod = 1;
            int cost_check = 0;
            u32 residue_check = 1;

            for (int i = 0; i < static_cast<int>(pa.size()); ++i) {
                const int q = pa[i];
                const int e = expsA[i];
                for (int j = 0; j < e; ++j) {
                    prod *= q;
                    residue_check = mul_mod(residue_check, static_cast<u32>(q), p);
                    cost_check += q - 1;
                }
            }
            for (int i = 0; i < static_cast<int>(pb.size()); ++i) {
                const int q = pb[i];
                const int e = expsB[i];
                for (int j = 0; j < e; ++j) {
                    prod *= q;
                    residue_check = mul_mod(residue_check, static_cast<u32>(q), p);
                    cost_check += q - 1;
                }
            }

            if (cost_check != best_total || residue_check != target) {
                throw std::runtime_error("consistency check failed");
            }

            return {best_total, prod};
        }

        limit += 40;
        if (limit > 600) {
            throw std::runtime_error("failed to find solution within limit");
        }
    }
}

static int brute_min_total_cost(int p) {
    const int n = p - 1;
    const int full = (1 << n) - 1;

    std::vector<int> memo(1 << n, -1);

    std::function<int(int)> dfs = [&](int mask) -> int {
        if (mask == 0) {
            return 0;
        }
        int& ret = memo[mask];
        if (ret != -1) {
            return ret;
        }

        int i = 0;
        while (((mask >> i) & 1) == 0) {
            ++i;
        }

        int best = std::numeric_limits<int>::max();
        const int rem = mask ^ (1 << i);
        for (int j = i + 1; j < n; ++j) {
            if (((rem >> j) & 1) == 0) {
                continue;
            }
            const int c = static_cast<int>((1LL * (i + 1) * (j + 1)) % p);
            best = std::min(best, c + dfs(rem ^ (1 << j)));
        }

        ret = best;
        return ret;
    };

    return dfs(full);
}

int main() {
    {
        const auto r5 = solve_prime(5);
        assert(r5.extra_cost == 2);
        assert(r5.product == 4);
    }

    for (int p : {7, 11, 13, 17, 19}) {
        const auto r = solve_prime(static_cast<u32>(p));
        const int brute = brute_min_total_cost(p);
        assert(brute == r.extra_cost + (p - 1) / 2);
    }

    const auto ans = solve_prime(2'000'000'011u);
    std::cout << ans.product << '\n';
    return 0;
}

Python

def solve():
    p = 2000000011

    def mul_mod(a, b): return a * b % p
    def mod_pow(a, e):
        r = 1
        while e > 0:
            if e & 1: r = r * a % p
            a = a * a % p; e >>= 1
        return r

    def primes_upto(n):
        sieve = bytearray(b'\x01')*(n+1); sieve[0] = sieve[1] = 0
        for i in range(2, int(n**0.5)+1):
            if sieve[i]:
                for j in range(i*i, n+1, i): sieve[j] = 0
        return [i for i in range(2, n+1) if sieve[i]]

    limit = 80
    while True:
        pmax = min(p-1, limit+1)
        primes = primes_upto(pmax)
        split = min(4, len(primes))
        pa, pb = primes[:split], primes[split:]

        def enum_min(prs, idx, cost, res, lim, best):
            if idx == len(prs):
                if res not in best or cost < best[res]: best[res] = cost
                return
            q = prs[idx]; w = q-1; me = (lim-cost)//w; r = res
            for e in range(me+1):
                enum_min(prs, idx+1, cost+e*w, r, lim, best)
                r = r * q % p

        bestA = {}; bestB = {}
        enum_min(pa, 0, 0, 1, limit, bestA)
        enum_min(pb, 0, 0, 1, limit, bestB)

        target = p - 1
        # Batch invert keys of bestA
        keys = list(bestA.keys()); costsA = [bestA[k] for k in keys]
        if not keys: limit += 40; continue
        prefix = [1]*(len(keys)+1)
        for i in range(len(keys)): prefix[i+1] = prefix[i] * keys[i] % p
        inv_total = mod_pow(prefix[-1], p-2)
        invs = [0]*len(keys)
        for i in range(len(keys)-1, -1, -1):
            invs[i] = inv_total * prefix[i] % p
            inv_total = inv_total * keys[i] % p

        bt = None; bca = bcb = -1; bra = brb = 0
        for i in range(len(keys)):
            need = target * invs[i] % p
            if need in bestB:
                cb = bestB[need]; tot = costsA[i]+cb
                if bt is None or tot < bt:
                    bt = tot; bca = costsA[i]; bcb = cb; bra = keys[i]; brb = need

        if bt is not None and bt <= limit:
            def find_exact(prs, idx, rc, res, tgt, exps):
                if idx == len(prs): return rc == 0 and res == tgt
                q = prs[idx]; w = q-1; me = rc//w; r = res
                for e in range(me+1):
                    exps[idx] = e
                    if find_exact(prs, idx+1, rc-e*w, r, tgt, exps): return True
                    r = r * q % p
                return False

            exA = [0]*len(pa); exB = [0]*len(pb)
            find_exact(pa, 0, bca, 1, bra, exA)
            find_exact(pb, 0, bcb, 1, brb, exB)

            prod = 1
            for i in range(len(pa)):
                for _ in range(exA[i]): prod *= pa[i]
            for i in range(len(pb)):
                for _ in range(exB[i]): prod *= pb[i]
            return str(prod)

        limit += 40
        if limit > 600: raise ValueError("No solution found")

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;

public class Euler789 {

    static ArrayList<Integer> primesUpto(int n) {
        boolean[] sieve = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            sieve[i] = true;

        for (int i = 2; i * i <= n; i++) {
            if (sieve[i]) {
                for (int j = i * i; j <= n; j += i) {
                    sieve[j] = false;
                }
            }
        }

        ArrayList<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (sieve[i])
                primes.add(i);
        }
        return primes;
    }

    static void enumerateMinCosts(
            ArrayList<Integer> primes, int idx, int cost, long residue,
            int limit, long mod, HashMap<Long, Integer> best) {
        if (idx == primes.size()) {
            Integer currentCost = best.get(residue);
            if (currentCost == null || cost < currentCost) {
                best.put(residue, cost);
            }
            return;
        }

        int q = primes.get(idx);
        int w = q - 1;
        int maxE = (limit - cost) / w;

        long r = residue;
        for (int e = 0; e <= maxE; e++) {
            enumerateMinCosts(primes, idx + 1, cost + e * w, r, limit, mod, best);
            r = (r * q) % mod;
        }
    }

    static boolean findComboExact(
            ArrayList<Integer> primes, int idx, int remCost, long residue,
            long target, long mod, int[] exps) {
        if (idx == primes.size()) {
            return remCost == 0 && residue == target;
        }

        int q = primes.get(idx);
        int w = q - 1;
        int maxE = remCost / w;

        long r = residue;
        for (int e = 0; e <= maxE; e++) {
            exps[idx] = e;
            if (findComboExact(primes, idx + 1, remCost - e * w, r, target, mod, exps)) {
                return true;
            }
            r = (r * q) % mod;
        }

        return false;
    }

    static long modPow(long a, long e, long mod) {
        long r = 1;
        while (e > 0) {
            if ((e & 1) == 1) {
                r = (r * a) % mod;
            }
            a = (a * a) % mod;
            e >>= 1;
        }
        return r;
    }

    static BigInteger solvePrime(long p) {
        int limit = 80;

        while (true) {
            int pmax = (int) Math.min(p - 1, limit + 1);
            ArrayList<Integer> primes = primesUpto(pmax);
            int split = Math.min(4, primes.size());

            ArrayList<Integer> pa = new ArrayList<>(primes.subList(0, split));
            ArrayList<Integer> pb = new ArrayList<>(primes.subList(split, primes.size()));

            HashMap<Long, Integer> bestA = new HashMap<>();
            HashMap<Long, Integer> bestB = new HashMap<>();

            enumerateMinCosts(pa, 0, 0, 1L, limit, p, bestA);
            enumerateMinCosts(pb, 0, 0, 1L, limit, p, bestB);

            long target = p - 1;

            long[] keys = new long[bestA.size()];
            int[] costsA = new int[bestA.size()];
            int idx = 0;
            for (Map.Entry<Long, Integer> entry : bestA.entrySet()) {
                keys[idx] = entry.getKey();
                costsA[idx] = entry.getValue();
                idx++;
            }

            long[] prefix = new long[keys.length + 1];
            prefix[0] = 1;
            for (int i = 0; i < keys.length; i++) {
                prefix[i + 1] = (prefix[i] * keys[i]) % p;
            }

            long invTotal = modPow(prefix[keys.length], p - 2, p);
            long[] invs = new long[keys.length];
            for (int i = keys.length - 1; i >= 0; i--) {
                invs[i] = (invTotal * prefix[i]) % p;
                invTotal = (invTotal * keys[i]) % p;
            }

            int bestTotal = Integer.MAX_VALUE;
            int bestCa = -1, bestCb = -1;
            long bestRa = 0, bestRb = 0;

            for (int i = 0; i < keys.length; i++) {
                long ra = keys[i];
                int ca = costsA[i];
                long need = (target * invs[i]) % p;

                Integer cb = bestB.get(need);
                if (cb == null)
                    continue;

                int total = ca + cb;
                if (total < bestTotal) {
                    bestTotal = total;
                    bestCa = ca;
                    bestCb = cb;
                    bestRa = ra;
                    bestRb = need;
                }
            }

            if (bestTotal != Integer.MAX_VALUE && bestTotal <= limit) {
                int[] expsA = new int[pa.size()];
                int[] expsB = new int[pb.size()];

                boolean okA = findComboExact(pa, 0, bestCa, 1L, bestRa, p, expsA);
                boolean okB = findComboExact(pb, 0, bestCb, 1L, bestRb, p, expsB);

                BigInteger prod = BigInteger.ONE;
                for (int i = 0; i < pa.size(); i++) {
                    BigInteger bq = BigInteger.valueOf(pa.get(i));
                    prod = prod.multiply(bq.pow(expsA[i]));
                }
                for (int i = 0; i < pb.size(); i++) {
                    BigInteger bq = BigInteger.valueOf(pb.get(i));
                    prod = prod.multiply(bq.pow(expsB[i]));
                }

                return prod;
            }

            limit += 40;
            if (limit > 600) {
                throw new RuntimeException("failed to find solution within limit");
            }
        }
    }

    public static String solve() {
        return solvePrime(2000000011L).toString();
    }

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