Problem 784: Reciprocal Pairs

Project Euler Problem 784 Solution

EulerSolve provides an optimized solution for Project Euler Problem 784, Reciprocal Pairs, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a given bound $n$, we sum $p+q$ over all triples of positive integers $(p,q,r)$ satisfying $$q > p > r,\qquad pr \equiv 1 \pmod q,\qquad qr \equiv 1 \pmod p,\qquad p \le n.$$ A naive search over all triples is far too expensive. The key observation is that, once $r$ is fixed, every valid pair $(p,q)$ comes from a divisor of $r^2-1$, so the problem becomes a structured divisor-enumeration task instead of a three-dimensional brute-force scan. Mathematical Approach Let $S(n)$ denote the required total: $$S(n)=\sum (p+q),$$ where the sum ranges over all triples $(p,q,r)$ satisfying the conditions above. The derivation below explains why the implementations only need to loop over $r$ and divisors of $r^2-1$. Step 1: Turn the congruences into exact equations Because $pr \equiv 1 \pmod q$, there is an integer $a$ such that $$pr-1=aq.$$ Likewise, from $qr \equiv 1 \pmod p$, there is an integer $b$ with $$qr-1=bp.$$ Since $q>r$, we have $$0 \le a=\frac{pr-1}{q} < p.$$ Reduce the first equation modulo $p$. Because $qr \equiv 1 \pmod p$, the number $q$ is the inverse of $r$ modulo $p$....

Detailed mathematical approach

Problem Summary

For a given bound $n$, we sum $p+q$ over all triples of positive integers $(p,q,r)$ satisfying

$$q > p > r,\qquad pr \equiv 1 \pmod q,\qquad qr \equiv 1 \pmod p,\qquad p \le n.$$

A naive search over all triples is far too expensive. The key observation is that, once $r$ is fixed, every valid pair $(p,q)$ comes from a divisor of $r^2-1$, so the problem becomes a structured divisor-enumeration task instead of a three-dimensional brute-force scan.

Mathematical Approach

Let $S(n)$ denote the required total:

$$S(n)=\sum (p+q),$$

where the sum ranges over all triples $(p,q,r)$ satisfying the conditions above. The derivation below explains why the implementations only need to loop over $r$ and divisors of $r^2-1$.

Step 1: Turn the congruences into exact equations

Because $pr \equiv 1 \pmod q$, there is an integer $a$ such that

$$pr-1=aq.$$

Likewise, from $qr \equiv 1 \pmod p$, there is an integer $b$ with

$$qr-1=bp.$$

Since $q>r$, we have

$$0 \le a=\frac{pr-1}{q} < p.$$

Reduce the first equation modulo $p$. Because $qr \equiv 1 \pmod p$, the number $q$ is the inverse of $r$ modulo $p$. Multiplying $pr-1=aq$ by $r$ modulo $p$ gives

$$-r \equiv aqr \equiv a \pmod p.$$

The unique integer in $[0,p-1]$ congruent to $-r$ modulo $p$ is $p-r$, so

$$a=p-r.$$

Therefore

$$pr-1=(p-r)q,\qquad q=\frac{pr-1}{p-r}.$$

By symmetry, the second congruence yields

$$qr-1=(q-r)p,\qquad p=\frac{qr-1}{q-r}.$$

Step 2: Measure how far $p$ and $q$ are from $r$

Set

$$d=p-r,\qquad e=q-r,$$

so that

$$p=r+d,\qquad q=r+e,\qquad d>0,\ e>0.$$

Insert $p=r+d$ into the formula for $q$:

$$q=\frac{(r+d)r-1}{d}=r+\frac{r^2-1}{d}.$$

Hence

$$e=\frac{r^2-1}{d},\qquad de=r^2-1.$$

This is the core reduction. For a fixed $r$, every valid triple corresponds to a factor pair of $r^2-1$.

Step 3: Determine the exact divisor range

The condition $p \le n$ becomes

$$r+d \le n,\qquad d \le n-r.$$

The ordering $q>p$ is equivalent to $e>d$. Since $de=r^2-1$, this means

