817
HDU 3978 斐波那契循環節
題意：給出f(f(f...f(n)...)) 總共嵌套k次。問最後模p的值是多少。
首先應該明白的是這個題有循環節的。一個數模N的循環節就是這個數分解成素因子乘積的形式p1^a1*p2^a2*p3^a3...後，斐波那契模pi^ai的循環節的最大公約數。
那麼一個素數的k次冪的循環節=斐波那契模上這個素數的循環節乘上p^(k-1)。
而一個素數p的循環節 如果p>5並且是5的二次剩餘，那麼循環節就是(p-1)的因子，否則就是2*(p+1)的因子。所以2 3 5 的時候需要特判一下。
知道這些就能求每一次嵌套的循環節了，通過矩陣連乘即可得出答案。
#include <iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
#define maxn 20050
bool isprime[maxn];
int prime[maxn],nprime,mod[maxn],M;
int primeloop[maxn];
void getprime()
{
long long i,j;
memset(isprime,1,sizeof(isprime));
nprime=0;
for(i=2; i<maxn; i++)
if(isprime[i])
{
prime[nprime++]=i;
for(j=i*i; j<maxn; j+=i) isprime[j]=0;
}
}
int factor[100][2],tol,fac[100][2],numfac;
void findfac(int n,int f[][2],int &t)
{
int x=n,l=(int)sqrt(1.0*n);
t=0;
for(int i=0; i<50; i++) f[i][0]=f[i][1]=0;
for(int i=0; prime[i]<=l; i++)
if(x%prime[i]==0)
{
f[t][0]=prime[i];
while(x%prime[i]==0) f[t][1]++,x/=prime[i];
t++;
}
if(x>1) f[t][0]=x,f[t++][1]++;
}
const int MAX=2;
typedef struct
{
long long m[MAX][MAX];
} Matrix;
Matrix P,I;
Matrix matrixmul(Matrix a,Matrix b) //矩陣乘法
{
int i,j,k;
Matrix c;
for (i=0; i<MAX; i++)
for (j=0; j<MAX; j++)
{
c.m[i][j]=0;
for(k=0; k<MAX; k++)
c.m[i][j]+=((a.m[i][k]%M)*(b.m[k][j]%M))%M;
c.m[i][j]%=M;
}
return c;
}
void quickpow(int n,int &x,int &y)
{
Matrix m=P,b=I;
while(n)
{
if(n&1) b=matrixmul(b,m);
n>>=1,m=matrixmul(m,m);
}
x=b.m[0][0],y=b.m[1][0];
}
int gcd(int a,int b)
{
return b?gcd(b,a%b):a;
}
int exp_mod(int a,int b,int c)
{
int ans=1;
a%=c;
while(b)
{
if(b&1) ans=ans*a%c;
b>>=1,a=a*a%c;
}
return ans;
}
int minloop,ff[maxn],numff;
void dfs(int num,int s=1)
{
if(num==numfac)
{
ff[numff++]=s;
return;
}
for(int i=0; i<=fac[num][1]; i++)
dfs(num+1,s),s*=fac[num][0];
}
int getPrimeLoop(int p)
{
if(p==2) return 3;
if(p==3) return 8;
if(p==5) return 20;
M=p;
if(exp_mod(5,(p-1)>>1,p)==1) p--;
else p=2*p+2;
findfac(p,fac,numfac);
minloop=1e9;
numff=0;
dfs(0,1);
sort(ff,ff+numff);
int x,y;
for(int i=1; i<numff; i++)
{
quickpow(ff[i]-1,x,y);
if(x==0&&y==1) return ff[i];
}
}
int getLoop(int p,int k)
{
int ret;
if(p>19583)
ret=getPrimeLoop(p);
else
ret=primeloop[p];
for(int i=0; i<k-1; i++) ret*=p;
return ret;
}
int getmod(int n)
{
findfac(n,factor,tol);
int ret=1,tem;
for(int i=0; i<tol; i++)
tem=getLoop(factor[i][0],factor[i][1]),ret=ret/gcd(ret,tem)*tem;
return ret;
}
int main()
{
int t,ca=0,p,n,k,x,y;
getprime();
P.m[0][0]=P.m[0][1]=P.m[1][0]=1,P.m[1][1]=0;
I.m[0][0]=I.m[1][1]=1,I.m[1][0]=I.m[0][1]=0;
for(int i=0; i<2250; i++)
primeloop[prime[i]]=getPrimeLoop(prime[i]);
scanf("%d",&t);
while(t--)
{
scanf("%d%d%d",&n,&k,&p);
mod[0]=p;
for(int i=1; i<=k; i++) mod[i]=getmod(mod[i-1]);
for(int i=k; i>=0; i--)
{
if(i<k) n%=mod[i+1];
M=mod[i],quickpow(n,x,y),n=x;
}
printf("Case #%d: %d\n",++ca,n);
}
return 0;
}
最後更新：2017-04-03 15:22:13