Problem 330: Euler's Number

Project Euler Problem 330 Solution

EulerSolve provides an optimized solution for Project Euler Problem 330, Euler's Number, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The sequence $a(n)$ is defined by $$a(n)=\sum_{i=1}^{\infty}\frac{a(n-i)}{i!}\qquad(n\ge 0),$$ together with the boundary rule $$a(n)=1\qquad(n<0).$$ The problem states that every term can be written in the form $$a(n)=\frac{A(n)e+B(n)}{n!},$$ where $A(n)$ and $B(n)$ are integers, and asks for $$A(10^9)+B(10^9)\pmod{77{,}777{,}777}.$$ The modulus factors as $$77{,}777{,}777=7\cdot 11\cdot 73\cdot 101\cdot 137.$$ Mathematical Approach 1) Split the infinite sum into a finite part plus an $e$-tail. For fixed $n\ge 0$, the terms with $i\le n$ refer to already-defined values $a(n-i)$, while the terms with $i>n$ hit the boundary region $n-i<0$, where $a(n-i)=1$. Therefore $$a(n)=\sum_{i=1}^{n}\frac{a(n-i)}{i!}+\sum_{i=n+1}^{\infty}\frac{1}{i!}.$$ Using $e=\sum_{i=0}^{\infty}1/i!$, the tail becomes $$\sum_{i=n+1}^{\infty}\frac{1}{i!}=e-\sum_{k=0}^{n}\frac{1}{k!}.$$ So the recurrence is really an inhomogeneous linear relation with one $e$-term and one rational term. 2) Insert the decomposition $a(n)=\frac{A(n)e+B(n)}{n!}$. Write $A_n=A(n)$, $B_n=B(n)$....

Detailed mathematical approach

Problem Summary

The sequence $a(n)$ is defined by

$$a(n)=\sum_{i=1}^{\infty}\frac{a(n-i)}{i!}\qquad(n\ge 0),$$

together with the boundary rule

$$a(n)=1\qquad(n<0).$$

The problem states that every term can be written in the form

$$a(n)=\frac{A(n)e+B(n)}{n!},$$

where $A(n)$ and $B(n)$ are integers, and asks for

$$A(10^9)+B(10^9)\pmod{77{,}777{,}777}.$$

The modulus factors as

$$77{,}777{,}777=7\cdot 11\cdot 73\cdot 101\cdot 137.$$

Mathematical Approach

1) Split the infinite sum into a finite part plus an $e$-tail.

For fixed $n\ge 0$, the terms with $i\le n$ refer to already-defined values $a(n-i)$, while the terms with $i>n$ hit the boundary region $n-i<0$, where $a(n-i)=1$. Therefore

$$a(n)=\sum_{i=1}^{n}\frac{a(n-i)}{i!}+\sum_{i=n+1}^{\infty}\frac{1}{i!}.$$

Using $e=\sum_{i=0}^{\infty}1/i!$, the tail becomes

$$\sum_{i=n+1}^{\infty}\frac{1}{i!}=e-\sum_{k=0}^{n}\frac{1}{k!}.$$

So the recurrence is really an inhomogeneous linear relation with one $e$-term and one rational term.

2) Insert the decomposition $a(n)=\frac{A(n)e+B(n)}{n!}$.

Write $A_n=A(n)$, $B_n=B(n)$. For $k=n-i$, the finite part becomes

$$\sum_{i=1}^{n}\frac{a(n-i)}{i!} =\sum_{k=0}^{n-1}\frac{A_k e+B_k}{k!(n-k)!}.$$

Multiply the whole identity by $n!$:

$$A_n e+B_n =\sum_{k=0}^{n-1}\binom{n}{k}(A_k e+B_k)+e\,n!-\sum_{k=0}^{n}\frac{n!}{k!}.$$

Now the coefficient of $e$ and the constant part must match separately.

3) This gives two exact integer recurrences.

