斐波那契数模p的循环节+矩阵快速幂
#include <iostream>
#include <string.h>
#include <algorithm>
#include <stdio.h>
#include <math.h>
using namespace std;
typedef long long LL;
const int M = 2;
struct Matrix
{
LL m[M][M];
};
Matrix A;
Matrix I = {1,0,0,1};
LL loop[10005];
Matrix multi(Matrix a,Matrix b,LL MOD)
{
Matrix c;
c.m[0][0]=(a.m[0][0]*b.m[0][0]%MOD+a.m[0][1]*b.m[1][0]%MOD)%MOD;
c.m[0][1]=(a.m[0][0]*b.m[0][1]%MOD+a.m[0][1]*b.m[1][1]%MOD)%MOD;
c.m[1][0]=(a.m[1][0]*b.m[0][0]%MOD+a.m[1][1]*b.m[1][0]%MOD)%MOD;
c.m[1][1]=(a.m[1][0]*b.m[0][1]%MOD+a.m[1][1]*b.m[1][1]%MOD)%MOD;
return c;
}
Matrix power(Matrix a,LL k,LL MOD)
{
Matrix ans = I,p = a;
while(k)
{
if(k & 1)
{
ans = multi(ans,p,MOD);
k--;
}
k >>= 1;
p = multi(p,p,MOD);
}
return ans;
}
LL gcd(LL a,LL b)
{
return b==0? a:gcd(b,a%b);
}
const int N = 10005;
const int NN = 5005;
LL num[NN],pri[NN];
LL fac[NN];
int cnt,c;
bool prime[N];
int p[NN];
int k;
void isprime()
{
k = 0;
memset(prime,true,sizeof(prime));
for(int i=2;i<N;i++)
{
if(prime[i]) {
p[k++]=i;
}
for(int j=0;j<k;j++)
{
if(i*p[j]>N) break;
prime[i*p[j]]=false;
if(i%p[j]==0) break;
}
}
}
LL quick_mod(LL a,LL b,LL m)
{
LL ans = 1;
a %= m;
while(b)
{
if(b & 1)
{
ans = ans * a % m;
b--;
}
b >>= 1;
a = a * a % m;
}
return ans;
}
//判断勒让德符号,判断a在模p意义下是否为二次剩余
LL legendre(LL a,LL p)
{
if(quick_mod(a,(p-1)>>1,p)==1) return 1;
else return -1;
}
//O(sqrt(n))求解唯一分解式
void Solve(LL n)
{
cnt = 0;
LL t = (LL)sqrt(1.0*n);
int f=lower_bound(p,p+k,n)-p;
if(f<k&&p[f]==n) {
pri[cnt]=n;
num[cnt]=1;
cnt++;
return;
}
for(int i=0; p[i]<=t; i++)
{
if(n%p[i]==0)
{
int a = 0;
pri[cnt] = p[i];
while(n%p[i]==0)
{
a++;
n /= p[i];
}
num[cnt] = a;
cnt++;
}
if(n==1) break;
}
if(n > 1)
{
pri[cnt] = n;
num[cnt] = 1;
cnt++;
}
}
//寻找因子
void Work(LL n)
{
c = 0;
LL t = (LL)sqrt(1.0*n);
for(int i=1; i<=t; i++)
{
if(n % i == 0)
{
if((LL)i * i == n) fac[c++] = i;
else
{
fac[c++] = i;
fac[c++] = n / i;
}
}
}
}
LL find_loop(LL n)
{
Solve(n);
LL ans=1;
for(int i=0; i<cnt; i++)
{
LL record=1;
if(pri[i]==2)
record=3;
else if(pri[i]==3)
record=8;
else if(pri[i]==5)
record=20;
else
{
if(legendre(5,pri[i])==1)
Work(pri[i]-1);
else
Work(2*(pri[i]+1));
sort(fac,fac+c);
//找到因子之后,再一个一个因子进行矩阵快速幂求解
for(int k=0; k<c; k++)
{
Matrix a = power(A,fac[k]-1,pri[i]);
LL x = (a.m[0][0]%pri[i]+a.m[0][1]%pri[i])%pri[i];
LL y = (a.m[1][0]%pri[i]+a.m[1][1]%pri[i])%pri[i];
if(x==1 && y==0)
{
record = fac[k];
break;
}
}
}
for(int k=1; k<num[i]; k++)
record *= pri[i];
ans = ans/gcd(ans,record)*record;
}
return ans;
}
void Init()
{
A.m[0][0] = 1;
A.m[0][1] = 1;
A.m[1][0] = 1;
A.m[1][1] = 0;
}
int main()
{
LL n,k1,p1;
Init();
isprime();
int T;
scanf("%d",&T);
int kase=0;
while(T--)
{
scanf("%lld%lld%lld",&n,&k1,&p1);
printf("Case #%d: ",++kase);
loop[1]=p1;
for(int i=2;i<=k1+1;i++)
loop[i]=find_loop(loop[i-1]);
LL solve=n;
Matrix ans;
bool flag=false;
for(int i=k1+1;i>=1;i--)
{
if(flag) break;
if(solve>=4) {
ans=power(A,solve-1,loop[i]);
solve=(ans.m[0][0]+ans.m[0][1])%loop[i];
}
else if(solve==0) {
printf("1\n");
flag=true;
}
else if(solve==1) {
printf("1\n");
flag=true;
}
else if(solve==2) {
printf("2\n");
flag=true;
}
else if(solve==3) {
printf("3\n");
flag=true;
}
}
if(!flag) printf("%lld\n",solve);
}
return 0;
}
需要先找到各种循环节,再进行矩阵快速幂
//求循环节+矩阵快速幂
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod1=1e9+7;
const LL mod2=222222224;
const LL mod3=183120;
struct Matrix{
LL m[2][2];
};
Matrix mul(Matrix a,Matrix b,LL MD)
{
Matrix c;
c.m[0][0]=((a.m[0][0]*b.m[0][0])%MD+(a.m[0][1]*b.m[1][0]%MD))%MD;
c.m[0][1]=((a.m[0][0]*b.m[0][1])%MD+(a.m[0][1]*b.m[1][1]%MD))%MD;
c.m[1][0]=((a.m[1][0]*b.m[0][0])%MD+(a.m[1][1]*b.m[1][0]%MD))%MD;
c.m[1][1]=((a.m[1][0]*b.m[0][1])%MD+(a.m[1][1]*b.m[1][1]%MD))%MD;
return c;
}
Matrix pow(Matrix A,LL k,LL MD)
{
Matrix ans;
ans.m[0][0]=1,ans.m[0][1]=0,ans.m[1][0]=0,ans.m[1][1]=1;
while(k)
{
if(k&1) {
ans=mul(ans,A,MD);
k--;
}
k>>=1;
A=mul(A,A,MD);
}
return ans;
}
int main()
{
LL n;
while(scanf("%lld",&n)==1)
{
if(n==0) {
printf("0\n");
continue;
}
if(n==1) {
printf("1\n");
continue;
}
Matrix A,ans;
A.m[0][0]=3,A.m[0][1]=1,A.m[1][0]=1,A.m[1][1]=0;
LL solve=n;
if(solve>=2) {
ans=pow(A,solve-1,mod3);
solve=ans.m[0][0]%mod3;
}
if(solve>=2) {
ans=pow(A,solve-1,mod2);
solve=ans.m[0][0]%mod2;
}
if(solve>=2) {
ans=pow(A,solve-1,mod1);
solve=ans.m[0][0]%mod1;
}
cout<<solve<<endl;
}
return 0;
}