2018.09.27【BZOJ5221】偏题（矩阵快速幂）（数学推理）

【问题描述】

斐波那契数列是一个经典递推数列，即$$Fi{b}_n=Fi{b}_{n-1}+Fi{b}_{n-2}$$
在这个问题中，定义一个新数列，对于$n\ge 2$有$$F_n=F_{n-1}+F_{n-2}+\sqrt{3+F_{n-1}{F}_{n-2}}$$
对于给出的四个整数$F_0,F_1,M,n$，求$F_{n}modM$。
数据保证$\sqrt{3+F_{n-1}{F}_{n-2}}$总是整数。

【输入格式】

从文件$fib.in$中读入数据。
一行四个整数$F_0,F_1,M,n$。

【输出格式】

输出到文件$fib.out$中。
一行一个整数，表示$F_n$对$M$取模后的值。

【样例输入1】

1 1 10 5

【样例输出1】

4

【样例解释1】

∙ $F_0=1$
∙ $F_1=1$
∙ $F_2=1+1+2\sqrt{3+1\times 1}=6$
∙ $F_3=6+1+2\sqrt{3+6\times 1}=13$
∙ $F_4=13+6+2\sqrt{3+13\times 6}=37$
∙ $F_5=37+13+2\sqrt{3+37\times 13}=94$

【样例输入2】

2 3 100 6

【样例输出2】

82

【数据规模】

对于30%的数据，$n$ ≤ 20。
对于60%的数据，$F_0=F_1=1$。
对于80%的数据，$n\le 10^5$。
对于100%的数据，$0\le n\le 10^9$，$1\le M\le 10^9$，$1\le F_0\le F_1\le 10^6$。

---

解析：

看到以为不是线性递推，整个人就方了。

然而听到有dalao说是二次剩余，貌似并没有$A$掉此题。

思路：

由上面的递推式我们有结论$$\sqrt{3+F_n{F}_{n-1}}=\sqrt{3+F_{n-1}{F}_{n-2}}+F_{n-1}$$

直接构造矩阵，快速幂推一下就行了。

---

代码：

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define re register
#define gc getchar
#define pc putchar
#define cs const

inline
ll getint(){
	re ll num;
	re char c;
	while(!isdigit(c=gc()));num=c^48;
	while(isdigit(c=gc()))num=(num<<1)+(num<<3)+(c^48);
	return num;
}

ll mod;
ll n,F1,F0;

struct matrix{
	ll a[3][3];
	matrix(){memset(a,0,sizeof a);}
	
	void origin(){
		a[0][0]=1,a[0][1]=0,a[0][2]=0;
		a[1][0]=0,a[1][1]=1,a[1][2]=0;
		a[2][0]=0,a[2][1]=0,a[2][2]=1;
	}
	
	void init(){
		a[0][0]=1,a[0][1]=1,a[0][2]=1;
		a[1][0]=1,a[1][1]=0,a[1][2]=0;
		a[2][0]=2,a[2][1]=0,a[2][2]=1;
	}
	
	ll *const operator[](cs int &offset){
		return a[offset];
	}
	
	friend matrix operator*(matrix t,matrix b){
		matrix c;
		(c[0][0]=t[0][0]*b[0][0]%mod+t[0][1]*b[1][0]%mod+t[0][2]*b[2][0]%mod)%=mod;
		(c[0][1]=t[0][0]*b[0][1]%mod+t[0][1]*b[1][1]%mod+t[0][2]*b[2][1]%mod)%=mod;
		(c[0][2]=t[0][0]*b[0][2]%mod+t[0][1]*b[1][2]%mod+t[0][2]*b[2][2]%mod)%=mod;
		(c[1][0]=t[1][0]*b[0][0]%mod+t[1][1]*b[1][0]%mod+t[1][2]*b[2][0]%mod)%=mod;
		(c[1][1]=t[1][0]*b[0][1]%mod+t[1][1]*b[1][1]%mod+t[1][2]*b[2][1]%mod)%=mod;
		(c[1][2]=t[1][0]*b[0][2]%mod+t[1][1]*b[1][2]%mod+t[1][2]*b[2][2]%mod)%=mod;
		(c[2][0]=t[2][0]*b[0][0]%mod+t[2][1]*b[1][0]%mod+t[2][2]*b[2][0]%mod)%=mod;
		(c[2][1]=t[2][0]*b[0][1]%mod+t[2][1]*b[1][1]%mod+t[2][2]*b[2][1]%mod)%=mod;
		(c[2][2]=t[2][0]*b[0][2]%mod+t[2][1]*b[1][2]%mod+t[2][2]*b[2][2]%mod)%=mod;
		return c;
	}
}ss;

inline
matrix quickpow(ll b){
	matrix a;
	a.init();
	matrix ans;
	ans.origin();
	while(b){
		if(b&1){
			ans=ans*a;
		}
		a=a*a;
		b>>=1;
	}
	return ans;
}

signed main(){
	F0=getint();
	F1=getint();
	mod=getint();
	n=getint();
	if(n==0){
		cout<<F0%mod;
		return 0;
	}
	ss[0][0]=F1%mod;
	ss[0][1]=F0%mod;
	ss[0][2]=sqrt(3ll+F1*F0);
	ss[0][2]%=mod;
	ss=ss*quickpow(n-1);
	cout<<ss[0][0]<<endl;
	pc('\n');
	/*循环展开打表代码
	freopen("task.txt","w",stdout);
	for(int re i=0;i<3;++i)
	for(int re j=0;j<3;++j)
	cout<<"(c["<<i<<"]["<<j<<"]=a["<<i<<"][0]*b[0]["<<j<<"]%mod+a["<<i<<"][1]*b[1]["<<j<<"]%mod+a["<<i<<"][2]*b[2]["<<j<<"]%mod)%=mod;"<<endl; 
	*/
	return 0;
}