原题链接：P4000 斐波那契数列
（一开始以为是道水题，没想到竟然如此复杂，呜呜呜）
题目：求 f_n\bmod P 的值。
注： 原题的模数是小写的 p，此处为了区分下文的质数 p，改为大写的 P。
首先，n 达到了 10^{30000000} 的级别，矩阵快速幂肯定 TLE。
于是我们想到了要通过一些方法降低 n 的大小，最先想到的是对 n 模以某个数，然后再矩阵快速幂不就好了嘛。
那么先手算一下取模的这个数和 P 的关系吧，假设 P=4，那么斐波那契数列长这样：
1,1,2,3,1,0,1,1,2,3,1,0,...
好像确实有规律，针对 P=4 时我们只需求 f_{n\bmod 6} 就好了！
可是对于其他的 P 要怎么求？为了解决这个问题，我们可能需要一点前置知识……
前置知识：二次剩余
可能需先掌握「原根」的知识。
二次剩余是指这样一类问题：
已知 n,p 的值，解方程 x^2\equiv n\pmod p 。
洛谷上有模板题：P5491 【模板】二次剩余。此处不做详述，只讲与下文有关的，即 p 是奇素数的情况。
首先 p\mid n 时，x^2\equiv n\equiv 0\pmod p，故 x=0。
当 p\nmid n 时，x 可能无解，如 n=2,p=3 时。那么 n,p 满足什么条件时 x 有解呢？
结论一：n^{\frac{p-1}{2}}\equiv 1\pmod p 是 x 有不为 0 的解的充要条件。以下证明：
当 x 有不为 0 的解时，x^{p-1}\equiv 1。因此 n^{\frac{p-1}{2}}\equiv(x^2)^{\frac{p-1}{2}}\equiv x^{p-1}\equiv 1，必要性证毕。
当 n^{\frac{p-1}{2}}\equiv 1 时，设 g 是 p 的原根，则存在 k 使得 n=g^k，故 g^{\frac{k}{2}(p-1)}\equiv 1。因为 g 的阶是 p-1，所以 (p-1)|\frac{k}{2}(p-1)，即 k 是偶数。那么解得 x\equiv g^{\frac{k}{2}}，有不为 0 的解。充分性证毕。
那什么时候无解呢？因为 n^{p-1}\equiv 1，(n^{\frac{p-1}{2}}-1)(n^{\frac{p-1}{2}}+1)\equiv 0，即 n^{\frac{p-1}{2}}\equiv \pm1。由结论一推理得，n^{\frac{p-1}{2}}\equiv -1 是 x 无解的充要条件，记为结论二。
事实上，上面这些对下文有用的只有这句话：
------------
## 斐波那契循环节
------------
接下来就是主线情节了。
**结论三**：斐波那契数列 $\{f_i\}$ 对 $p$ 取模后必定是周期数列。
显然，只要存在 $i\neq j$ 满足 $f_i\equiv f_j$ 且 $f_{i+1}\equiv f_{j+1}$。那么由递推关系可以得到 $f_{i+k}\equiv f_{j+k}$，命题就成立了。
要证明存在这样的 $i,j$ 也不难，因为不可能有无穷多对不同的 $(f_i,f_{i+1})$，至多 $p^2$ 对，因此一定存在 $i\neq j$，满足 $f_i\equiv f_j$ 且 $f_{i+1}\equiv f_{j+1}$。故循环节一定存在，且长度不超过 $p^2$，证毕。
设 $\pi(p)$ 是这个数列的周期（注意，此处并不要求 $\pi(p)$ 是最小正周期）。那要怎样求 $\pi(p)$？
不难发现，若存在 $k$ 满足 $\forall i, f_{i}\equiv f_{i+k}\pmod p$，那 $k$ 就是 $\{f_i\bmod p\}$ 的周期，取 $\pi(p)=k$ 就好了。
------------
### Part 1：
让我们从最简单的情况开始：当模数 $P$ 等于**某个质数** $p$ 的时候。
先搬出斐波那契数列的通项公式：$f_i=\dfrac{(\frac{1+\sqrt 5}{2})^i-(\frac{1-\sqrt 5}{2})^i}{\sqrt 5}$。
因为公式里有 $\sqrt 5,\frac1 2$，所以我们先单独讨论 $p=2,5$ 的情况：经计算取 $\pi(2)=3,\pi(5)=20$。
撇除这些特殊情况后，$p$ 便是一个不等于 $5$ 的奇素数。但根据二次剩余的知识，$\sqrt 5$ 模 $p$ 时不一定能表示成整数，需分类讨论。
- $x^2\equiv 5\pmod p$ 有非 $0$ 解
这时 $x$ 等价于 $\sqrt 5$ ，那么 $\sqrt 5,\frac{1}{\sqrt 5},\frac 1 2$ 模 $p$ 就都有意义，即 $\frac{1+\sqrt 5}{2},\frac{1-\sqrt 5}{2}$ 能用整数表示，适用于欧拉定理。所以：
$$\begin{aligned}f_{i+\varphi(p)}&\equiv \frac{(\frac{1+\sqrt 5}{2})^{i+\varphi(p)}-(\frac{1-\sqrt 5}{2})^{i+\varphi(p)}}{\sqrt 5}\\&\equiv\frac{(\frac{1+\sqrt 5}{2})^i-(\frac{1-\sqrt 5}{2})^i}{\sqrt 5}=f_i\end{aligned}$$
故 $\varphi(p)$ 是数列的一个周期，取 $\pi(p)= \varphi(p)=p-1$。
- $x^2\equiv 5\pmod p$ 无解
无解意味着找不到整数替代 $\sqrt 5$，需要**扩域**。
扩域指的是令某个数 $w$ 满足 $w^2\equiv K$，其中 $K^{\frac{p-1}{2}}\equiv -1$。因为 $x^2\equiv K$ 无解，所以这里的 $w$ 类似于 $i^2=-1$ 的 $i$，并非整数。同时，由 $w$ 构成的数 $a+bw$ 也不满足欧拉定理等，需要新的结论。
**结论四**：对于奇素数 $p$，有 $(a+bw)^{p+1}\equiv a^2-b^2K\pmod p$。
以下证明：由二项式定理可得：
$$(a+bw)^{p+1}=\sum_{i=0}^{p+1}C_{p+1}^ia^{p+1-i}(bw)^i$$
聚焦于组合数 $C_{p+1}^i\equiv\frac{p!}{i!(p+1-i)!}$，可以发现仅在 $i=0,1,p,p+1$ 这四个值时组合数模 $p$ 后不为 $0$。因此原式等于：
$$C_{p+1}^0a^{p+1-0}(bw)^0+C_{p+1}^1a^{p+1-1}(bw)^1+C_{p+1}^{p}a^{p+1-p}(bw)^p+C_{p+1}^{p+1}a^{p+1-p-1}(bw)^{p+1}$$
即 $a^{p+1}+a^{p}bw+ab^pw^p+b^{p+1}w^{p+1}\equiv a^2+abw+abw^p+b^2w^{p+1}$。
$\because w^2\equiv K,K^{\frac{p-1}{2}}\equiv -1
\therefore (w^2)^{\frac{p-1}{2}}\equiv w^{p-1}\equiv -1
代入原式得：a^2-b^2w^2\equiv a^2-b^2K，证毕。
接下来，设 w^2\equiv 5，则通项公式可写成：f_i=\dfrac{(\frac{p+1}{2}+\frac{p+1}{2}w)^i-(\frac{p+1}{2}-\frac{p+1}{2}w)^i}{w}。
先看 \frac{p+1}{2}+\frac{p+1}{2}w，由结论四得：(\frac{p+1}{2}+\frac{p+1}{2}w)^{p+1}\equiv (\frac{p+1}{2})^2-5(\frac{p+1}{2})^2\equiv -1。
所以 (\frac{p+1}{2}+\frac{p+1}{2}w)^{2p+2}\equiv 1，同理 (\frac{p+1}{2}-\frac{p+1}{2}w)^{2p+2}\equiv 1。
因此 f_{i+2p+2}\equiv f_i，取 \pi(p)=2p+2。
Part 2：
接下来，讨论模数 P 是某个质数的幂的情况，记 P=p^k。
注：Part 2 的证明参考了部分资料（链接见文末）。
结论五：对于质数 p，已知 a\equiv 1\pmod p，则 a^{p^k}\equiv 1\pmod {p^{k+1}}。
以下证明：考虑数学归纳法，假设 k 时等式成立，证明 k+1 时等式也成立。
设 a^{p^{k}}=sp^{k+1}+1，那么 a^{p^{k+1}}=(a^{p^k})^p=(sp^{k+1}+1)^p。
由二项式定理：
(sp^{k+1}+1)^p=\sum_{i=0}^pC_p^i(sp^{k+1})^i
当 i\ge 2 时，p^{(k+1)i}\bmod p^{k+2}=0，只需算出 i=0,1 的值。代入得 1+sp^{k+2}\equiv 1\pmod {p^{k+2}}，证毕。
方便起见，记 \frac{1+\sqrt 5}{2}=u,\frac{1-\sqrt 5}{2}=v。
结论六：对于质数 p，有 u^{\pi(p)}\equiv v^{\pi(p)}\equiv 1\pmod p。
\because f_{\pi(p)}\equiv 0\pmod p
\therefore u^{\pi(p)}\equiv v^{\pi(p)}
\because f_{\pi(p)+1}\equiv 1
\therefore {u^{\pi(p)}u-v^{\pi(p)}v}\equiv \sqrt 5
代入 u^{\pi(p)}\equiv v^{\pi(p)} 化简得：u^{\pi(p)}\equiv v^{\pi(p)}\equiv 1，证毕。
结合结论五和结论六可得：u^{\pi(p)p^{k-1}}\equiv v^{\pi(p)p^{k-1}}\equiv 1\pmod {p^k}。
因此 f_{\pi(p)p^{k-1}}\equiv \frac{1}{\sqrt 5}(u^{\pi(p)p^{k-1}}-v^{\pi(p)p^{k-1}})\equiv 0\pmod {p^k}。
且 f_{1+\pi(p)p^{k-1}}\equiv \frac{1}{\sqrt 5}(u^{\pi(p)p^{k-1}}u-v^{\pi(p)p^{k-1}}v)\equiv \frac{u-v}{\sqrt 5}\equiv 1\pmod {p^k}。
也就是说：f_{i+\pi(p)p^{k-1}}\equiv f_i\pmod {p^k}，取 \pi(p^k)=\pi(p)p^{k-1}。
Part 3：
最后一步了，也是最简单的情况。
结论七：设模数 P 的质因数分解式为 \prod p^k，则可取 \pi(P)=\operatorname{lcm}(\pi(p_i^{k_i}))。
以下证明：
\because f_{\pi(P)}\equiv 0\pmod{p_i^{k_i}},p_i^{k_i}\mid f_{\pi(P)}
\therefore {\pi(p_i^{k_i})}\mid {\pi(P)}
取 \pi(P) 为 \pi(p_i^{k_i}) 的最小公倍数即可，证毕。
代码实现
回到原题，开始整理思路。
-
对 P 进行质因数分解，时间复杂度 O(\sqrt P)；
-
对于 P 的每个质因数分解 p_i^{k_i}，求出对应的周期长 \pi(p_i^{k_i})；
-
求出每个周期的最小公倍数，作为 \pi(P) 的值；
-
把读入的 n 模以 \pi(P)，然后矩阵快速幂解决。
就这样，代码在此：
#include <iostream>
#include <cstdio>
#include <cstring>
#define ll long long
using namespace std;
ll FastPow(ll a, ll b, ll mod){
ll res = 1;
while(b){
if(b & 1) res = res * a % mod;
b >>= 1, a = a * a % mod;
}
return res;
}
ll FastPowNoMod(ll a, ll b){
ll res = 1;
while(b){
if(b & 1) res = res * a;
b >>= 1, a = a * a;
}
return res;
}
const int Mx = 3e7 + 10;
int Fac[30], Tim[30], tot;
char N[Mx];
void Divide(int p){ // 分解质因数
for(int i = 2; i * i <= p; ++i) if(p % i == 0){
Fac[++tot] = i;
while(p % i == 0) p /= i, ++Tim[tot];
}
if(p != 1) Fac[++tot] = p, Tim[tot] = 1;
}
ll PrimeLoop(ll p){ // 求 pi(p)
if(p == 2) return 3;
if(p == 5) return 20;
if(FastPow(5, (p - 1) >> 1, p) == 1) return p - 1;
return 2 * p + 2;
}
ll PrimePow(int id){ return FastPowNoMod(Fac[id], Tim[id] - 1) * PrimeLoop(Fac[id]);} // 求 pi(p^k)
ll Gcd(ll a, ll b){ return b == 0 ? a : Gcd(b, a % b);}
ll GetLoop(){ // 求 pi(P)
ll res = PrimePow(1), temp;
for(int i = 2; i <= tot; ++i) temp = PrimePow(i), res = res / Gcd(res, temp) * temp;
return res;
}
int P;
struct Matrix{
ll val[3][3];
Matrix operator * (const Matrix& a) const {
Matrix Res = (*this);
for(int i = 1; i <= 2; ++i) for(int j = 1; j <= 2; ++j)
Res.val[i][j] = (val[i][1] * a.val[1][j] % P + val[i][2] * a.val[2][j] % P) % P;
return Res;
}
Matrix operator *= (const Matrix& a) { return (*this) = a * (*this);}
Matrix operator ^ (ll k) const {
Matrix Res, Temp = (*this);
Res.val[1][1] = Res.val[2][2] = 1, Res.val[1][2] = Res.val[2][1] = 0;
while(k){
if(k & 1) Res *= Temp;
k >>= 1, Temp *= Temp;
}
return Res;
}
} Fibo, Cell;
void Solve(ll n){
if(n == 0){
printf("0");
return;
}
Fibo.val[1][1] = 0, Fibo.val[2][1] = Fibo.val[1][2] = Fibo.val[2][2] = 1;
Cell.val[1][1] = 1, Cell.val[1][2] = 1;
Matrix Res = Cell * (Fibo ^ (n - 1));
printf("%lld", Res.val[1][1] % P);
}
int main(){
scanf(" %s %d", N, &P);
if(P == 1){
printf("0");
return 0;
}
Divide(P);
ll Loop = GetLoop(), n = 0;
int Len = strlen(N);
for(int i = 0; i < Len; ++i) n = (n * 10 + (N[i] ^ 48)) % Loop;
Solve(n);
return 0;
}
参考资料：