$$d^2 < r^2-1.$$

For integers $r \ge 2$, that is equivalent to

$$d \le r-1.$$

So the admissible divisors are exactly those satisfying

$$1 \le d \le m_r,\qquad m_r=\min(r-1,n-r),\qquad d \mid (r^2-1).$$

Once such a divisor is chosen, the corresponding pair is forced:

$$p=r+d,\qquad q=r+\frac{r^2-1}{d}.$$

Step 4: Prove the converse direction

The reduction is not just necessary; it is exact. Suppose $r \ge 2$ and $d$ is a divisor of $r^2-1$ with $1 \le d \le m_r$. Define

$$e=\frac{r^2-1}{d},\qquad p=r+d,\qquad q=r+e.$$

Then

$$pr-1=r(r+d)-1=r^2+rd-1=d(r+e)=dq=(p-r)q,$$

so $pr \equiv 1 \pmod q$. Similarly,

$$qr-1=r(r+e)-1=r^2+re-1=e(r+d)=ep=(q-r)p,$$

so $qr \equiv 1 \pmod p$. The bound $d \le n-r$ guarantees $p \le n$, and $d \le r-1$ guarantees $q>p$. Therefore every admissible divisor produces one valid triple, and every valid triple produces one admissible divisor. This is a bijection.

Step 5: Final summation formula

For a fixed $r$, each valid divisor contributes

$$p+q=(r+d)+\left(r+\frac{r^2-1}{d}\right)=2r+d+\frac{r^2-1}{d}.$$

Thus

$$S(n)=\sum_{r=2}^{n-1}\ \sum_{d \mid (r^2-1),\ 1 \le d \le m_r}\left(2r+d+\frac{r^2-1}{d}\right),$$

where $m_r=\min(r-1,n-r)$. This is exactly the quantity computed by the implementations.

Worked Example: $n=5$

Only $r=2,3,4$ can contribute.

For $r=2$, we have $r^2-1=3$ and $m_r=\min(1,3)=1$. The only admissible divisor is $d=1$, giving

$$p=2+1=3,\qquad q=2+\frac{3}{1}=5,$$

so the contribution is $3+5=8$.

For $r=3$, we have $r^2-1=8$ and $m_r=\min(2,2)=2$. The admissible divisors are $d=1,2$:

$$d=1 \Rightarrow (p,q)=(4,11),\quad p+q=15,$$

$$d=2 \Rightarrow (p,q)=(5,7),\quad p+q=12.$$

For $r=4$, we have $r^2-1=15$ and $m_r=\min(3,1)=1$. Only $d=1$ is allowed, producing

$$ (p,q)=(5,19),\quad p+q=24. $$

Adding everything gives

$$S(5)=8+15+12+24=59.$$

How the Code Works

The C++, Python, and Java implementations all follow the same structure. They begin by building a smallest-prime-factor table up to slightly above $n$. That allows fast factorizations of $r-1$ and $r+1$ for every $r$.

For each $r$ from $2$ to $n-1$, the implementation computes $m_r=\min(r-1,n-r)$. If $m_r \le 0$, that $r$ cannot contribute. Otherwise, it factors $r-1$ and $r+1$, merges the prime exponents, and thereby obtains the prime factorization of

$$r^2-1=(r-1)(r+1).$$

Next, it recursively enumerates divisors from that factorization. During this search, branches are pruned as soon as the partial divisor already exceeds $m_r$, because such a divisor cannot yield a valid $p$.

Whenever a complete admissible divisor $d$ is produced, the implementation computes

$$p=r+d,\qquad q=r+\frac{r^2-1}{d},$$

and adds $p+q$ to the running total. The whole algorithm is therefore a careful divisor-generation procedure driven by the exact bijection proved above.

Complexity Analysis

Building the smallest-prime-factor table costs $O(n\log\log n)$ time and $O(n)$ memory. For each $r$, factoring $r-1$ and $r+1$ is fast because the factorization follows the precomputed table. The dominant work is enumerating admissible divisors of $r^2-1$, so the total running time is well described by

$$O\!\left(n\log\log n+\sum_{r=2}^{n-1}\tau(r^2-1)\right),$$

