题解

模型与性质

记 $g_i$ 为游戏规则下第 $i$ 个月的兔子对数，$F_i$ 为普通斐波那契数。根据题意：

1. $g_1 = g_2 = 1$；
2. $g_i = g_{i-1} + g_{i-2}$，若该和在模 $k$ 意义下等于 $1$，则需要把刚出生的一对兔子删掉，即再减去 $1$。

递推是线性的，因此“在第 $t$ 个月减去 $1$” 对后续所有月份的影响都是固定的，可以看作把基准的斐波那契解乘上一个系数再减去常数。设当前月份的真实数量满足 $g_i \equiv x \cdot F_i \pmod{k}$（初始 $x=1$，之后一段时间都保持不变），那么下一次需要减一的月份就是第一个满足

$$x \cdot F_t \equiv 1 \pmod{k}$$

的 $t$。只要预先把 $F_i \bmod k$ 的第一次出现位置存入 idx（Pisano 圈的长度最多 $6k$），就能在 $O(1)$ 时间得到这一段的长度 len=t。如果当前的 $x$ 在模 $k$ 下不可逆，或者 $1/x$ 对应的余数从未在斐波那契模数列中出现，则说明之后再也不会触发“减一”，后续月份就等价于普通斐波那契数列。

完成一次“减一”后，对后续月份的影响等价于把系数更新为

$$x \leftarrow x \cdot F_{t-1} \bmod k$$

（可由 $[g_{t-1}, g_t]^T = x \cdot [F_{t-1}, F_t]^T - [0,1]^T$ 推出）。于是整体过程被分解成若干段，其结构完全由当前的 $x$ 决定。持续迭代下去要么遇到前述“再无删减”的终止情形，要么某个 $x$ 再次出现，从而形成循环。

快速求值

设 3×3 斐波那契转移矩阵为

$$A = \begin{bmatrix}1 & 1 & 0\\1 & 0 & 0\\0 & 0 & 1\end{bmatrix}, $$

向量 $[g_i, g_{i-1}, 1]^T$ 即可通过乘以 $A$ 完成一次转移。对长度为 len 的一段，直接把 $A$ 快速幂到 len 次即可一次性跳过，随后再把 $g$ 减去 $1$。这样就能在线性段之间跳跃。若存在循环，把循环内所有段的矩阵相乘得到 $B$，再用快速幂把 $B$ 提升到 $n / (\text{loop\_len})$ 次即可跳过若干完整循环，最后处理余下的段。

整个过程中：

• 预处理斐波那契模 $k$ 的最短出现位置及 get_fib 都只需 $O(k)$；
• 模拟段落、构造循环长度不超过 $k$，每段只涉及常数次矩阵乘法；
• 对 $3\times 3$ 的矩阵做快速幂，单次复杂度为 $O(\log n)$。

因此总复杂度约为 $O(k \log k + \log n)$，可以覆盖 $n \le 10^{18}, k \le 10^6$ 的数据范围。

参考代码

