Problem 773: Ruff Numbers

View on Project Euler

Project Euler Problem 773 Solution

EulerSolve provides an optimized solution for Project Euler Problem 773, Ruff Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For Problem 773, the implementations begin with the first \(k\) primes whose last digit is \(7\). The target quantity is then reduced to an inclusion-exclusion calculation in which a large multiplicative term is corrected by a small decimal-residue term. After that reduction, the whole answer depends only on the product of those primes, the totient of that product, and a short alternating binomial sum with period \(4\). Mathematical Approach Let $$p_1,p_2,\dots,p_k$$ be the first \(k\) primes satisfying $$p_i \equiv 7 \pmod{10}.$$ Define the two central products $$P_k=\prod_{i=1}^k p_i,\qquad \Phi_k=\prod_{i=1}^k (p_i-1).$$ Because the primes are distinct, \(\Phi_k\) is exactly Euler's totient of \(P_k\): $$\Phi_k=\varphi(P_k).$$ Step 1: Separate the Main Multiplicative Part The derived counting formula used by the implementations splits each inclusion-exclusion contribution into a large uniform part and a small decimal correction. The uniform part depends on which primes are excluded from a chosen subset, so the natural sum runs over all subsets \(S\subseteq \{1,\dots,k\}\): $$\sum_{S}(-1)^{|S|}\prod_{i\notin S} p_i.$$ This is a standard product expansion....

Detailed mathematical approach

Problem Summary

For Problem 773, the implementations begin with the first \(k\) primes whose last digit is \(7\). The target quantity is then reduced to an inclusion-exclusion calculation in which a large multiplicative term is corrected by a small decimal-residue term. After that reduction, the whole answer depends only on the product of those primes, the totient of that product, and a short alternating binomial sum with period \(4\).

Mathematical Approach

Let

$$p_1,p_2,\dots,p_k$$

be the first \(k\) primes satisfying

$$p_i \equiv 7 \pmod{10}.$$

Define the two central products

$$P_k=\prod_{i=1}^k p_i,\qquad \Phi_k=\prod_{i=1}^k (p_i-1).$$

Because the primes are distinct, \(\Phi_k\) is exactly Euler's totient of \(P_k\):

$$\Phi_k=\varphi(P_k).$$

Step 1: Separate the Main Multiplicative Part

The derived counting formula used by the implementations splits each inclusion-exclusion contribution into a large uniform part and a small decimal correction. The uniform part depends on which primes are excluded from a chosen subset, so the natural sum runs over all subsets \(S\subseteq \{1,\dots,k\}\):

$$\sum_{S}(-1)^{|S|}\prod_{i\notin S} p_i.$$

This is a standard product expansion. Each prime \(p_i\) contributes either \(p_i\) or \(-1\), so the whole sum factors as

$$\sum_{S}(-1)^{|S|}\prod_{i\notin S} p_i=\prod_{i=1}^k (p_i-1)=\Phi_k.$$

In the final formula this term appears multiplied by \(5\), which is why the dominant contribution is \(5\Phi_k\).

Step 2: Track the Decimal Residue Class

The problem-specific correction depends on the last digit of the product of the selected primes. If a subset has size \(t\), then every chosen prime contributes a factor congruent to \(7 \pmod{10}\), so the subset product has last digit

$$7^t \pmod{10}.$$

The powers of \(7\) modulo \(10\) are periodic with period \(4\):

$$7^0\equiv 1,\quad 7^1\equiv 7,\quad 7^2\equiv 9,\quad 7^3\equiv 3 \pmod{10}.$$

To force the relevant multiple back into the decimal class ending in \(7\), the multiplier must satisfy

$$m \equiv 7\cdot (7^t)^{-1} \pmod{10}.$$

Since the invertible residues modulo \(10\) are \(1,3,7,9\), their inverses are

$$1^{-1}\equiv 1,\qquad 7^{-1}\equiv 3,\qquad 9^{-1}\equiv 9,\qquad 3^{-1}\equiv 7 \pmod{10}.$$

Therefore the correction sequence is

$$w_0=7,\qquad w_1=1,\qquad w_2=3,\qquad w_3=9,$$

and then it repeats every four steps.

Step 3: Compress the Correction with Binomial Coefficients

The key simplification is that the correction depends only on the subset size \(t\), not on which particular primes were selected. There are exactly \(\binom{k}{t}\) subsets of size \(t\), so the full correction becomes

$$S_k=\sum_{t=0}^{k}(-1)^t\binom{k}{t}w_t,$$

where \(w_t\) is understood periodically:

