Problem 952: Order Modulo Factorial

Project Euler Problem 952 Solution

EulerSolve provides an optimized solution for Project Euler Problem 952, Order Modulo Factorial, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The quantity in this problem is the multiplicative order of $$a=10^9+7$$ modulo $$n!=10{,}000{,}000!.$$ In other words, we want the smallest positive integer $R$ such that $$a^R\equiv 1 \pmod{n!},$$ and then we report $R \bmod a$. Because $a$ is much larger than $n$, it does not divide $n!$, so the order is well defined. The brute-force viewpoint is hopeless. The number $n!$ is enormous, and the order itself is far beyond direct enumeration. The solution works only because the modulus $n!$ has a very structured prime-power factorization, and multiplicative orders behave cleanly on prime powers. Mathematical Approach The key is to replace one huge congruence modulo $n!$ by many local congruences modulo prime powers, solve each local problem, and then recombine them through an LCM. Decomposing $n!$ into prime powers For each prime $q\le n$, let $$e_q=v_q(n!)=\sum_{j\ge 1}\left\lfloor\frac{n}{q^j}\right\rfloor.$$ This is Legendre's formula, and it gives the exact exponent of $q$ in the factorization of $n!$. Therefore $$n!=\prod_{q\le n} q^{e_q},$$ where the product runs over primes....

Detailed mathematical approach

Problem Summary

The quantity in this problem is the multiplicative order of

$$a=10^9+7$$

modulo

$$n!=10{,}000{,}000!.$$

In other words, we want the smallest positive integer $R$ such that

$$a^R\equiv 1 \pmod{n!},$$

and then we report $R \bmod a$. Because $a$ is much larger than $n$, it does not divide $n!$, so the order is well defined.

The brute-force viewpoint is hopeless. The number $n!$ is enormous, and the order itself is far beyond direct enumeration. The solution works only because the modulus $n!$ has a very structured prime-power factorization, and multiplicative orders behave cleanly on prime powers.

Mathematical Approach

The key is to replace one huge congruence modulo $n!$ by many local congruences modulo prime powers, solve each local problem, and then recombine them through an LCM.

Decomposing $n!$ into prime powers

For each prime $q\le n$, let

$$e_q=v_q(n!)=\sum_{j\ge 1}\left\lfloor\frac{n}{q^j}\right\rfloor.$$

This is Legendre's formula, and it gives the exact exponent of $q$ in the factorization of $n!$. Therefore

$$n!=\prod_{q\le n} q^{e_q},$$

where the product runs over primes.

Since the prime-power factors are pairwise coprime, the Chinese remainder theorem implies that the order modulo the full factorial is the least common multiple of the local orders:

$$\operatorname{ord}_{n!}(a)=\operatorname{lcm}_{q\le n} \operatorname{ord}_{q^{e_q}}(a).$$

So the problem is reduced to computing $\operatorname{ord}_{q^{e_q}}(a)$ for every prime $q\le n$.

The odd-prime local order and its lifting depth

Assume first that $q$ is odd. Let

$$t_q=\operatorname{ord}_q(a).$$

By Fermat's little theorem, $t_q\mid (q-1)$. That is why the implementations start from $q-1$, factor it, and repeatedly test whether a prime factor can be removed while keeping the congruence $a^{t_q}\equiv 1 \pmod q$. After all removable factors are stripped away, the remaining value is the exact order modulo $q$.

The next question is how this order changes when the modulus is lifted from $q$ to $q^2,q^3,\dots,q^{e_q}$. Define

$$s_q=v_q(a^{t_q}-1).$$

This is the first level at which the congruence $a^{t_q}\equiv 1$ stops being automatic. For odd $q$, the lifting-the-exponent principle gives

$$v_q\!\left(a^{t_q q^m}-1\right)=s_q+m \qquad (m\ge 0).$$

Hence the minimal exponent that reaches modulus $q^{e_q}$ is

$$\operatorname{ord}_{q^{e_q}}(a)=t_q\,q^{\max(0,e_q-s_q)}.$$

This formula is exactly what the implementations use. They compute $t_q$, then test whether $a^{t_q}\equiv 1$ modulo $q^2,q^3,\dots$ until the lift fails, and the last successful level is $s_q$.

The special $2$-adic branch

The prime $q=2$ is different because the unit group modulo $2^e$ does not behave like the odd-prime case. The implementations therefore use explicit $2$-adic formulas.

$$u=v_2(a-1)$$

and $a\equiv 1 \pmod 4$, then $a$ is already very close to 1 in the $2$-adic sense, and

