Problem 381: (prime-k) factorial

Mathematical Foundation

Theorem (Wilson’s Theorem). For any prime $p$,

Proof. In $\mathbb{Z}/p\mathbb{Z}$, every nonzero element has a unique multiplicative inverse. The only elements equal to their own inverse are $\pm 1$ (roots of $x^2 \equiv 1$). In the product $(p-1)! = 1 \cdot 2 \cdots (p-1)$, all other elements cancel in inverse pairs, leaving $(p-1)! \equiv 1 \cdot (-1) = -1 \pmod{p}$. $\square$

Lemma (Backward Factorial Recurrence). For a prime $p \ge 5$ and $1 \le k \le 5$:

where $\operatorname{inv}(a) = a^{p-2} \bmod p$ denotes the modular inverse by Fermat’s little theorem.

Proof. We use the identity $(p-k)! = (p-1)! / \prod_{j=p-k+1}^{p-1} j$ modulo $p$, combined with Wilson’s theorem.

• $k=1$: $(p-1)! \equiv -1 \pmod{p}$, so $(p-1)! \bmod p = p - 1$.

• $k=2$: $(p-2)! = (p-1)! / (p-1) \equiv (-1)/(-1) = 1 \pmod{p}$.

• $k=3$: $(p-3)! = (p-2)! / (p-2) \equiv 1 / (-2) \equiv -(p+1)/2 \pmod{p}$.

  Since $-1 \equiv p-1 \pmod{p}$, this equals $(p-1) \cdot \operatorname{inv}(2) \bmod p$.

• $k=4$: $(p-4)! = (p-3)! / (p-3) \equiv \frac{-\operatorname{inv}(2)}{-3} = \operatorname{inv}(6) \pmod{p}$.

• $k=5$: $(p-5)! = (p-4)! / (p-4) \equiv \frac{\operatorname{inv}(6)}{-4} = -\operatorname{inv}(24) \equiv (p-1)\cdot\operatorname{inv}(24) \pmod{p}$.

$\square$

Theorem (Closed Form for $S(p)$). For any prime $p \ge 5$:

which simplifies to

Concretely, $S(p) = \bigl(p - 3 \cdot \operatorname{inv}(8, p)\bigr) \bmod p$, where $\operatorname{inv}(8, p) = 8^{p-2} \bmod p$.

Proof. Multiply through by 24:

Hence $S(p) \equiv -9 \cdot \operatorname{inv}(24) = -9 \cdot \operatorname{inv}(24) \equiv -3 \cdot \operatorname{inv}(8) \pmod{p}$, since $9/24 = 3/8$. $\square$

Editorial

For prime p, S(p) = sum of (p-k)! mod p for k=1..5. Using Wilson’s theorem to derive closed forms, then sieve primes.

Pseudocode

    primes = sieve_of_eratosthenes(10^8)
    total = 0
    For each p in primes where p >= 5:
        inv8 = power_mod(8, p - 2, p)
        sp = (p - 3 * inv8 % p) % p
        total += sp
    Return total

Optimization: Instead of computing $8^{p-2} \bmod p$ via fast exponentiation for each prime, note that $8^{-1} \equiv (p+1)/8 \pmod{p}$ when $8 \mid (p+1)$, or use the extended Euclidean algorithm in $O(\log p)$ per prime. Alternatively, since $8 = 2^3$, compute $\operatorname{inv}(2) = (p+1)/2$ and cube it.

## Complexity Analysis

- **Time:** $O(N / \ln N \cdot \log p)$ for modular inverse per prime via fast exponentiation, or $O(N)$ with sieve. With the Sieve of Eratosthenes taking $O(N \log \log N)$ and $O(1)$ per prime using the algebraic inverse, total is $O(N \log \log N)$ where $N = 10^8$.
- **Space:** $O(N)$ for the prime sieve.

## Answer

$$
\boxed{139602943319822}
$$

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 381: (prime-k) factorial
 *
 * For prime p, S(p) = sum of (p-k)! mod p for k=1..5.
 * Using Wilson's theorem to derive closed forms, then sieve primes.
 *
 * Answer: 139602943319822
 */