#include<bits/stdc++.h>
using namespace std;
const int K=1000005;
int idx[K],vis[K];
struct AAA{
	int len,x;
}hhy[K];
int gcd(int a,int b){return (b==0?a:gcd(b,a%b));}
void ex_gcd(int a,int b,int& x,int& y){
	if(b==0){
		x=1;
		y=0;
		return;
	}
	ex_gcd(b,a%b,y,x);
	y-=a/b*x;
}
int inv(int a,int p){
	if(gcd(a,p)==1){
		int x,y;
		ex_gcd(a,p,x,y);
		return (x%p+p)%p;
	} 
	return -1;
}
class mat{
	private:
		int a[3][3];
	public:
		mat():a(){memset(a,0,sizeof(a));}
		int* operator[](int i){return a[i];}
}ara,brb;
mat mul(mat a,mat b,int mod){
	mat c;
	for(int i=0;i<3;++i){
		for(int j=0;j<3;++j){
			for(int k=0;k<3;++k){
				c[i][j]=(c[i][j]+1ll*a[i][k]*b[k][j]%mod)%mod;
			}
		}
	}
	return c;
}
mat qpow(mat a,long long b,int mod){
	mat ans;
	ans[0][0]=ans[1][1]=ans[2][2]=1;
	while(b>0){
		if(b&1) ans=mul(ans,a,mod);
		b>>=1;
		a=mul(a,a,mod);
	}
	return ans;
}
int get_fib(long long pos,int mod){
	if(pos==1 || pos==2) return 1;
	mat b;
	b[0][0]=b[1][0]=1;
	b[2][0]=mod-1;
	return (mul(qpow(ara,pos-2,mod),b,mod))[0][0];
}
int main(){
	ios::sync_with_stdio(false);cin.tie(nullptr);cout.tie(nullptr);
	memset(idx,255,sizeof(idx));
	memset(vis,255,sizeof(vis));
	ara[0][0]=ara[0][1]=ara[1][0]=ara[2][2]=1;
	brb[0][0]=brb[0][2]=brb[1][1]=brb[2][2]=1;
	long long n;
	int k,p;
	cin >> n >> k >> p;
	if(n==1 || n==2){
		cout << "1\n";
		return 0;
	}
	int f1=1,f2=1;
	for(int i=3;;++i){
		int f3=(f1+f2)%k;
		if(idx[f3]==-1) idx[f3]=i;
		if(f1==0 && f2==1 && f3==1) break;
		f1=f2;
		f2=f3;
	}
	int x=1;
	bool tag=false;
	int tot=0,cir=-1;
	while(true){
		const int Inv=inv(x,k);
		if(Inv==-1){
			tag=true;
			break;
		}
		int pos=idx[Inv];
		if(pos==-1){
			tag=true;
			break;
		}
		hhy[tot]=AAA{pos,x};
		vis[x]=tot++;
		x=1ll*x*get_fib(pos-1,k)%k;
		if(vis[x]!=-1){
			cir=vis[x];
			break;
		}
	}
	if(tag){
		mat aa;
		aa[0][0]=aa[1][0]=1;
		aa[2][0]=p-1;
		for(int i=0;i<tot;++i){
			if(hhy[i].len>=n){
				if(i==0) aa=mul(qpow(ara,n-2,p),aa,p);
				else aa=mul(qpow(ara,n,p),aa,p);
				if(n==hhy[i].len) aa[0][0]=(aa[0][0]+p-1)%p;
				cout << aa[0][0] << endl;
				return 0;
			}
			n-=hhy[i].len;
			if(i==0) aa=mul(qpow(ara,hhy[i].len-2,p),aa,p);
			else aa=mul(qpow(ara,hhy[i].len,p),aa,p);
			aa[0][0]=(aa[0][0]+p-1)%p;
		}
		cout << mul(qpow(ara,n,p),aa,p)[0][0] << endl;
	}else{
		mat aa;
		aa[0][0]=aa[1][0]=1;
		aa[2][0]=p-1;
		for(int i=0;i<tot;++i){
			if(hhy[i].len>=n){
				if(i==0) aa=mul(qpow(ara,n-2,p),aa,p);
				else aa=mul(qpow(ara,n,p),aa,p);
				if(n==hhy[i].len) aa[0][0]=(aa[0][0]+p-1)%p;
				cout << aa[0][0] << endl;
				return 0;
			}
			n-=hhy[i].len;
			if(i==0) aa=mul(qpow(ara,hhy[i].len-2,p),aa,p);
			else aa=mul(qpow(ara,hhy[i].len,p),aa,p);
			aa[0][0]=(aa[0][0]+p-1)%p;
		}
		long long sum=0;
		for(int i=cir;i<tot;++i) sum+=hhy[i].len;
		if(n<sum){
			for(int i=cir;i<tot;++i){
				if(hhy[i].len>=n){
					if(i==0) aa=mul(qpow(ara,n-2,p),aa,p);
					else aa=mul(qpow(ara,n,p),aa,p);
					if(n==hhy[i].len) aa[0][0]=(aa[0][0]+p-1)%p;
					cout << aa[0][0] << endl;
					return 0;
				}
				n-=hhy[i].len;
				if(i==0) aa=mul(qpow(ara,hhy[i].len-2,p),aa,p);
				else aa=mul(qpow(ara,hhy[i].len,p),aa,p);
				aa[0][0]=(aa[0][0]+p-1)%p;
			}
			return 0;
		}
		mat bb;
		bb[0][0]=bb[1][1]=bb[2][2]=1;
		for(int i=cir;i<tot;++i){
			if(i==0) bb=mul(qpow(ara,hhy[i].len-2,p),bb,p);
			else bb=mul(qpow(ara,hhy[i].len,p),bb,p);
			bb=mul(brb,bb,p);
		}
		aa=mul(qpow(bb,n/sum,p),aa,p);
		n%=sum;
		for(int i=cir;i<tot;++i){
			if(hhy[i].len>=n){
				if(i==0) aa=mul(qpow(ara,n-2,p),aa,p);
				else aa=mul(qpow(ara,n,p),aa,p);
				if(n==hhy[i].len) aa[0][0]=(aa[0][0]+p-1)%p;
				cout << aa[0][0] << endl;
				return 0;
			}
			n-=hhy[i].len;
			if(i==0) aa=mul(qpow(ara,hhy[i].len-2,p),aa,p);
			else aa=mul(qpow(ara,hhy[i].len,p),aa,p);
			aa[0][0]=(aa[0][0]+p-1)%p;
		}
	}
	return 0;
}