$$\operatorname{ord}_{2^e}(a)= \begin{cases} 1, & e\le u,\\ 2^{e-u}, & e>u. \end{cases}$$

If instead

$$w=v_2(a+1)$$

and $a\equiv 3 \pmod 4$, then the order must be even once $e\ge 2$, and the correct formula is

$$\operatorname{ord}_{2^e}(a)= \begin{cases} 1, & e=1,\\ 2^{\max(1,e-w)}, & e\ge 2. \end{cases}$$

For this problem, $a=10^9+7\equiv 3 \pmod 4$ and $v_2(a+1)=3$, so the contribution from the $2$-power inside $n!$ is governed by $2^{\max(1,e_2-3)}$.

Reassembling the global order from prime exponents

Once every local order $\operatorname{ord}_{q^{e_q}}(a)$ is known, taking their LCM is still too large to do with ordinary integers. The clean way is to collect prime exponents.

For each prime $r\le n$, define

$$E_r=\max_{q\le n} v_r\!\left(\operatorname{ord}_{q^{e_q}}(a)\right).$$

Then

$$\operatorname{ord}_{n!}(a)=\prod_{r\le n} r^{E_r}.$$

This is why the implementations maintain only a table of maximal exponents. For odd $q$, the factors coming from $t_q$ are among the prime divisors of $q-1$, and the lifting step contributes an extra power of $q$ itself. No other primes can appear.

Worked Example: $n=12$

The smaller case $n=12$ shows the whole mechanism in a concrete way. First,

$$12!=2^{10}\cdot 3^5\cdot 5^2\cdot 7\cdot 11.$$

So we need the orders modulo $2^{10}$, $3^5$, $5^2$, $7$, and $11$.

For $q=2$, we have $a\equiv 3 \pmod 4$ and $v_2(a+1)=3$, so

$$\operatorname{ord}_{2^{10}}(a)=2^{\max(1,10-3)}=2^7=128.$$

For $q=3$, $a\equiv -1 \pmod 3$, hence $t_3=2$. Also $a\equiv -1 \pmod 9$ but not modulo $27$, so $s_3=2$. Therefore

$$\operatorname{ord}_{3^5}(a)=2\cdot 3^{5-2}=54.$$

For $q=5$, $a\equiv 2 \pmod 5$, whose order modulo 5 is $t_5=4$. Since $2^4-1=15$, we get $s_5=1$, and therefore

$$\operatorname{ord}_{5^2}(a)=4\cdot 5=20.$$

The orders modulo 7 and 11 divide 6 and 10 respectively, so they introduce no new prime powers beyond what is already present in 128, 54, and 20. Hence

$$\operatorname{ord}_{12!}(a)=\operatorname{lcm}(128,54,20)=2^7\cdot 3^3\cdot 5=17280.$$

This is exactly the kind of local-to-global reconstruction used for the full $n=10^7$ case.

How the Code Works

Prime preprocessing and factorial valuations

The C++, Python, and Java implementations first build a smallest-prime-factor sieve up to $n$. That single preprocessing step serves two purposes: it produces the full prime list $q\le n$, and it allows fast factorizations of numbers like $q-1$, which are needed when reducing candidate orders.

For each prime $q$, the implementation evaluates Legendre's formula to obtain $e_q=v_q(n!)$. This identifies the exact prime-power modulus $q^{e_q}$ that must be handled.

Computing each local order

When $q$ is odd, the implementation starts from the candidate $q-1$, factors it, and removes one prime factor at a time whenever modular exponentiation shows that the smaller exponent still works modulo $q$. That produces $t_q=\operatorname{ord}_q(a)$.

Next it determines the lift depth $s_q$ by testing the same exponent $t_q$ modulo $q^2,q^3,\dots$. If the congruence survives to a higher power, the order does not grow yet; once it fails, every additional power of $q$ in the modulus multiplies the order by $q$. The implementation records both parts of the factorization: the prime factors of $t_q$, and the extra $q$-power coming from $e_q-s_q$.

For $q=2$, the implementation bypasses the odd-prime logic and applies the closed $2$-adic formulas directly using $v_2(a-1)$ or $v_2(a+1)$, depending on whether $a\equiv 1$ or $3 \pmod 4$.

Accumulating the LCM and reducing modulo $a$

Instead of forming gigantic local orders and then taking a literal least common multiple, the implementation updates a global table of maximum prime exponents. After all primes $q\le n$ are processed, the answer is reconstructed as