$$w_t=w_{t\bmod 4}\in \{7,1,3,9\}.$$

This is why the code only needs one short period-4 table instead of handling subsets individually.

Step 4: Assemble the Closed Formula

Combining the multiplicative part and the residue correction yields the exact quantity computed by the implementations:

$$R_k \equiv P_k\left(5\Phi_k + S_k\right)\pmod{10^9+7}.$$

Substituting the explicit definitions gives

$$\boxed{R_k \equiv \left(\prod_{i=1}^k p_i\right)\left(5\prod_{i=1}^k (p_i-1)+\sum_{t=0}^{k}(-1)^t\binom{k}{t}w_t\right)\pmod{10^9+7}.}$$

So once the relevant primes are known, the remaining work is purely modular arithmetic.

Step 5: Worked Example for \(k=3\)

The first three primes ending in \(7\) are

$$7,\ 17,\ 37.$$

Hence

$$P_3=7\cdot 17\cdot 37=4403,$$

and

$$\Phi_3=(7-1)(17-1)(37-1)=6\cdot 16\cdot 36=3456.$$

The correction sum uses the pattern \(7,1,3,9\):

$$\begin{aligned} S_3&=\binom{3}{0}7-\binom{3}{1}1+\binom{3}{2}3-\binom{3}{3}9\\ &=7-3+9-9\\ &=4. \end{aligned}$$

Therefore

$$R_3=4403\left(5\cdot 3456+4\right)=4403\cdot 17284=76101452,$$

which matches the checkpoint built into the implementations.

How the Code Works

The C++, Python, and Java implementations all follow the same sequence. First they scan the arithmetic progression \(7,17,27,\dots\) and retain exactly those terms that are prime. The primality test is straightforward trial division up to \(\sqrt{n}\), which is sufficient because the required \(k\) is small.

After the prime list is built, the implementation multiplies those primes modulo \(10^9+7\) to obtain \(P_k\), and simultaneously multiplies the factors \(p_i-1\) modulo \(10^9+7\) to obtain \(\Phi_k\).

Next it evaluates the correction sum \(S_k\). Instead of recomputing each binomial coefficient from factorials, it updates them iteratively via

$$\binom{k}{t+1}=\binom{k}{t}\frac{k-t}{t+1}.$$

Division modulo the prime modulus is performed with a modular inverse:

$$a^{-1}\equiv a^{M-2}\pmod{M},\qquad M=10^9+7.$$

The alternating sign is handled term by term, the period-4 residue table supplies \(w_t\), and finally the program multiplies \(P_k\) by \(5\Phi_k+S_k\) modulo \(10^9+7\). One implementation also includes the \(k=3\) checkpoint shown above before evaluating the final \(k=97\) case.

Complexity Analysis

Let \(B\) be the largest candidate examined while collecting the first \(k\) primes ending in \(7\). With trial division, primality testing costs \(O(\sqrt{n})\) per candidate, so the prime-generation phase is the dominant part and is bounded by \(O(B^{3/2})\) in the simplest worst-case estimate. The modular accumulation of \(P_k\) and \(\Phi_k\) is \(O(k)\), and the binomial correction loop is \(O(k\log M)\) because each modular inverse is computed by fast exponentiation. Memory usage is \(O(k)\) for the stored prime list.

Footnotes and References

  1. Problem page: Project Euler 773
  2. Euler's totient function: Wikipedia - Euler's totient function
  3. Inclusion-exclusion principle: Wikipedia - Inclusion-exclusion principle
  4. Binomial coefficient: Wikipedia - Binomial coefficient
  5. Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse

Problem 773 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using i64 = long long;

constexpr i64 MOD = 1'000'000'007LL;

bool is_prime(int n) {
    if (n < 2) {
        return false;
    }
    if (n % 2 == 0) {
        return n == 2;
    }
    for (int d = 3; static_cast<i64>(d) * d <= n; d += 2) {
        if (n % d == 0) {
            return false;
        }
    }
    return true;
}

std::vector<int> first_primes_ending7(int k) {
    std::vector<int> out;
    for (int x = 7; static_cast<int>(out.size()) < k; x += 10) {
        if (is_prime(x)) {
            out.push_back(x);
        }
    }
    return out;
}

i64 mod_pow(i64 a, i64 e) {
    i64 r = 1;
    i64 cur = a % MOD;
    while (e > 0) {
        if (e & 1LL) {
            r = static_cast<i64>((static_cast<__int128>(r) * cur) % MOD);
        }
        cur = static_cast<i64>((static_cast<__int128>(cur) * cur) % MOD);
        e >>= 1LL;
    }
    return r;
}