typedef long long ll;

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    const int N = 100000000;
    vector<bool> is_prime(N, true);
    is_prime[0] = is_prime[1] = false;

    for (int i = 2; (ll)i * i < N; i++) {
        if (is_prime[i]) {
            for (int j = i * i; j < N; j += i)
                is_prime[j] = false;
        }
    }

    // For prime p >= 5:
    // (p-1)! = p-1 mod p
    // (p-2)! = 1 mod p
    // (p-3)! = (p-1) * inv(2) mod p   [since (p-2)!/(p-2) = 1/(p-2) = -inv(2) ... let's be careful]
    // Actually:
    // (p-1)! = p-1
    // (p-2)! = (p-1)!/(p-1) = (p-1)/(p-1) = 1 ... wait, (p-1) mod p = -1
    // So (p-2)! = -1 / (p-1) = -1 * (-1) = 1 mod p. Good.
    // (p-3)! = (p-2)!/(p-2) = 1/(p-2) = inv(p-2) = -inv(2) mod p = (p - (p+1)/2) ...
    // Better: inv(p-2) = inv(-2) = -inv(2) = -(p+1)/2 mod p = (p - (p+1)/2) = (p-1)/2
    // (p-4)! = (p-3)!/(p-3) = (p-1)/2 * inv(p-3) = (p-1)/2 * (-inv(3)) = -(p-1)/(2*3) = -(p-1)*inv(6)
    // (p-5)! = (p-4)!/(p-4) = [-(p-1)*inv(6)] * inv(p-4) = [-(p-1)*inv(6)] * (-inv(4))
    //        = (p-1)*inv(24) mod p = (p-1)*inv(24)

    // S(p) = (p-1) + 1 + (p-1)/2 + (-(p-1)*inv(6)) + (p-1)*inv(24) mod p
    // S(p) = p + (p-1)*(1 + 1/2 - 1/6 + 1/24) mod p
    // Wait, let me recompute:
    // S(p) = (p-1) + 1 + (p-1)*inv(2) - (p-1)*inv(6) + (p-1)*inv(24)
    // = p + (p-1) * (inv(2) - inv(6) + inv(24))
    // = p + (p-1) * (12 - 4 + 1) * inv(24)
    // = p + (p-1) * 9 * inv(24)
    // = p + (p-1) * 3 * inv(8)
    // mod p: = 0 + (-1)*3*inv(8) = -3*inv(8) mod p = p - 3*inv(8) mod p

    // inv(8) mod p:  use Fermat's little theorem or direct
    // For p=5: inv(8) = inv(3) = 2. S = 5 - 6 = -1 = 4 mod 5.
    // Check: (4!+3!+2!+1!+0!) mod 5 = (24+6+2+1+1) mod 5 = 34 mod 5 = 4. Correct!

    ll total = 0;
    for (int p = 5; p < N; p++) {
        if (!is_prime[p]) continue;
        ll pp = p;
        // Compute inv(8) mod p
        // 8 * inv8 = 1 mod p => inv8 = (p+1)/8 if 8|(p+1), else use pow
        // Use Fermat: inv8 = pow(8, p-2, p)
        // But that's slow for every prime. Better: use the formula directly.
        // S(p) = (-3 * inv(8)) mod p
        // inv(8) = inv(2)^3. inv(2) = (p+1)/2.
        ll inv2 = (pp + 1) / 2;
        ll inv8 = inv2 % pp;
        inv8 = inv8 * inv2 % pp;
        inv8 = inv8 * inv2 % pp;
        ll s = (pp - 3 * inv8 % pp) % pp;
        total += s;
    }

    cout << total << endl;

    return 0;
}

"""
Problem 381: (prime-k) factorial

For prime p, S(p) = sum of (p-k)! mod p for k=1..5.
Using Wilson's theorem to derive closed forms, then sieve primes.

Answer: 139602943319822
"""

def solve():
    """
    By Wilson's theorem, (p-1)! = -1 mod p.

    Working backwards:
    (p-1)! = p-1 mod p
    (p-2)! = 1 mod p
    (p-3)! = (p-1)/2 mod p  (i.e., -inv(2) mod p)
    (p-4)! = -(p-1)*inv(6) mod p
    (p-5)! = (p-1)*inv(24) mod p

    Summing: S(p) = -3 * inv(8) mod p

    We sieve primes up to 10^8 and compute S(p) for each.
    """
    N = 10**8
    # Sieve of Eratosthenes
    sieve = bytearray([1]) * N
    sieve[0] = sieve[1] = 0
    for i in range(2, int(N**0.5) + 1):
        if sieve[i]:
            sieve[i*i::i] = bytearray(len(sieve[i*i::i]))

    total = 0
    for p in range(5, N):
        if sieve[p]:
            inv2 = (p + 1) // 2
            inv8 = pow(inv2, 3, p)
            s = (-3 * inv8) % p
            total += s

    print(total)

if __name__ == "__main__":
    solve()