$$\prod_{r\le n} r^{E_r}\pmod a.$$

The C++ version uses 128-bit arithmetic for ordinary modular products and switches to multiprecision arithmetic when larger prime powers must be tested exactly. Python gets arbitrary precision automatically, and the Java version uses built-in big integers for the deeper lifting checks.

Complexity Analysis

The sieve up to $n=10^7$ is $O(n)$ time and $O(n)$ memory. This is the dominant storage cost, since both the smallest-prime-factor table and the global exponent table are indexed up to $n$.

For each prime $q\le n$, evaluating $e_q$ costs $O(\log_q n)$. Factoring $q-1$ is fast because the sieve is already available, and reducing $q-1$ to the true order modulo $q$ requires only a small number of modular exponentiation tests, one for each prime-factor occurrence that might be removed. The lifting phase adds a few more modular exponentiations until the congruence fails at the first unattainable prime power.

So the implementation is best viewed as linear preprocessing plus a per-prime amount of arithmetic that is polylogarithmic in the modulus size. The crucial point is that no part of the algorithm ever tries to build $n!$ itself or search for the order by repeated multiplication.

Footnotes and References

Problem page: https://projecteuler.net/problem=952
Multiplicative order: Wikipedia - Multiplicative order
Legendre's formula: Wikipedia - Legendre's formula
Chinese remainder theorem: Wikipedia - Chinese remainder theorem
Lifting-the-exponent lemma: Wikipedia - Lifting-the-exponent lemma
Modular exponentiation: Wikipedia - Modular exponentiation

Problem 952 source code

C++

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

#include <boost/multiprecision/cpp_int.hpp>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
using boost::multiprecision::cpp_int;

u64 mod_pow_u64(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    base %= mod;
    while (exp > 0) {
        if (exp & 1ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * base) % mod);
        }
        base = static_cast<u64>((static_cast<u128>(base) * base) % mod);
        exp >>= 1ULL;
    }
    return result;
}

int v_p_factorial(int n, int p) {
    int e = 0;
    while (n > 0) {
        n /= p;
        e += n;
    }
    return e;
}

int v2_u64(u64 x) {
    int v = 0;
    while ((x & 1ULL) == 0ULL) {
        x >>= 1ULL;
        ++v;
    }
    return v;
}

std::pair<std::vector<int>, std::vector<int>> build_spf_and_primes(int n) {
    std::vector<int> spf(n + 1, 0);
    std::vector<int> primes;
    primes.reserve(n / 10);

    for (int i = 2; i <= n; ++i) {
        if (spf[i] == 0) {
            spf[i] = i;
            primes.push_back(i);
        }
        for (int p : primes) {
            const i64 v = static_cast<i64>(p) * i;
            if (v > n || p > spf[i]) {
                break;
            }
            spf[static_cast<int>(v)] = p;
        }
    }

    return {std::move(spf), std::move(primes)};
}

std::vector<std::pair<int, int>> factorize_int(int x, const std::vector<int>& spf) {
    std::vector<std::pair<int, int>> fac;
    while (x > 1) {
        int p = spf[x];
        int e = 0;
        do {
            x /= p;
            ++e;
        } while (x > 1 && spf[x] == p);
        fac.push_back({p, e});
    }
    return fac;
}

