Project Euler Problem 952

Project Euler Problem 952

Solution

Answer: 794394453

Let

$$M = 10^7!, \qquad p = 10^9+7.$$

We need the multiplicative order

$$R(p,10^7)=\operatorname{ord}_{M}(p).$$

Because $M=\prod q^{a_q}$ over primes $q\le 10^7$,

$$\operatorname{ord}{M}(p) = \operatorname{lcm}{q\le 10^7} \operatorname{ord}_{q^{a_q}}(p),$$

where

$$a_q=v_q(10^7!).$$

For odd primes $q\neq p$, if

$$d_q=\operatorname{ord}_q(p),$$

then the standard lifting theorem gives

$$\operatorname{ord}_{q^{a_q}}(p) = d_q, q^{\max(0,a_q-s_q)},$$

where

$$s_q=v_q(p^{d_q}-1).$$

Almost always $s_q=1$; the only exceptions up to $10^7$ are

$$q\in{3,5,43,109},$$

for which $s_q=2$.

For $q=2$, since $p=10^9+7\equiv 7\pmod 8$,

$$\operatorname{ord}_{2^a}(p)=2^{a-3}\qquad (a\ge 3).$$

So the task reduces to:

1. sieve primes up to $10^7$,
2. compute $d_q\mid (q-1)$,
3. assemble the LCM prime-exponent-wise,
4. reduce modulo $10^9+7$.

A clean Python implementation is:

from array import array

N = 10_000_000
MOD = 1_000_000_007
P = MOD

# ---------- smallest prime factor sieve ----------
spf = array('I', [0]) * (N + 1)
primes = array('I')

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 p > spf[i] or v > N:
            break
        spf[v] = p

# ---------- build lcm prime exponents ----------
maxexp = {}

# q = 2
a2 = 0
t = N
while t:
    t //= 2
    a2 += t

maxexp[2] = a2 - 3

# odd primes
for q in primes[1:]:

    # exponent of q in N!
    a = 0
    t = N
    while t:
        t //= q
        a += t

    # d = ord_q(P), starting from q-1
    d = q - 1
    x = q - 1
    last = 0

    while x > 1:
        r = spf[x]
        if r != last:
            while d % r == 0 and pow(P, d // r, q) == 1:
                d //= r
            last = r
        x //= r

    # factor d into global lcm
    x = d
    while x > 1:
        r = spf[x]
        e = 0
        while x % r == 0:
            x //= r
            e += 1
        if e > maxexp.get(r, 0):
            maxexp[r] = e

    # lifting exponent s_q
    s = 1
    if pow(P, d, q * q) == 1:
        s = 2

    e = max(0, a - s)

    if e > maxexp.get(q, 0):
        maxexp[q] = e

# ---------- reconstruct answer mod MOD ----------
ans = 1
for r, e in maxexp.items():
    ans = (ans * pow(r, e, MOD)) % MOD

print(ans)

This evaluates to:

Answer: 794394453