i64 solve(int k) {
    const std::vector<int> primes = first_primes_ending7(k);

    i64 p_mod = 1;
    i64 phi_mod = 1;
    for (const int p : primes) {
        p_mod = static_cast<i64>((static_cast<__int128>(p_mod) * p) % MOD);
        phi_mod = static_cast<i64>((static_cast<__int128>(phi_mod) * (p - 1)) % MOD);
    }

    const int pattern[4] = {7, 1, 3, 9};
    i64 s = 0;
    i64 comb = 1;

    for (int t = 0; t <= k; ++t) {
        i64 term = static_cast<i64>((static_cast<__int128>(comb) * pattern[t & 3]) % MOD);
        if ((t & 1) == 0) {
            s += term;
            if (s >= MOD) {
                s -= MOD;
            }
        } else {
            s -= term;
            if (s < 0) {
                s += MOD;
            }
        }

        if (t < k) {
            const i64 num = k - t;
            const i64 inv = mod_pow(t + 1, MOD - 2);
            comb = static_cast<i64>((static_cast<__int128>(comb) * num) % MOD);
            comb = static_cast<i64>((static_cast<__int128>(comb) * inv) % MOD);
        }
    }

    i64 inner = (5LL * phi_mod + s) % MOD;
    return static_cast<i64>((static_cast<__int128>(p_mod) * inner) % MOD);
}

}  // namespace

int main() {
    assert(solve(3) == 76'101'452LL);
    std::cout << solve(97) << '\n';
    return 0;
}

Python

import math

MOD = 1000000007

def is_prime(n):
    if n < 2: return False
    if n % 2 == 0: return n == 2
    for d in range(3, int(math.isqrt(n)) + 1, 2):
        if n % d == 0: return False
    return True

def first_primes_ending7(k):
    out = []
    x = 7
    while len(out) < k:
        if is_prime(x):
            out.append(x)
        x += 10
    return out

def mod_pow(a, e):
    return pow(a, e, MOD)

def solve_k(k):
    primes = first_primes_ending7(k)
    
    p_mod = 1
    phi_mod = 1
    for p in primes:
        p_mod = (p_mod * p) % MOD
        phi_mod = (phi_mod * (p - 1)) % MOD
        
    pattern = [7, 1, 3, 9]
    s = 0
    comb = 1
    
    for t in range(k + 1):
        term = (comb * pattern[t & 3]) % MOD
        if (t & 1) == 0:
            s = (s + term) % MOD
        else:
            s = (s - term + MOD) % MOD
            
        if t < k:
            num = k - t
            inv = mod_pow(t + 1, MOD - 2)
            comb = (comb * num) % MOD
            comb = (comb * inv) % MOD
            
    inner = (5 * phi_mod + s) % MOD
    return (p_mod * inner) % MOD

def solve():
    return str(solve_k(97))

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

Java

import java.util.ArrayList;

public class Euler773 {

    static final long MOD = 1000000007L;

    static boolean isPrime(int n) {
        if (n < 2)
            return false;
        if (n % 2 == 0)
            return n == 2;
        for (int d = 3; (long) d * d <= n; d += 2) {
            if (n % d == 0)
                return false;
        }
        return true;
    }

    static int[] firstPrimesEnding7(int k) {
        ArrayList<Integer> out = new ArrayList<>();
        int x = 7;
        while (out.size() < k) {
            if (isPrime(x)) {
                out.add(x);
            }
            x += 10;
        }
        int[] result = new int[out.size()];
        for (int i = 0; i < out.size(); i++) {
            result[i] = out.get(i);
        }
        return result;
    }

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

    static long solveK(int k) {
        int[] primes = firstPrimesEnding7(k);

        long pMod = 1;
        long phiMod = 1;
        for (int p : primes) {
            pMod = (pMod * p) % MOD;
            phiMod = (phiMod * (p - 1)) % MOD;
        }

        int[] pattern = { 7, 1, 3, 9 };
        long s = 0;
        long comb = 1;

        for (int t = 0; t <= k; ++t) {
            long term = (comb * pattern[t & 3]) % MOD;
            if ((t & 1) == 0) {
                s = (s + term) % MOD;
            } else {
                s = (s - term + MOD) % MOD;
            }

            if (t < k) {
                long num = k - t;
                long inv = modPow(t + 1, MOD - 2);
                comb = (comb * num) % MOD;
                comb = (comb * inv) % MOD;
            }
        }

        long inner = (5L * phiMod + s) % MOD;
        return (pMod * inner) % MOD;
    }

    public static String solve() {
        return Long.toString(solveK(97));
    }

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