up to the usual small-factor savings from pruning divisors larger than $m_r$. In practice, this is dramatically smaller than brute-force enumeration over all candidate triples.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=784
Modular arithmetic: Wikipedia — Modular arithmetic
Divisor function: Wikipedia — Divisor function
Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic
Sieve of Eratosthenes: Wikipedia — Sieve of Eratosthenes

Problem 784 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <functional>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

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

static std::vector<int> build_spf(int n) {
    std::vector<int> spf(n + 1);
    for (int i = 0; i <= n; ++i) {
        spf[i] = i;
    }
    for (int i = 2; 1LL * i * i <= n; ++i) {
        if (spf[i] != i) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            if (spf[j] == j) {
                spf[j] = i;
            }
        }
    }
    return spf;
}

static std::vector<std::pair<int, int>> factorize_small(int x, const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> out;
    while (x > 1) {
        const int p = spf[x];
        int c = 0;
        while (x % p == 0) {
            x /= p;
            ++c;
        }
        out.push_back({p, c});
    }
    return out;
}

static std::vector<std::pair<int, int>> merge_factorizations(
    const std::vector<std::pair<int, int>>& a,
    const std::vector<std::pair<int, int>>& b
) {
    std::vector<std::pair<int, int>> f = a;
    for (const auto& [p, c] : b) {
        bool found = false;
        for (auto& [q, d] : f) {
            if (q == p) {
                d += c;
                found = true;
                break;
            }
        }
        if (!found) {
            f.push_back({p, c});
        }
    }
    return f;
}