Comparing the coefficient of $e$:

$$A_n=n!+\sum_{k=0}^{n-1}\binom{n}{k}A_k.$$

Comparing the constant term:

$$B_n=-\sum_{k=0}^{n}\frac{n!}{k!}+\sum_{k=0}^{n-1}\binom{n}{k}B_k.$$

These are exactly the two formulas implemented in the checkpoint routine compute_exact_ab.

4) Small values show the pattern.

From the recurrences:

$$A_0=1,\qquad B_0=-1,$$

$$A_1=2,\qquad B_1=-3,$$

$$A_2=7,\qquad B_2=-12.$$

$$a(0)=e-1,\qquad a(1)=2e-3,\qquad a(2)=\frac{7e-12}{2!},$$

which matches the problem statement. The code also checks the much larger checkpoint

$$A_{10}=328161643,\qquad B_{10}=-652694486.$$

5) The target $A_n+B_n$ collapses to a much simpler expression.

Define

$$C_n=A_n+B_n.$$

Add the two recurrences:

$$C_n=n!-\sum_{k=0}^{n}\frac{n!}{k!}+\sum_{k=0}^{n-1}\binom{n}{k}C_k.$$

Now use the induction hypothesis $C_k=k!-A_k$. Then

$$\sum_{k=0}^{n-1}\binom{n}{k}C_k =\sum_{k=0}^{n-1}\frac{n!}{(n-k)!}-\sum_{k=0}^{n-1}\binom{n}{k}A_k =\sum_{j=1}^{n}\frac{n!}{j!}-\sum_{k=0}^{n-1}\binom{n}{k}A_k.$$

The two factorial sums cancel, leaving

$$C_n=-\sum_{k=0}^{n-1}\binom{n}{k}A_k=n!-A_n.$$

So the quantity asked by the problem is

$$A_n+B_n=n!-A_n.$$

6) Modulo a prime $p$, the factorial term disappears after $n\ge p$.

For each prime factor $p\in\{7,11,73,101,137\}$, Wilson-style reasoning is not even needed: once $n\ge p$, the product $n!$ contains $p$, hence

$$n!\equiv 0\pmod p.$$

Therefore, for large $n$,

$$A_n+B_n\equiv -A_n\pmod p.$$

The code keeps the unified formula

$$r_p(n)\equiv n!-A_{n'}\pmod p,$$

which also handles the tiny case $n<p$.

7) The key compression is eventual periodicity modulo each prime.

The implementation exploits the modular fact that, for each prime $p$, the sequence $A_n\bmod p$ becomes periodic from index $p$ onward with period

$$\pi_p=p(p-1).$$

So instead of working at $n=10^9$, the code reduces the index to

$$n'= \begin{cases} n, & n<p,\\ p+\bigl((n-p)\bmod p(p-1)\bigr), & n\ge p. \end{cases}$$

This is the entire reason the huge input becomes easy. The program also spot-checks this periodicity numerically for each prime factor on one full band of sample offsets.

8) Computing $A_n \bmod p$ up to the reduced index is straightforward.

For a fixed prime $p$, the code computes all values

$$A_0,A_1,\dots,A_{n'}\pmod p$$

sequentially from the recurrence

$$A_n\equiv n!+\sum_{k=0}^{n-1}\binom{n}{k}A_k\pmod p.$$

It keeps one rolling row of Pascal's triangle in the array comb, updating the binomial coefficients in place from right to left. At the same time it updates $n!\bmod p$, which instantly becomes $0$ once $n\ge p$.

9) Recombine the five prime residues by CRT.

For $n=10^9$, the per-prime residues found by the code are

$$r_7=1,\qquad r_{11}=3,\qquad r_{73}=66,\qquad r_{101}=44,\qquad r_{137}=117.$$

The Chinese remainder theorem then reconstructs the unique residue modulo

$$M=77{,}777{,}777.$$

