Project Euler 288

题目

An enormous factorial

For any prime \(p\) the number \(N(p,q)\) is defined by \(N(p,q) = \sum_{n=0\text{ to }q} T_n*p^n\) with \(T_n\) generated by the following random number generator:

\[\begin{aligned} S_0 &= 290797\\ S_{n+1} &= S_n^2 \bmod 50515093\\ T_n &= S_n \bmod p \end{aligned}\]

Let \(\text{Nfac}(p,q)\) be the factorial of \(N(p,q)\).

Let \(\text{NF}(p,q)\) be the number of factors \(p\) in \(\text{Nfac}(p,q)\).

You are given that \(\text{NF}(3,10000) \bmod 3^{20}=624955285\).

Find \(\text{NF}(61,10^7) \bmod 61^{10}\).

解决方案

设\(f(n, p)\)是质因子\(p\)在\(n!\)中的次数，那么\(f(n,p)=\left\lfloor\dfrac{n}{p}\right\rfloor+\left\lfloor\dfrac{n}{p^2}\right\rfloor+\left\lfloor\dfrac{n}{p^3}\right\rfloor+\dots\)，每一项分别表示\(1\sim n\)中有多少个数是\(p,p^2,p^3\dots\)的倍数。

那么，\(\text{NF}(p,q)=f(N(p,q),p)\)。

不过，恰巧的是，\(N(p,q)\)是一个关于\(p\)的\(q\)项式。因此，以\(f(n,p)\)的前两项为例子，有：

\(\begin{aligned} &\left\lfloor\dfrac{N(p,q)}{p}\right\rfloor=\sum_{i=1}^{q} T_ip^{i-1}\\ &\left\lfloor\dfrac{N(p,q)}{p^2}\right\rfloor=\sum_{i=2}^{q} T_ip^{i-2}\\ &\dots\\ &\left\lfloor\dfrac{N(p,q)}{p^q}\right\rfloor=T_q\\\\ \end{aligned}\)

那么，可以写出

\[\text{NF}(p,q)=\sum_{k=1}^q\sum_{i=k}^qT_ip^{i-k}=\sum_{i=1}^qT_i\sum_{k=0}^{i-1}p^k\]

那么，式子两边可以分别同步计算\(T_i\)和\(\sum_{k=0}^{i-1}p^k\)的值，最终重新相乘即可。

一个可以优化的地方：由于取模数恰好是一个\(p^K\)，其中\(K=10\)，因此计算\(\sum_{k=0}^{i-1}p^i\% p^K\)时，会发现到了第\(K\)项之后，值都是一样的。最终，计算过程节省了取模这一步骤。

代码

P = 61
Q = 10 ** 7
mod = 61 ** 10
s = 290797
sp = 1
ans = 0
for i in range(Q):
    s = s * s % 50515093
    t = s % P
    ans = (ans + t * sp) % mod
    sp = (sp * P + 1) % mod
print(ans)

P = 61
Q = 10 ** 7
N = 10
mod = P ** N
s = 290797
sp = 1
ans = 0
for i in range(Q):
    s = s * s % 50515093
    t = s % P
    ans = (ans + t * sp) % mod
    if i < N:
        sp = sp * P + 1
print(ans)