求 Fib 数模 n 的循环节:
1. 对 n 做因数分解:
n=p1^e1 * p2^e2 * … * pt^et;
2. 求出每个素数 pi 对应 Fib 数模 pi 的循环节mi0 ,则 pi^ei 对应的 Fib 数模 pi^ei 的 循环节 mi=mi0 * pi^(ei-1);
3. Fib 数模 n 的循环节就等于 lcm(mi)。
关键在于如何求Fib 数模素数 pi 的循环节mi0,有以下结论,如果 pi 是5的二次剩余, mi0 是 pi-1 的约数; 如果 pi 不是,则 mi0 是 2(pi+1) 的约数。
如果想了解上述方法证明过程可以看一下这篇文章
Charles W. Campbell II. The Period of the Fibonacci Sequence Modulo j. 2007.
以下是扩展
广义斐波那契数列
但这道题数据量非常大,很有可能会超时,所以需要一些技巧。
1. 将取模的操作尽量用加减法代替:
比如说 k = (a%p + b%p)%p 就可以写成
if((k = a%p + b%p)>=p) k -= p;
2. 筛出一些质数分解n的复杂度可以由O(sqrt(n))变成O(sqrt(n)/ln(n)):
3. 算模质数p意义下的循环节时,由于斐波那契数列的”循环节的倍数”次项在模意义下相等,所以可以用类似计算欧拉函数的形式将时间复杂度由O(sqrt(n))变成不到O(log(n)log(n)):
首先: 设 L 是模质数p意义下的循环节; 如果 p 是5的二次剩余,M = L 是 p-1 的约数,否则 M = L 是 2(p+1) 的约数。
可以设 L = p1^k1 * p2^k2 * … * pt^kt
M= p1^e1 * p2^e2 * … * pt^et * … * ps^es
由于如果 L|k ,则 k 也是循环节(不一定是最小,但求循环节默认是求最小) ,所以检查 M/(p1^j) 是否满足 (用 4 改进后的快速幂求) f(p1^j)=0 , f(p1^j +1) =1,如果 j 满足而 j - 1 不满足,则一定有 k1=j。依次类推可以将 k1,k2,…,kt求出。最差的情况是要试 (e1+e2+…+es+s) 次,为 log(M)级别,检查时也是 (log) 级别(或者更小)。
4. 斐波那契数列是线性递推数列,采用的公式Fib(n+m) = Fib(n-1) * Fib(m)+Fib(n) * Fib(m+1):
A={
{1,1},{1,0}}, {Fib(n+1),Fib(n)} = A^n * {Fib(1),Fib(0)}。
可以发现 其实 A^n = {
{Fib(n+1),Fib(n)},{Fib(n),Fib(n-1)}} = {
{Fib(n) + Fib(n-1),Fib(n) },{Fib(n),Fib(n-1)}} ,也就是说每次快速幂做矩阵乘法时没必要像普通矩阵乘法那样算, 比如算 A^(n+m) 本来是 由矩阵 A^n 和 A^m 相乘,现在就可以简化为 求 Fib(n+m) = Fib(n-1) * Fib(m)+Fib(n) * Fib(m+1) 和 Fib(n+m-1) = Fib(n-1) * Fib(m-1)+Fib(n) * Fib(m),(这些值都在 矩阵 A^n 和 A^m 中) 然后 A^(n+m) ={
{Fib(n+m) + Fib(n+m-1),Fib(n+m) },{Fib(n+m),Fib(n+m-1)}}
原先利用2 * 2的矩阵进行快速幂,每次的矩阵乘法是8次加法、8次乘法,但是使用这种改进后,3次加法、4次乘法,可以节省一半时间。
5. 打表(我打了1000以内的素数的情况)减少小质数情况的计算,其实最后这个技巧不加也已经能过了。
查阅博客 http://blog.csdn.net/acdreamers/article/details/10983813
http://blog.csdn.net/zhuangmezhuang/article/details/52627308
并感谢糖老师的讲解。
(这里打表不是必要的)
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <stack>
#include <vector>
#include <cstring>
#include <queue>
#define msc(X) memset(X,-1,sizeof(X))
#define ms(X) memset(X,0,sizeof(X))
typedef long long LL;
using namespace std;
int prime[3440];
bool notprime[32000];
void getPrime(void)
{
ms(notprime);
for(int i=2;i<32000;i++)
{
if(!notprime[i]) prime[++prime[0]]=i;
for(int j=1;j<=prime[0]&&prime[j]*i<32000;j++)
{
notprime[prime[j]*i]=true;
if(i%prime[j]==0) break;
}
}
}
int factor[30][2];
int cnt;
void getFactors(int n)
{
cnt=0;
int tmp=n;
for(int i=1;prime[i]<=tmp/prime[i];i++)
{
factor[cnt][1]=0;
if(tmp%prime[i]==0){
factor[cnt][0]=prime[i];
while(tmp%prime[i]==0){
factor[cnt][1]++;
tmp/=prime[i];
}
cnt++;
}
}
if(tmp!=1){
factor[cnt][0]=tmp;
factor[cnt++][1]=1;
}
}
LL gcd(LL a,LL b)
{
return b?gcd(b,a%b):a;}
inline void multiply(LL a[2][2],LL b[2][2],int p)
{
LL tmp01,tmp11;
tmp01=a[0][0]*b[0][1]%p+a[0][1]*b[1][1]%p;
if(tmp01>=p) tmp01-=p;
tmp11=a[1][0]*b[0][1]%p+a[1][1]*b[1][1]%p;
if(tmp11>=p) tmp11-=p;
if((a[0][0]=tmp01+tmp11)>=p) a[0][0]-=p;
a[0][1]=a[1][0]=tmp01;
a[1][1]=tmp11;
}
bool check(int n,int p)
{
LL a[2][2]={
{
1,1},{
1,0}},res[2][2]={
{
1,0},{
0,1}};
while(n){
if(n&1) multiply(res,a,p);
n>>=1;
multiply(a,a,p);
}
return res[0][0]==1&&res[1][0]==0;
}
int biao[1001]={
0,3,8,20,16,10,28,36,18,48,14,30,76,40,88,32,108,58,60,136,70,148,78,168,44,196,50,208,