Algorithm

1) For each prime factor $p$, reduce $n$ to $n'$ using the period $p(p-1)$ after index $p$.

2) Compute $A_0,\dots,A_{n'}\pmod p$ by the binomial-convolution recurrence.

3) Return the prime residue $r_p(n)\equiv n!-A_{n'}\pmod p$.

4) Combine the five residues with the Chinese remainder theorem.

Complexity Analysis

For a fixed prime $p$, the reduced index satisfies $n'<p^2$, and the direct convolution computes each $A_n$ by summing over all smaller $k$. So the cost per prime is

$$O((n')^2),$$

which is still tiny here because the largest reduced index is only on the order of $10^4$. The memory cost is

$$O(n').$$

Since there are only five prime factors, the whole computation is comfortably fast.

Checks And Final Result

The code verifies

$$A_{10}=328161643,\qquad B_{10}=-652694486,$$

and also checks for $0\le n\le 12$ that the modular solver agrees with the exact integer recurrence.

After CRT recombination, the final answer is

$$\boxed{15955822}.$$

Problem 330 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

constexpr u64 kDefaultN = 1'000'000'000ULL;
constexpr int kModulus = 77'777'777;
constexpr int kPrimeFactors[5] = {7, 11, 73, 101, 137};

struct Options {
    u64 n = kDefaultN;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.n = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

int normalize_mod(i64 x, const int mod) {
    const i64 r = x % static_cast<i64>(mod);
    return static_cast<int>(r < 0 ? r + mod : r);
}

int pow_mod(int base, int exp, const int mod) {
    i64 result = 1;
    i64 b = normalize_mod(base, mod);

    while (exp > 0) {
        if ((exp & 1) != 0) {
            result = (result * b) % mod;
        }
        b = (b * b) % mod;
        exp >>= 1;
    }

    return static_cast<int>(result);
}

int inverse_mod_prime(const int a, const int p) {
    // Fermat inverse, valid because p is prime and gcd(a, p)=1.
    return pow_mod(a, p - 2, p);
}

int factorial_mod_small(const u64 n, const int p) {
    int f = 1 % p;
    for (u64 i = 1ULL; i <= n; ++i) {
        f = static_cast<int>((static_cast<i64>(f) * static_cast<i64>(i % static_cast<u64>(p))) % p);
    }
    return f;
}

std::vector<int> compute_A_mod_prime(const int p, const u64 max_n) {
    std::vector<int> values(static_cast<std::size_t>(max_n) + 1ULL, 0);
    std::vector<int> comb(static_cast<std::size_t>(max_n) + 1ULL, 0);
    comb[0] = 1;

    int fact = 1 % p;

    for (u64 n = 0ULL; n <= max_n; ++n) {
        if (n > 0ULL) {
            comb[static_cast<std::size_t>(n)] = 1;
            for (u64 k = n - 1ULL; k >= 1ULL; --k) {
                const std::size_t idx = static_cast<std::size_t>(k);
                comb[idx] += comb[idx - 1ULL];
                if (comb[idx] >= p) {
                    comb[idx] -= p;
                }
            }

            if (n < static_cast<u64>(p)) {
                fact = static_cast<int>((static_cast<i64>(fact) * static_cast<i64>(n)) % p);
            } else {
                fact = 0;
            }
        }

        i64 current = fact;
        for (u64 k = 0ULL; k < n; ++k) {
            current += static_cast<i64>(comb[static_cast<std::size_t>(k)]) *
                       static_cast<i64>(values[static_cast<std::size_t>(k)]);
            current %= p;
        }

        values[static_cast<std::size_t>(n)] = static_cast<int>(current);
    }

    return values;
}

u64 reduced_index_for_prime(const u64 n, const int p) {
    if (n < static_cast<u64>(p)) {
        return n;
    }

    const u64 period = static_cast<u64>(p) * static_cast<u64>(p - 1);
    return static_cast<u64>(p) + ((n - static_cast<u64>(p)) % period);
}

int solve_mod_prime(const int p, const u64 n) {
    const u64 index = reduced_index_for_prime(n, p);
    const std::vector<int> A = compute_A_mod_prime(p, index);

    const int n_factorial_mod = (n >= static_cast<u64>(p)) ? 0 : factorial_mod_small(n, p);
    return normalize_mod(static_cast<i64>(n_factorial_mod) - static_cast<i64>(A[static_cast<std::size_t>(index)]), p);
}

int crt_combine(const std::vector<int>& residues, const std::vector<int>& moduli) {
    i64 x = 0;
    i64 current_modulus = 1;

    for (std::size_t i = 0; i < residues.size(); ++i) {
        const int m = moduli[i];
        const int a = residues[i];

        const int inv = inverse_mod_prime(static_cast<int>(current_modulus % m), m);
        const int t = static_cast<int>((static_cast<i64>(normalize_mod(static_cast<i64>(a) - x, m)) * inv) % m);
        x += current_modulus * static_cast<i64>(t);
        current_modulus *= static_cast<i64>(m);
    }

    return normalize_mod(x, kModulus);
}

int solve_mod(const u64 n, const bool allow_multithreading, const unsigned requested_threads) {
    const std::vector<int> moduli(std::begin(kPrimeFactors), std::end(kPrimeFactors));
    std::vector<int> residues(moduli.size(), 0);

    const unsigned thread_count = choose_thread_count(allow_multithreading, requested_threads, moduli.size());
    std::atomic<std::size_t> next_idx(0ULL);

    auto worker = [&]() {
        while (true) {
            const std::size_t idx = next_idx.fetch_add(1ULL, std::memory_order_relaxed);
            if (idx >= moduli.size()) {
                break;
            }
            residues[idx] = solve_mod_prime(moduli[idx], n);
        }
    };

    std::vector<std::thread> pool;
    pool.reserve(thread_count);
    for (unsigned t = 0U; t < thread_count; ++t) {
        pool.emplace_back(worker);
    }
    for (std::thread& th : pool) {
        th.join();
    }

    return crt_combine(residues, moduli);
}

struct ExactAB {
    std::vector<i64> A;
    std::vector<i64> B;
    std::vector<i64> fact;
};

ExactAB compute_exact_ab(const int max_n) {
    ExactAB result;
    result.A.assign(static_cast<std::size_t>(max_n) + 1ULL, 0);
    result.B.assign(static_cast<std::size_t>(max_n) + 1ULL, 0);
    result.fact.assign(static_cast<std::size_t>(max_n) + 1ULL, 1);

    std::vector<std::vector<i64>> binom(
        static_cast<std::size_t>(max_n) + 1ULL,
        std::vector<i64>(static_cast<std::size_t>(max_n) + 1ULL, 0));

    for (int n = 0; n <= max_n; ++n) {
        binom[static_cast<std::size_t>(n)][0] = 1;
        binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(n)] = 1;
        for (int k = 1; k < n; ++k) {
            binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] =
                binom[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k - 1)] +
                binom[static_cast<std::size_t>(n - 1)][static_cast<std::size_t>(k)];
        }
    }

    for (int n = 1; n <= max_n; ++n) {
        result.fact[static_cast<std::size_t>(n)] = result.fact[static_cast<std::size_t>(n - 1)] * static_cast<i64>(n);
    }

    for (int n = 0; n <= max_n; ++n) {
        i64 a_n = result.fact[static_cast<std::size_t>(n)];
        for (int k = 0; k < n; ++k) {
            a_n += binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] * result.A[static_cast<std::size_t>(k)];
        }

        i64 s_n = 0;
        for (int k = 0; k <= n; ++k) {
            s_n += result.fact[static_cast<std::size_t>(n)] / result.fact[static_cast<std::size_t>(k)];
        }

        i64 b_n = -s_n;
        for (int k = 0; k < n; ++k) {
            b_n += binom[static_cast<std::size_t>(n)][static_cast<std::size_t>(k)] * result.B[static_cast<std::size_t>(k)];
        }

        result.A[static_cast<std::size_t>(n)] = a_n;
        result.B[static_cast<std::size_t>(n)] = b_n;
    }

    return result;
}