static std::string to_string_u128(u128 x) {
    if (x == 0) {
        return "0";
    }
    std::string s;
    while (x > 0) {
        const int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

static u128 solve_fast(int n) {
    const auto spf = build_spf(n + 1);
    u128 ans = 0;

    for (int r = 2; r < n; ++r) {
        const int m = std::min(r - 1, n - r);
        if (m <= 0) {
            continue;
        }

        const auto fa = factorize_small(r - 1, spf);
        const auto fb = factorize_small(r + 1, spf);
        auto factors = merge_factorizations(fa, fb);

        const u64 sq_minus_1 = static_cast<u64>(r) * static_cast<u64>(r) - 1ULL;

        std::function<void(int, u64)> dfs = [&](int idx, u64 cur) {
            if (idx == static_cast<int>(factors.size())) {
                if (cur == 0 || cur > static_cast<u64>(m)) {
                    return;
                }
                const u64 d = cur;
                const u64 p = static_cast<u64>(r) + d;
                const u64 q = static_cast<u64>(r) + sq_minus_1 / d;
                ans += static_cast<u128>(p) + static_cast<u128>(q);
                return;
            }

            const auto [prime, exp] = factors[idx];
            u64 val = cur;
            for (int e = 0; e <= exp; ++e) {
                if (val > static_cast<u64>(m)) {
                    break;
                }
                dfs(idx + 1, val);
                if (e == exp || val > static_cast<u64>(m) / static_cast<u64>(prime)) {
                    break;
                }
                val *= static_cast<u64>(prime);
            }
        };

        dfs(0, 1);
    }

    return ans;
}

static u128 solve_bruteforce(int n) {
    u128 ans = 0;
    for (int p = 2; p <= n; ++p) {
        for (int r = 1; r < p; ++r) {
            const int den = p - r;
            const int num = p * r - 1;
            if (num % den != 0) {
                continue;
            }
            const int q = num / den;
            if (q <= p) {
                continue;
            }
            if ((static_cast<long long>(p) * r) % q != 1) {
                continue;
            }
            if ((static_cast<long long>(q) * r) % p != 1) {
                continue;
            }
            ans += static_cast<u128>(p + q);
        }
    }
    return ans;
}

int main() {
    assert(to_string_u128(solve_fast(5)) == "59");
    assert(to_string_u128(solve_fast(100)) == "697317");
    assert(solve_fast(250) == solve_bruteforce(250));

    const u128 ans = solve_fast(2'000'000);
    std::cout << to_string_u128(ans) << '\n';
    return 0;
}

Python

def solve():
    n = 2000000
    spf = list(range(n+2))
    for i in range(2, int((n+1)**0.5)+1):
        if spf[i] == i:
            for j in range(i*i, n+2, i):
                if spf[j] == j: spf[j] = i

    def factorize(x):
        fac = []
        while x > 1:
            p = spf[x]; c = 0
            while x % p == 0: x //= p; c += 1
            fac.append((p, c))
        return fac

    ans = 0
    for r in range(2, n):
        m = min(r-1, n-r)
        if m <= 0: continue
        fa = factorize(r-1); fb = factorize(r+1)
        fac = dict(fa)
        for p, c in fb: fac[p] = fac.get(p, 0) + c
        factors = list(fac.items())
        sq1 = r*r - 1
        def dfs(idx, cur):
            nonlocal ans
            if idx == len(factors):
                if cur == 0 or cur > m: return
                d = cur; p_ = r + d; q_ = r + sq1 // d
                ans += p_ + q_; return
            pr, exp = factors[idx]; val = cur
            for e in range(exp+1):
                if val > m: break
                dfs(idx+1, val)
                if e == exp or val > m // pr: break
                val *= pr
        dfs(0, 1)

    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.HashMap;

public class Euler784 {

    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; i++) {
            spf[i] = i;
        }
        for (int i = 2; i * i <= n; i++) {
            if (spf[i] == i) {
                for (int j = i * i; j <= n; j += i) {
                    if (spf[j] == j) {
                        spf[j] = i;
                    }
                }
            }
        }
        return spf;
    }

    static class Factor {
        int p, c;

        Factor(int p, int c) {
            this.p = p;
            this.c = c;
        }
    }

    static ArrayList<Factor> factorizeSmall(int x, int[] spf) {
        ArrayList<Factor> factors = new ArrayList<>();
        while (x > 1) {
            int p = spf[x];
            int c = 0;
            while (x % p == 0) {
                x /= p;
                c++;
            }
            factors.add(new Factor(p, c));
        }
        return factors;
    }

    static ArrayList<Factor> mergeFactorizations(ArrayList<Factor> a, ArrayList<Factor> b) {
        HashMap<Integer, Integer> merged = new HashMap<>();
        for (Factor f : a) {
            merged.put(f.p, merged.getOrDefault(f.p, 0) + f.c);
        }
        for (Factor f : b) {
            merged.put(f.p, merged.getOrDefault(f.p, 0) + f.c);
        }
        ArrayList<Factor> res = new ArrayList<>();
        for (HashMap.Entry<Integer, Integer> entry : merged.entrySet()) {
            res.add(new Factor(entry.getKey(), entry.getValue()));
        }
        return res;
    }

    static long solveFast(int n) {
        int[] spf = buildSpf(n + 1);
        long ans = 0;

        for (int r = 2; r < n; r++) {
            int m = Math.min(r - 1, n - r);
            if (m <= 0)
                continue;

            ArrayList<Factor> fa = factorizeSmall(r - 1, spf);
            ArrayList<Factor> fb = factorizeSmall(r + 1, spf);
            ArrayList<Factor> factors = mergeFactorizations(fa, fb);

            long sqMinus1 = (long) r * r - 1L;
            ans += dfs(0, 1L, factors, m, r, sqMinus1);
        }

        return ans;
    }

    static long dfs(int idx, long cur, ArrayList<Factor> factors, int m, long r, long sqMinus1) {
        if (idx == factors.size()) {
            if (cur == 0 || cur > m)
                return 0;
            long p = r + cur;
            long q = r + sqMinus1 / cur;
            return p + q;
        }

        Factor f = factors.get(idx);
        long val = cur;
        long sum = 0;

        for (int e = 0; e <= f.c; e++) {
            if (val > m)
                break;
            sum += dfs(idx + 1, val, factors, m, r, sqMinus1);
            if (e == f.c || val > m / f.p)
                break;
            val *= f.p;
        }

        return sum;
    }

    public static String solve() {
        return Long.toString(solveFast(2000000));
    }

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