Project Euler 320

题目

Factorials divisible by a huge integer

Let \(N(i)\) be the smallest integer \(n\) such that \(n!\) is divisible by \((i!)^{1234567890}\)

Let \(S(u)=\sum N(i)\) for \(10 \le i \le u\).

\(S(1000)=614538266565663\).

Find \(S(1 000 000) \bmod 10^{18}\).

解决方案

对素数 \(p\)，记 \(p\) 在整数 \(x\) 中的指数为 \(\nu_p(x)=\max\{e:p^e\mid x\}.\)

那么\((i!)^M\mid n!\Longleftrightarrow \forall p,\nu_p(n!)\ge M\cdot \nu_p(i!).\)

因此\(\displaystyle{N(i)=\max_{p\le i}\{\operatorname{inv}_p(M\cdot \nu_p(i!))\}}\)，其中\(\displaystyle{\operatorname{inv}_p(K)=\min\{n:\nu_p(n!)\ge K\}.}\)

接下来问题变成两件事：

1. 快速维护所有需要的 \(\nu_p(i!)\)；
2. 快速计算 \(\operatorname{inv}_p(K)\)。

令\(E_p(i)=\nu_p(i!)\)，那么有\(E_p(i)=E_p(i-1)+\nu_p(i).\)若 \(p\nmid i\)，则 \(\nu_p(i)=0\)，所以 \(E_p(i)=E_p(i-1)\)。

另一方面，\(N(i)\) 随 \(i\) 单调不减，因为 \((i!)^M\) 随 \(i\) 增大只会更难整除。

假设我们已经有 \(n=N(i-1)\)，则对所有素数 \(p\) 已满足\(\nu_p(n!)\ge M\cdot E_p(i-1).\)

当从 \(i-1\) 走到 \(i\) 时，只有 那些 \(p\mid i\) 的 \(E_p\) 才会增加；其它 \(p\) 的约束不变，原来的 \(n\) 仍然满足。

因此更新时只需对 \(i\) 的质因数 \(p\) 计算新约束：\(n\leftarrow \max(n,\operatorname{inv}_p(M\cdot E_p(i)))\)。

Legendre 公式可用于计算\(n!\)有多少个质因子\(p\)：\(\displaystyle{\nu_p(n!)=\sum_{t\ge 1}\left\lfloor \frac{n}{p^t}\right\rfloor}\)。

它可用于二分反求 \(n\)。下面给出一个直接构造最小 \(n\) 的方法。

把 \(n\) 写成 \(p\) 进制：\(\displaystyle{n=\sum_{t\ge 0} d_t p^t, 0\le d_t<p}.\)记数位和\(\displaystyle{s_p(n)=\sum_{t\ge 0} d_t.}\)

对任意 \(k\ge1\)，可以得到\(\displaystyle{\left\lfloor\dfrac{n}{p^k}\right\rfloor=\sum_{t\ge k} d_t p^{t-k}}\)，毕竟可以看成这个数在\(p\)进制下右移了\(k\)位。

那么经过一系列变化，得到：

\(\begin{aligned} \nu_p(n!)&=\sum_{k\ge1}\sum_{t\ge k} d_t p^{t-k}\\ &=\sum_{t\ge1} d_t\sum_{k=1}^{t} p^{t-k}\\ &=\sum_{t\ge1} d_t\cdot \dfrac{p^t-1}{p-1}\\ &=\dfrac{(n-d_0)-(s_p(n)-d_0)}{p-1}\\ &=\dfrac{n-s_p(n)}{p-1} \end{aligned}\)

定义权重\(L_t=\dfrac{p^{t+1}-1}{p-1}=1+p+p^2+\cdots+p^t (t\ge 0).\)

如果我们让 \(n\) 的最低位 \(d_0=0\)，并把 \(d_{t+1}\) 记作 \(q_t\)，则\(\displaystyle{n=\sum_{t\ge 0} q_t p^{t+1}}\)且 \(\displaystyle{\nu_p(n!)=\sum_{t\ge 0} q_t L_t.}\)

这说明：想让 \(\nu_p(n!)\ge K\)，等价于要选一些非负整数 \(q_t\) 使得\(\displaystyle{\sum_{t\ge 0} q_t L_t\ge K}\)，同时让\(\displaystyle{n=\sum_{t\ge 0} q_t p^{t+1}}\)尽量小。

由于\(L_{t+1}=pL_t+1\)，这说明\(L_t\) 增长很快，因此可以从大到小贪心地“凑” \(K\)：

• 先找最大的 \(t\) 使得 \(L_t\le K\)；
• 令\(q_t=\left\lfloor\dfrac{K}{L_t}\right\rfloor, K\leftarrow K-q_tL_t\)，同时把 \(n\) 加上 \(q_t p^{t+1}\)；
• 然后对 \(t-1,t-2,\dots,0\) 继续。

可见，由此得到的\(n\)一定是最小的。

代码

from tools import Factorizer

N= 10 ** 6
M = 1234567890
mod = 10 ** 18

def nu(p: int, n: int) -> int:
    s = 0
    while n:
        n //= p
        s += n
    return s


def cal(p: int, k: int) -> int:
    low = []
    pw_p = []
    l = 1
    pw = p
    while l <= k:
        low.append(l)
        pw_p.append(pw)
        l = p * l + 1
        pw *= p
    rem = k
    res = 0
    for l, pw in zip(reversed(low), reversed(pw_p)):
        q, rem = divmod(rem, l)
        res += q * pw
    return res

factorizer = Factorizer(N)
n = 0
ans = 0
for i in range(1, N + 1):
    tmp = [v[0] for v in factorizer.factor_small(i)]
    for p in tmp:
        kp = nu(p, i)
        cand = cal(p, M * kp)
        if cand > n:
            n = cand
    if i >= 10:
        ans = (ans + n) % mod
print(ans)