bool run_checkpoints(const Options& options) {
    const ExactAB exact = compute_exact_ab(12);

    if (exact.A[10] != 328161643LL || exact.B[10] != -652694486LL) {
        std::cerr << "Checkpoint failed: a(10) decomposition mismatch. got A(10)="
                  << exact.A[10] << ", B(10)=" << exact.B[10] << '\n';
        return false;
    }

    for (int n = 0; n <= 12; ++n) {
        const i64 expected_sum = exact.A[static_cast<std::size_t>(n)] + exact.B[static_cast<std::size_t>(n)];
        const int expected_mod = normalize_mod(expected_sum, kModulus);
        const int got_mod = solve_mod(static_cast<u64>(n), false, 1U);
        if (got_mod != expected_mod) {
            std::cerr << "Checkpoint failed at n=" << n << ": expected " << expected_mod
                      << ", got " << got_mod << '\n';
            return false;
        }
    }

    for (const int p : kPrimeFactors) {
        const u64 period = static_cast<u64>(p) * static_cast<u64>(p - 1);
        constexpr u64 kSpotChecks = 80ULL;
        const u64 max_index = static_cast<u64>(p) + period + kSpotChecks;
        const std::vector<int> A = compute_A_mod_prime(p, max_index);

        for (u64 t = 0ULL; t < kSpotChecks; ++t) {
            const u64 lhs_idx = static_cast<u64>(p) + t;
            const u64 rhs_idx = lhs_idx + period;
            if (A[static_cast<std::size_t>(lhs_idx)] != A[static_cast<std::size_t>(rhs_idx)]) {
                std::cerr << "Checkpoint failed: period spot-check mismatch for p=" << p
                          << " at offset " << t << '\n';
                return false;
            }
        }
    }

    const unsigned threads = choose_thread_count(options.allow_multithreading,
                                                 options.requested_threads,
                                                 std::size_t{5});
    std::cout << "Checkpoints passed (threads=" << threads << ").\n";
    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(options)) {
        return 1;
    }

    const int answer = solve_mod(options.n, options.allow_multithreading, options.requested_threads);
    std::cout << "Answer: " << answer << '\n';
    return 0;
}

