P4000 斐波那契数列 另类做法
Itst
2020-01-29 20:13:53
同步发布于[cnblogs](https://www.cnblogs.com/Itst/p/12241180.html)
## 前言
> 你打开了“P4000 斐波那契数列”一题；
> 你发现是已经写过 $\mathrm{998244853}$ 遍的求 $\mathrm{Fib}_n$；
> 你熟练地写出矩阵快速幂并提交；
> 你得到了一版的 $\mathrm{TLE}$ ，因为 $n \leq 10^{30000000}$；
> 你点开了题解，发现第一篇题解一大片公式和定理；
> 你向下滑动界面，其他的题解都直接摆结论。
相信经历了这么多的你心中满满的“我太难了”。
那么这篇题解将拯救你于水火之中（？）
## 前置芝士
- 主角：生日悖论；
- 辅助结论：模数为 $p$ 的斐波那契循环节长度 $\pi(p) \leq 6p$；
- $O(\sqrt p)-O(1)$ 快速幂，[LOJ有模板题](https://loj.ac/problem/162)。
关于辅助结论的证明，因为前置结论比较多 ~~（实际上是我不会）~~ ，故在本篇题解中被省略。如果你有条件，可以参考 $\mathrm{Wikipedia}$ 中的 $\mathrm{Pisano\ Period}$ 条目，其中有证明。
所以实际上这个做法也有很多前置。但是相比纯数论做法，结论只有这一个，而且这个结论很好记。
## 算法流程
一件显然的事情是需要计算斐波那契数列的循环节。暴力做法基本没有优化空间，考虑：设循环节长度为 $\pi(p)$，对于 $i \neq j$ ，如果 $\mathrm{Fib}_i = \mathrm{Fib}_j$ 且 $\mathrm{Fib}_{i+1} = \mathrm{Fib}_{j+1}$，那么$\pi(p) \mid j-i$。也就是说如果找到两个位置 $i,j$ 满足条件就可以得到循环节长度的某个倍数。
因为只需求值，所以求出$\pi(p)$的真实值没有必要，求出其倍数也可以接受。所以可以利用随机化得到一个新算法：每一次随机两个数 $i,j$，计算 $\mathrm{Fib}_i , \mathrm{Fib}_{i+1} , \mathrm{Fib}_j , \mathrm{Fib}_{j+1}$ 的值判断是否相等。
显然这个做法是期望 $O(p)$ 的，跟暴力没有区别。但我们可以换一个角度描述这个问题：每一次随机两个位置 $i,j$，相当于随机两个数 $x = i \mod \pi(p)$ 和 $y = j \mod \pi(p)$，如果 $x=y$ 就可以找到循环节的倍数。
可以发现这个问题与“生日悖论”问题基本一致。利用生日悖论思想，每一次随机两个命中概率很小，但随机三个、四个、很多个位置，随着随机的数量增长，命中概率是以平方级别增长的。
所以可以得到一个更优秀的算法：使用一个哈希表记录之前随机到的所有 $(\mathrm{Fib}_i , \mathrm{Fib}_{i+1} , i)$ 三元组，每次随机位置 $j$ 并得到 $(\mathrm{Fib}_j , \mathrm{Fib}_{j+1} , j)$，如果在哈希表中可以找到满足 $\mathrm{Fib}_i = \mathrm{Fib}_j$ 且 $\mathrm{Fib}_{i+1} = \mathrm{Fib}_{j+1}$ 的三元组 $(\mathrm{Fib}_i , \mathrm{Fib}_{i+1} , i)$，就可以得到循环节的某个倍数。
根据生日悖论，上述算法在期望 $O(\sqrt p)$ 次内可以完成寻找。使用 $O(\sqrt p) - O(1)$ 矩阵/扩域快速幂实现求斐波那契数，可以做到 $O(\sqrt p)$ 求出循环节。
## 后话
这个做法是在复习`Pollard-rho`算法时 ~~打隔膜~~ 翻到有博客记录 $\pi(p) \leq 6p$ 时浮现的。如果在考场上遇到了类似问题，只需要信仰猜想循环节长度为模数级别、把 $\frac{\text{循环节长度}}{\text{模数}}$ 这个常数估计大一些然后套用本题做法就可以了，拓展性比较高。
如果想求 $\pi(p)$ 的真实值，只需要枚举算出来的循环节的约数然后逐个判断即可。
## Show me the Code.
两个实现细节：
- 建议随机位置上界大于 $12p$，否则期望次数可能会退化。代码中的随机上界是 $2^{36}$，选择二的次幂作为模数（或许）可以最小化因为取模导致的`mt19937_64`随机不均匀问题。
- $p < 2^{31}$ 所以可以直接用一个`long long`存下 $(\mathrm{Fib}_i,\mathrm{Fib}_{i+1})$ 二元组，然后就可以用 `unordered_map`实现哈希表了。
与此同时这份代码常数比较大，用扩域常数会小一半。
```cpp
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
unordered_map < ull , ll > circ; ll len;
int MOD , MX = 1 << 18; mt19937_64 rnd(time(0));
struct matrix{
ll arr[2][2];
matrix(){memset(arr , 0 , sizeof(arr));}
ll* operator [](int x){return arr[x];}
friend matrix operator *(matrix p , matrix q){
matrix x;
for(int i = 0 ; i < 2 ; ++i)
for(int j = 0 ; j < 2 ; ++j)
for(int k = 0 ; k < 2 ; ++k)
x[i][k] += p[i][j] * q[j][k];
for(int i = 0 ; i < 2 ; ++i)
for(int j = 0 ; j < 2 ; ++j) x[i][j] %= MOD;
return x;
}
}G , T[2][1 << 18 | 1];
signed main(){
static char str[300000003]; scanf("%s %d" , str + 1 , &MOD);
T[0][0][0][0] = T[0][0][1][1] = T[1][0][0][0] = T[1][0][1][1] = 1;
T[0][1][0][1] = T[0][1][1][0] = T[0][1][1][1] = 1;
for(int i = 2 ; i <= MX ; ++i) T[0][i] = T[0][i - 1] * T[0][1];
T[1][1] = T[0][MX]; for(int i = 2 ; i <= MX ; ++i) T[1][i] = T[1][i - 1] * T[1][1];
while(1){
ll x = (rnd() << 28 >> 28); matrix C = T[0][x & (MX - 1)] * T[1][x >> 18];
ull val = ((1ull * C[0][0]) << 32) | C[0][1];
if(circ.find(val) != circ.end()){len = abs(circ[val] - x); break;}
circ[val] = x;
}
ll sum = 0;
for(int i = 1 ; str[i] ; ++i) sum = (sum * 10 + str[i] - '0') % len;
cout << (T[0][sum & (MX - 1)] * T[1][sum >> 18])[0][1];
return 0;
}
```