u64 compute_R_mod(u64 p, int n, u64 out_mod) {
    const auto [spf, primes] = build_spf_and_primes(n);
    std::vector<int> max_exp(n + 1, 0);

    for (int q : primes) {
        if (q > n) {
            break;
        }

        const int e = v_p_factorial(n, q);

        if (q == 2) {
            int exp2 = 0;
            if (e >= 2) {
                if ((p & 3ULL) == 1ULL) {
                    const int u = v2_u64(p - 1ULL);
                    exp2 = (e <= u) ? 0 : (e - u);
                } else {
                    const int v = v2_u64(p + 1ULL);
                    exp2 = (e <= v) ? 1 : (e - v);
                }
            }
            if (exp2 > max_exp[2]) {
                max_exp[2] = exp2;
            }
            continue;
        }

        u64 t = static_cast<u64>(q - 1);
        const auto fac_qm1 = factorize_int(q - 1, spf);
        const u64 base_q = p % static_cast<u64>(q);

        for (const auto& [r, cnt] : fac_qm1) {
            for (int i = 0; i < cnt; ++i) {
                if (t % static_cast<u64>(r) == 0ULL &&
                    mod_pow_u64(base_q, t / static_cast<u64>(r), static_cast<u64>(q)) == 1ULL) {
                    t /= static_cast<u64>(r);
                } else {
                    break;
                }
            }
        }

        u64 tmp = t;
        for (const auto& [r, _] : fac_qm1) {
            int cnt = 0;
            while (tmp % static_cast<u64>(r) == 0ULL) {
                tmp /= static_cast<u64>(r);
                ++cnt;
            }
            if (cnt > max_exp[r]) {
                max_exp[r] = cnt;
            }
        }
        assert(tmp == 1ULL);

        int s = 1;
        if (e > 1) {
            const u64 mod_q2 = static_cast<u64>(q) * static_cast<u64>(q);
            if (mod_pow_u64(p % mod_q2, t, mod_q2) == 1ULL) {
                s = 2;
                if (e > 2) {
                    cpp_int base = p;
                    cpp_int exp = t;
                    cpp_int mod = cpp_int(q) * q * q;
                    for (int k = 3; k <= e; ++k) {
                        if (boost::multiprecision::powm(base, exp, mod) == 1) {
                            s = k;
                            mod *= q;
                        } else {
                            break;
                        }
                    }
                }
            }
        }

        const int extra = e - s;
        if (extra > max_exp[q]) {
            max_exp[q] = extra;
        }
    }

    u64 ans = 1 % out_mod;
    for (int r = 2; r <= n; ++r) {
        if (max_exp[r] > 0) {
            ans = static_cast<u64>((static_cast<u128>(ans) * mod_pow_u64(static_cast<u64>(r), static_cast<u64>(max_exp[r]), out_mod)) % out_mod);
        }
    }

    return ans;
}