Python

def normalize_mod(x, mod_val):
    return x % mod_val

def pow_mod(base, exp, mod_val):
    return pow(base, exp, mod_val)

def inverse_mod_prime(a, p):
    return pow(a, p - 2, p)

def factorial_mod_small(n, p):
    f = 1 % p
    for i in range(1, n + 1):
        f = (f * (i % p)) % p
    return f

def compute_A_mod_prime(p, max_n):
    values = [0] * (max_n + 1)
    comb = [0] * (max_n + 1)
    comb[0] = 1
    
    fact = 1 % p
    
    for n in range(max_n + 1):
        if n > 0:
            comb[n] = 1
            for k in range(n - 1, 0, -1):
                comb[k] = (comb[k] + comb[k - 1]) % p
                
            if n < p:
                fact = (fact * n) % p
            else:
                fact = 0
                
        current = fact
        for k in range(n):
            current = (current + comb[k] * values[k]) % p
            
        values[n] = current
        
    return values

def reduced_index_for_prime(n, p):
    if n < p:
        return n
    period = p * (p - 1)
    return p + ((n - p) % period)

def solve_mod_prime(p, n):
    index = reduced_index_for_prime(n, p)
    A = compute_A_mod_prime(p, index)
    
    n_factorial_mod = 0 if n >= p else factorial_mod_small(n, p)
    return (n_factorial_mod - A[index]) % p

