傳送門：bzoj2432

題解

良心題，暴力有 $70pts$。。。

顯然是分別維護斐波那契數列 $\mathrm{m}$

$1,1,2,3,5,0$
$5,5,3,0$
$3,3,6,2,0$
$2,2,4,6,3,2,5,0$
$5,5,3,0$
$...$

觀察到一些特殊的性質：

• 每行開頭必然是兩個相同的數(設為 $st_i$)，且每一行形式相同： $st_i*f_1,st_i*f_2,st_i*f_3...$

• 如果有迴圈節的話，迴圈行數必然不超過 $k$。

• 迴圈節中每行結尾都是以 $st_i*f_{len}\equiv 1\mod k$結束的：
  $3,3,6,2,0$
  $2,2,4,6,3,2,5,0,5,5,3,0$
  這樣表示比較直觀。

只需要找到最小的 $f_i\equiv st_i^{-1}\mod k$就得到了第 $i$行的長度，且 $st_{i+1}=st_i*f_{len-1}\mod k$

具體來說，用 $vs_i$表示最小的模 $k$等於 $i$的斐波那契數項。而可以證明模 $k$意義下，斐波那契迴圈節長度不超過 $6*k$，且迴圈節中前三個數必然為 $0,1,1$(一個paper)

$exgcd$求解逆元，但 $k$不保證是質數，無解的情況就等同於沒有特殊的減一操作，直接矩乘。

所以我們只需要在不超過 $k$次內找出迴圈節然後矩陣快速冪即可。

行內的轉移矩陣：
$\left[ \begin{matrix} 1& 1& 0\\ 1&0&0\\ 0&0&1 \end{matrix} \right]$
行之間的轉移矩陣：
$\left[ \begin{matrix} 1& 0& 0\\ 0&1&0\\ -1&0&1 \end{matrix} \right]$
初始矩陣:
$\left[ \begin{matrix} f_0=0&f_{-1}=1 &1 \end{matrix} \right]$

程式碼

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e6+10;

ll n;bool fwg[N];
int K,mod,pr,vs[N],f[N*6],inv[N];

inline int ad(int x,int y){x+=y;return x>=mod?x-mod:x;}
inline int dc(int x,int y){x-=y;return x<0?x+mod:x;}

struct Mat{
	int g[3][3];
	inline int *operator[](int x){return g[x];}
	Mat operator *(const Mat&ky)const{
	    Mat re;register int i,j,k,v;
	    for(i=0;i<3;++i)
	     for(j=0;j<3;++j){
	     	for(v=k=0;k<3;++k) v=ad(v,(ll)g[i][k]*ky.g[k][j]%mod);
	     	re.g[i][j]=v;
	     }
	    return re;
	}
	inline void itia()
	{
		register int i,j;
		for(i=0;i<3;++i)
		 for(j=0;j<3;++j)
		  g[i][j]=(i==j);
	}
}A,B,C,D,ori,mt[N];

inline Mat fp(Mat a,ll y)
{
	ori.itia();
	for(;y;y>>=1,a=a*a)
	 if(y&1) ori=ori*a;
	return ori;
}


int exgcd(int a,int b,ll &x,ll &y)
{
	if(!b) {x=1;y=0;return a;}
	int re=exgcd(b,a%b,y,x);y-=x*(a/b);
	return re;
}

inline int gtiv(int a,int b)
{
	ll x,y;
	return exgcd(a,b,x,y)==1?((x%b+b)%b):(-1);
}

int main(){
    int i,x,z,len;ll cot;
    scanf("%lld%d%d",&n,&K,&mod);
    f[1]=f[2]=1;
	for(i=3;;++i){
		f[i]=(f[i-1]+f[i-2])%K;
		if(!vs[f[i]]) vs[f[i]]=i;
		if(f[i]==1 && f[i]==f[i-1]) break;
    }
    A[0][0]=A[0][1]=A[1][0]=A[2][2]=1;
    B[0][0]=B[1][1]=B[2][2]=1;B[2][0]=mod-1;
    C[0][1]=C[0][2]=1;
    for(z=1,pr=0;n;){
    	if(!inv[z]) inv[z]=gtiv(z,K);
    	if(inv[z]==-1) {C=C*fp(A,n);break;}
    	if(pr || (!fwg[z])){
    		len=vs[inv[z]];
			if(!len || n<len){C=C*fp(A,n);break;}
    		if(!fwg[z]){fwg[z]=true;mt[z]=fp(A,len)*B;}
    		n-=len;C=C*mt[z];z=(ll)z*f[len-1]%K;
    	}else{
    		cot=0;cot=vs[inv[z]];D=mt[z];
    		for(x=(ll)z*f[vs[inv[z]]-1]%K;x!=z;x=(ll)x*f[vs[inv[x]]-1]%K)
    		  D=D*mt[x],cot+=vs[inv[x]];
    		C=C*fp(D,n/cot);n%=cot;pr=1;
    	}
    }
    printf("%d",C[0][0]);
    return 0;
}