Project Euler 854

题目

Pisano Periods 2

For every positive integer \(n\) the Fibonacci sequence modulo \(n\) is periodic. The period depends on the value of \(n\). This period is called the Pisano period for \(n\), often shortened to \(\pi(n)\).

Define \(M(p)\) as the largest integer \(n\) such that \(\pi(n) = p\), and define \(M(p) = 1\) if there is no such \(n\).

For example, there are three values of \(n\) for which \(\pi(n)\) equals \(18\): \(19, 38, 76\). Therefore \(M(18) = 76\).

Let the product function \(P(n)\) be:

\[P(n)=\prod_{p = 1}^{n}M(p).\]

You are given: \(P(10)=264\).

Find \(P(1\,000\,000)\bmod 1\,234\,567\,891\).

解决方案

记\(F_0=0, F_1=1, F_{n+2}=F_{n+1}+F_n(n\ge 0),\)其中 \(F_n\) 表示第 \(n\) 项斐波那契数。\(L_0=2, L_1=1, L_{n+2}=L_{n+1}+L_n(n\ge 0),\)其中 \(L_n\) 表示第 \(n\) 项 Lucas 数。

斐波那契递推是二阶线性递推，因此模 \(m\) 下的整个序列由有序对 \((F_n,F_{n+1})\bmod m\) 决定。特别地，\(\pi(m)=\min\{k\ge1:(F_k,F_{k+1})\equiv(0,1)\pmod m\}.\)于是对任意 \(k\ge1\)，

\[ \pi(m)\mid k\iff (F_k,F_{k+1})\equiv(0,1)\pmod m \iff \begin{cases} F_k\equiv 0\pmod m,\\ F_{k+1}\equiv 1\pmod m, \end{cases} \]

等价于\(m\mid F_k\land m\mid(F_{k+1}-1)\iff m\mid \gcd(F_k,F_{k+1}-1).\)因此令\(G(k)=\gcd(F_k,F_{k+1}-1),\) 则“所有满足 \(\pi(m)\mid k\) 的模数 \(m\)”恰好是 \(G(k)\) 的所有正因子；特别地，满足 \(\pi(m)\mid k\) 的最大模数就是 \(G(k)\)。

论文的引理6.8给出，当 \(m>2\) 时，斐波那契与 Lucas 的 Pisano 周期都是偶数。并且论文给出了更精确的\(\pi\)的值域描述：\(\{3\}\cup\{n:n\in 2\mathbb Z,n\ge 6\}.\)

于是，\(\pi(1)=1\)，所以 \(M(1)=1\)；\(\pi(2)=3\)，所以 \(M(3)=2\)；对任意奇数 \(p>3\)，根本不存在 \(n\) 使得 \(\pi(n)=p\)，故 \(M(p)=1\)。

因此在积 \(P(N)\) 中，所有奇数 \(p\) 只有 \(p=3\) 会贡献一个因子 \(2\)，其余全为 \(1\)。

论文的定理6.9同样证明了给定周期 \(k\) 时产生该周期的最大模数。结论是：若 \(k\ge 6\) 为偶数，则产生 Pisano 周期 \(k\) 的最大模数为

\[ m_F= \begin{cases} F_{k/2}, & k\equiv 0\pmod 4,\\ L_{k/2}, & k\equiv 2\pmod 4, \end{cases} \]

综合奇偶两部分，可得：

\[ M(p)= \begin{cases} 2,& p=3,\\ 1,& p\equiv 1\pmod2\land p\ne3,\\ F_{p/2},& p\equiv 0\pmod4,\\ L_{p/2},& p\equiv 2\pmod4. \end{cases} \] 其中 \(p=2,4\) 的两条偶数公式也仍然给出 \(1\)（因为 \(F_1=L_1=1,F_2=1\)），与“这些周期实际上并不存在”并不冲突（此时 \(M(p)\) 规定为 \(1\)）。

因此最终有

\[ P(N)=2\cdot\prod_{t=1}^{\lfloor N/2\rfloor}g(t),g(t)=\begin{cases} F_t,& t\equiv 0\pmod 2,\\ L_t,& t\equiv 1\pmod 2, \end{cases} \]

其中前面的 \(2\) 来自 \(M(3)=2\)。

代码

N = 1_000_000
MOD = 1234567891

def solve(N):
    prod = 2 if N >= 3 else 1
    half = N // 2

    f_prev, f = 0, 1
    l_prev, l = 2, 1

    for t in range(1, half + 1):
        m = f if (t & 1) == 0 else l
        prod = (prod * m) % MOD
        f_prev, f = f, (f_prev + f) % MOD
        l_prev, l = l, (l_prev + l) % MOD

    return prod

print(solve(N))