def crt_combine(residues, moduli, full_mod):
    x = 0
    current_modulus = 1
    
    for i in range(len(residues)):
        m = moduli[i]
        a = residues[i]
        
        inv = inverse_mod_prime(current_modulus % m, m)
        t = ((a - x) % m * inv) % m
        x = (x + current_modulus * t) % full_mod
        current_modulus *= m
        
    return x % full_mod

def solve():
    n = 1000000000
    modulus = 77777777
    primes = [7, 11, 73, 101, 137]
    
    residues = []
    for p in primes:
        residues.append(solve_mod_prime(p, n))
        
    ans = crt_combine(residues, primes, modulus)
    return str(ans)

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

Java

import java.util.*;

public class Euler330 {
    static final long DEFAULT_N = 1000000000L;
    static final int MODULUS = 77777777;
    static final int[] PRIME_FACTORS = { 7, 11, 73, 101, 137 };

    static int normalizeMod(long x, int mod) {
        long r = x % mod;
        return (int) (r < 0 ? r + mod : r);
    }

    static int powMod(int base, int exp, int mod) {
        long result = 1;
        long b = normalizeMod(base, mod);
        while (exp > 0) {
            if ((exp & 1) != 0) {
                result = (result * b) % mod;
            }
            b = (b * b) % mod;
            exp >>= 1;
        }
        return (int) result;
    }

    static int inverseModPrime(int a, int p) {
        return powMod(a, p - 2, p);
    }

    static int factorialModSmall(long n, int p) {
        int f = 1 % p;
        for (long i = 1; i <= n; i++) {
            f = (int) ((f * (i % p)) % p);
        }
        return f;
    }

    static int[] computeAModPrime(int p, long maxN) {
        int sz = (int) maxN;
        int[] values = new int[sz + 1];
        int[] comb = new int[sz + 1];
        comb[0] = 1;

        int fact = 1 % p;

        for (int n = 0; n <= sz; n++) {
            if (n > 0) {
                comb[n] = 1;
                for (int k = n - 1; k >= 1; k--) {
                    comb[k] += comb[k - 1];
                    if (comb[k] >= p)
                        comb[k] -= p;
                }

                if (n < p) {
                    fact = (int) (((long) fact * n) % p);
                } else {
                    fact = 0;
                }
            }

            long current = fact;
            for (int k = 0; k < n; k++) {
                current += (long) comb[k] * values[k];
                current %= p;
            }

            values[n] = (int) current;
        }

        return values;
    }

    static long reducedIndexForPrime(long n, int p) {
        if (n < p)
            return n;
        long period = (long) p * (p - 1);
        return p + ((n - p) % period);
    }

    static int solveModPrime(int p, long n) {
        long index = reducedIndexForPrime(n, p);
        int[] A = computeAModPrime(p, index);

        int nFactorialMod = (n >= p) ? 0 : factorialModSmall(n, p);
        return normalizeMod((long) nFactorialMod - A[(int) index], p);
    }

    static int crtCombine(int[] residues, int[] moduli) {
        long x = 0;
        long currentModulus = 1;

        for (int i = 0; i < residues.length; i++) {
            int m = moduli[i];
            int a = residues[i];

            int inv = inverseModPrime((int) (currentModulus % m), m);
            int t = (int) ((normalizeMod(a - x, m) * (long) inv) % m);
            x += currentModulus * t;
            currentModulus *= m;
        }

        return normalizeMod(x, MODULUS);
    }

    public static String solve() {
        int[] residues = new int[PRIME_FACTORS.length];
        for (int i = 0; i < PRIME_FACTORS.length; i++) {
            residues[i] = solveModPrime(PRIME_FACTORS[i], DEFAULT_N);
        }
        return String.valueOf(crtCombine(residues, PRIME_FACTORS));
    }

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