Project Euler 853

题目

Pisano Periods 1

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)\).

There are three values of \(n\) for which \(\pi(n)\) equals \(18\): \(19\), \(38\) and \(76\). The sum of those smaller than \(50\) is \(57\).

Find the sum of the values of \(n\) smaller than \(1\,000\,000\,000\) for which \(\pi(n)\) equals \(120\).

解决方案

定义矩阵\(\displaystyle{A=\begin{pmatrix}1&1\\1&0\end{pmatrix}.}\)

引理 1：对所有整数 \(k\ge 1,A^k=\begin{pmatrix}F_{k+1} & F_k\\F_k & F_{k-1}\end{pmatrix}.\)

证明： \(k=1\) 时显然成立。假设对 \(k\) 成立，则 \(A^{k+1}=A^kA=\begin{pmatrix}F_{k+1} & F_k\\F_k & F_{k-1}\end{pmatrix}\begin{pmatrix}1&1\\1&0\end{pmatrix}=\begin{pmatrix}F_{k+1}+F_k & F_{k+1}\\F_k+F_{k-1} & F_k\end{pmatrix}=\begin{pmatrix}F_{k+2} & F_{k+1}\\F_{k+1} & F_k\end{pmatrix},\) 完成归纳。由此得到\(({F_{k+1}},{F_k})^T=A^k(1,0)^T.\)

记模 \(n\) 的状态对为 \((F_k,F_{k+1})\bmod n\)。由于递推是确定的，这个状态对回到 \((0,1)\) 的最小正步数正是 \(\pi(n)\)。

引理 2：对 \(t\ge 1\)，\((F_t,F_{t+1})\equiv(0,1)\pmod n\Longleftrightarrow A^t\equiv I\pmod n.\)

证明： 由引理 1，\(A^t\equiv I\pmod n\Longleftrightarrow\begin{pmatrix}F_{t+1} & F_t\\F_t & F_{t-1}\end{pmatrix}\equiv\begin{pmatrix}1&0\\0&1\end{pmatrix}\pmod n\Longleftrightarrow F_t\equiv 0,F_{t+1}\equiv 1\pmod n.\)而 \(F_t\equiv 0,F_{t+1}\equiv 1\) 正是状态对回到初态 \((0,1)\) 的条件。

因此\(\pi(n)=M\)当且仅当\(A^M\equiv I\pmod n\)且\(\forall d\mid M,0<d<M,A^d\not\equiv I\pmod n\)均成立。

由由引理 1，\(A^M\equiv I\pmod n\) 等价于\(F_M\equiv 0\pmod n, F_{M+1}\equiv 1\pmod n.\)利用 \(F_{M+1}=F_M+F_{M-1}\)，当 \(F_M\equiv 0\) 时有\(F_{M+1}\equiv 1 \iff F_{M-1}\equiv 1\pmod n.\)所以\(A^M\equiv I\pmod n\Longleftrightarrow n\mid F_M\land n\mid (F_{M-1}-1).\)

令\(G=\gcd(F_M,F_{M-1}-1),\)则必要且充分条件是\(A^M\equiv I\pmod n \Longleftrightarrow n\mid G.\)

这一步极其关键：所有满足 \(\pi(n)\mid M\) 的 \(n\) 必然属于 \(G\) 的因子集合。于是我们不必在 \([1,N)\) 上枚举 \(n\)，只需枚举 \(G\) 的所有因子即可。

现在要从 \(n\mid G\) 中筛出那些使矩阵阶恰好等于 \(M\)的 \(n\)。

设 \(t\) 是 \(A\) 在模 \(n\) 下的阶（即最小正 \(t\) 使 \(A^t\equiv I\pmod n\)）。若已知 \(A^M\equiv I\pmod n\)，则 \(t\mid M\)。我们想要 \(t=M\)。

引理3：设 \(t\mid M\)。则 \(t<M\) 当且仅当存在素数 \(p\mid M\) 使\(A^{M/p}\equiv I\pmod n.\)

证明： 若存在 \(p\mid M\) 且 \(A^{M/p}\equiv I\)，则 \(t\mid M/p<M\)，因此 \(t<M\)。反之若 \(t<M\)，则 \(q=M/t>1\)。取 \(q\) 的任意素因子 \(p\)，则 \(p\mid M\) 且 \(t\mid M/p\)，从而\(A^{M/p}=(A^t)^{(M/p)/t}\equiv I.\)

因此判定 \(\pi(n)=M\) 的策略为：先限制 \(n\mid G\)；对每个素因子 \(p\mid M\)，都要求\(A^{M/p}\not\equiv I\pmod n.\)

又由引理 2，\(A^{k}\equiv I\pmod n\) 等价于\(F_k\equiv 0\pmod n, F_{k+1}\equiv 1\pmod n.\)

最终，我们先求 \(F_M,F_{M-1},G=\gcd(F_M,F_{M-1}-1)\)；随后枚举 \(G\) 的所有因子；并取 \(M\) 的不同素因子集合；对每个因子 \(n\)，检查所有 \(M/p\) 都不能回到 \((0,1)\)，通过则累加。

代码

from tools import factorization, divisors, gcd

M = 120
N = 10 ** 9

def fib_pair_big(n):
    if n == 0:
        return 0, 1
    a, b = fib_pair_big(n >> 1)
    c = a * ((b << 1) - a)
    d = a * a + b * b
    if n & 1:
        return d, c + d
    return c, d

def fib_big(n):
    return fib_pair_big(n)[0]

def fib_pair_mod(n, mod):
    if n == 0:
        return 0, 1 % mod
    a, b = fib_pair_mod(n >> 1, mod)
    c = (a * ((2 * b - a) % mod)) % mod
    d = (a * a + b * b) % mod
    if n & 1:
        return d, (c + d) % mod
    return c, d

def is_identity_power(k, n):
    fk, fk1 = fib_pair_mod(k, n)
    return fk == 0 and fk1 == (1 % n)

def sum_pisano_equals(d, limit):
    G = gcd(fib_big(d), fib_big(d - 1) - 1)
    ps = [v[0] for v in factorization(d)]
    ans = 0
    for n in divisors(G):
        if n == 1 or n >= limit:
            continue
        for p in ps:
            if is_identity_power(d // p, n):
                break
        else:
            ans += n
    return ans

print(sum_pisano_equals(M,N))