void run_validations() {
    assert(compute_R_mod(7ULL, 4, 1'000'000'007ULL) == 2ULL);
    assert(compute_R_mod(1'000'000'007ULL, 12, 1'000'000'007ULL) == 17'280ULL);
}

}  // namespace

int main() {
    run_validations();
    constexpr u64 kP = 1'000'000'007ULL;
    constexpr int kN = 10'000'000;
    std::cout << compute_R_mod(kP, kN, kP) << '\n';
    return 0;
}

Python

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def v_p_factorial(n, p):
    e = 0
    while n > 0:
        n //= p
        e += n
    return e

def v2_int(x):
    v = 0
    while (x & 1) == 0:
        x >>= 1
        v += 1
    return v

def build_spf_and_primes(n):
    spf = [0] * (n + 1)
    primes = []
    
    for i in range(2, n + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            v = p * i
            if v > n or p > spf[i]:
                break
            spf[v] = p
            
    return spf, primes

def factorize_int(x, spf):
    fac = []
    while x > 1:
        p = spf[x]
        e = 0
        while x > 1 and spf[x] == p:
            x //= p
            e += 1
        fac.append((p, e))
    return fac

def compute_R_mod(p, n, out_mod):
    spf, primes = build_spf_and_primes(n)
    max_exp = [0] * (n + 1)
    
    for q in primes:
        if q > n:
            break
            
        e = v_p_factorial(n, q)
        
        if q == 2:
            exp2 = 0
            if e >= 2:
                if (p & 3) == 1:
                    u = v2_int(p - 1)
                    exp2 = 0 if e <= u else (e - u)
                else:
                    v = v2_int(p + 1)
                    exp2 = 1 if e <= v else (e - v)
            if exp2 > max_exp[2]:
                max_exp[2] = exp2
            continue
            
        t = q - 1
        fac_qm1 = factorize_int(q - 1, spf)
        base_q = p % q
        
        for r, cnt in fac_qm1:
            for _ in range(cnt):
                if t % r == 0 and mod_pow(base_q, t // r, q) == 1:
                    t //= r
                else:
                    break
                    
        tmp = t
        for r, _ in fac_qm1:
            cnt = 0
            while tmp % r == 0:
                tmp //= r
                cnt += 1
            if cnt > max_exp[r]:
                max_exp[r] = cnt
                
        s = 1
        if e > 1:
            mod_q2 = q * q
            if mod_pow(p % mod_q2, t, mod_q2) == 1:
                s = 2
                if e > 2:
                    base = p
                    exp = t
                    mod = q * q * q
                    for k in range(3, e + 1):
                        if mod_pow(base, exp, mod) == 1:
                            s = k
                            mod *= q
                        else:
                            break
                            
        extra = e - s
        if extra > max_exp[q]:
            max_exp[q] = extra
            
    ans = 1 % out_mod
    for r in range(2, n + 1):
        if max_exp[r] > 0:
            ans = (ans * mod_pow(r, max_exp[r], out_mod)) % out_mod
            
    return ans

def solve():
    return str(compute_R_mod(1000000007, 10000000, 1000000007))

if __name__ == "__main__":
    assert compute_R_mod(7, 4, 1000000007) == 2
    assert compute_R_mod(1000000007, 12, 1000000007) == 17280
    print(solve())

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;

public class Euler952 {

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

    static int vPFactorial(int n, int p) {
        int e = 0;
        while (n > 0) {
            n /= p;
            e += n;
        }
        return e;
    }

    static int v2Int(long x) {
        int v = 0;
        while ((x & 1) == 0) {
            x >>= 1;
            ++v;
        }
        return v;
    }

    static class SpfResult {
        int[] spf;
        List<Integer> primes;

        SpfResult(int[] spf, List<Integer> primes) {
            this.spf = spf;
            this.primes = primes;
        }
    }

    static SpfResult buildSpfAndPrimes(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 10);

        for (int i = 2; i <= n; ++i) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > n || p > spf[i]) {
                    break;
                }
                spf[(int) v] = p;
            }
        }
        return new SpfResult(spf, primes);
    }

    static class Factor {
        int p;
        int e;

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

    static List<Factor> factorizeInt(int x, int[] spf) {
        List<Factor> fac = new ArrayList<>();
        while (x > 1) {
            int p = spf[x];
            int e = 0;
            do {
                x /= p;
                ++e;
            } while (x > 1 && spf[x] == p);
            fac.add(new Factor(p, e));
        }
        return fac;
    }

    static long computeRMod(long p, int n, long outMod) {
        SpfResult res = buildSpfAndPrimes(n);
        int[] spf = res.spf;
        List<Integer> primes = res.primes;
        int[] maxExp = new int[n + 1];

        for (int q : primes) {
            if (q > n) {
                break;
            }

            int e = vPFactorial(n, q);

            if (q == 2) {
                int exp2 = 0;
                if (e >= 2) {
                    if ((p & 3) == 1) {
                        int u = v2Int(p - 1);
                        exp2 = (e <= u) ? 0 : (e - u);
                    } else {
                        int v = v2Int(p + 1);
                        exp2 = (e <= v) ? 1 : (e - v);
                    }
                }
                if (exp2 > maxExp[2]) {
                    maxExp[2] = exp2;
                }
                continue;
            }

            long t = q - 1;
            List<Factor> facQm1 = factorizeInt(q - 1, spf);
            long baseQ = p % q;

            for (Factor f : facQm1) {
                for (int i = 0; i < f.e; ++i) {
                    if (t % f.p == 0 && modPow(baseQ, t / f.p, q) == 1) {
                        t /= f.p;
                    } else {
                        break;
                    }
                }
            }

            long tmp = t;
            for (Factor f : facQm1) {
                int cnt = 0;
                while (tmp % f.p == 0) {
                    tmp /= f.p;
                    ++cnt;
                }
                if (cnt > maxExp[f.p]) {
                    maxExp[f.p] = cnt;
                }
            }

            int s = 1;
            if (e > 1) {
                long modQ2 = (long) q * q;
                if (modPow(p % modQ2, t, modQ2) == 1) {
                    s = 2;
                    if (e > 2) {
                        BigInteger base = BigInteger.valueOf(p);
                        BigInteger exp = BigInteger.valueOf(t);
                        BigInteger mod = BigInteger.valueOf(q).pow(3);
                        BigInteger bigQ = BigInteger.valueOf(q);
                        for (int k = 3; k <= e; ++k) {
                            if (base.modPow(exp, mod).equals(BigInteger.ONE)) {
                                s = k;
                                mod = mod.multiply(bigQ);
                            } else {
                                break;
                            }
                        }
                    }
                }
            }

            int extra = e - s;
            if (extra > maxExp[q]) {
                maxExp[q] = extra;
            }
        }

        long ans = 1 % outMod;
        for (int r = 2; r <= n; ++r) {
            if (maxExp[r] > 0) {
                ans = (ans * modPow(r, maxExp[r], outMod)) % outMod;
            }
        }

        return ans;
    }

    public static String solve() {
        return Long.toString(computeRMod(1000000007L, 10000000, 1000000007L));
    }

    public static void main(String[] args) {
        if (computeRMod(7L, 4, 1000000007L) != 2L || computeRMod(1000000007L, 12, 1000000007L) != 17280L) {
            System.out.println("Validation failed");
            return;
        }
        System.out.println(solve());
    }
}