题解 P4000 【斐波那契数列】
飞雨烟雁
2022-04-11 13:34:37
#### 原题链接：[P4000 斐波那契数列](https://www.luogu.com.cn/problem/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$ 要怎么求？为了解决这个问题，我们可能需要一点前置知识……
------------
## 前置知识：二次剩余
------------
可能需先掌握「[原根](https://www.luogu.com.cn/problem/P6091)」的知识。
二次剩余是指这样一类问题：
> 已知 $n,p$ 的值，解方程 $x^2\equiv n\pmod p$ 。
洛谷上有模板题：[P5491 【模板】二次剩余](https://www.luogu.com.cn/problem/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$ 无解的充要条件，记为**结论二**。
事实上，上面这些对下文有用的只有这句话：
$x^2=n\pmod p$ 无解等价于 $n^{\frac{p-1}{2}}\equiv -1\pmod p$。
------------
## 斐波那契循环节
------------
接下来就是主线情节了。
**结论三**：斐波那契数列 $\{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})$ 的最小公倍数即可，证毕。
------------
## 代码实现
------------
回到原题，开始整理思路。
1. 对 $P$ 进行质因数分解，时间复杂度 $O(\sqrt P)$；
2. 对于 $P$ 的每个质因数分解 $p_i^{k_i}$，求出对应的周期长 $\pi(p_i^{k_i})$；
3. 求出每个周期的最小公倍数，作为 $\pi(P)$ 的值；
4. 把读入的 $n$ 模以 $\pi(P)$，然后矩阵快速幂解决。
就这样，代码在此：
```cpp
#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;
}
```
------------
#### 参考资料：
- [Pisano Period](https://www.cnblogs.com/cjoierShiina-Mashiro/p/11385287.html)
- [The Period of the Fibonacci Sequence Modulo j
](https://gradprogram.math.arizona.edu/~ura-reports/071/Campbell.Charles/Final.pdf)