LG4000 斐波那契数列 矩阵快速幂、随机化

前言

你打开了“P4000 斐波那契数列”一题;

你发现是已经写过 \(\mathrm{998244853}\) 遍的求 \(\mathrm{Fib}_n\)

你熟练地写出矩阵快速幂并提交;

你得到了一版的 \(\mathrm{TLE}\) ,因为 \(n \leq 10^{30000000}\)

你点开了题解,发现第一篇题解一大片公式和定理;

你向下滑动界面,其他的题解都直接摆结论。

相信经历了这么多的你心中满满的“我太难了”。

那么这篇题解将拯救你于水火之中(?)

前置芝士

  • 主角:生日悖论;
  • 辅助结论:模数为 \(p\) 的斐波那契循环节长度 \(\pi(p) \leq 6p\)
  • \(O(\sqrt p)-O(1)\) 快速幂,LOJ有模板题

关于辅助结论的证明,因为前置结论比较多 (实际上是我不会) ,故在本篇题解中被省略。如果你有条件,可以参考 \(\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实现哈希表了。

与此同时这份代码常数比较大。

#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;
}
posted @ 2020-01-29 20:17  cjoier_Itst  阅读(751)  评论(0编辑  收藏  举报