【bzoj4002】[JLOI2015]有意义的字符串【数论】【矩阵快速幂】

题意：求$\lfloor (\frac{b+\sqrt{d}}{2}{)}^n\rfloor mod7528443412579576937$的值。
题解：get到了好多的新姿势qwq。
首先get到了特征多项式相关姿势。orz zjr and yeh dark♂佬，直接秒了特征多项式。
我们设有一个数列$x$，满足$x_{n+2}=c_1{x}_{n+1}+c_2{x}_n$。
我们设存在$r$，$s$，满足$x_{n+2}-r{x}_{n+1}=s(x_{n+1}-r{x}_n)$。
移项得$x_{n+2}=(s+r)x_{n+1}-sr\times x_n$
可以得到$c_1=s+r(1)$，$c_2=-sr(2)$。
我们让$(1)\times r+(2)$，得到$r{c}_1+c_2=r^2$
这就是一个一元二次方程，可以直接解出。
我们设这个方程的两个根为$p,q$。
存在$A,B$使得$x_n=A\times p^n+B\times q^n$。
只需要带入两个$n$的值，就可以把$A$和$B$求出，也就是得出了通项公式。
我们现在要求的就是$(\frac{b+\sqrt{d}}{2}{)}^n$。
我们构造一个式子$((\frac{b+\sqrt{d}}{2}{)}^n+(\frac{b-\sqrt{d}}{2}{)}^n)-(\frac{b-\sqrt{d}}{2}{)}^n$
前面这些东西$((\frac{b+\sqrt{d}}{2}{)}^n+(\frac{b-\sqrt{d}}{2}{)}^n)$，我们可以把这个看成一个递推式的通项，得到$A=1,B=1,p=\frac{b+\sqrt{d}}{2},q=\frac{b-\sqrt{d}}{2}$。所以现在我们要做的就是由通项求出递推式。
因为$p,q$是方程$r{c}_1+c_2=r^2$的两个根，我们可以尝试把他们带入，得到
$$(\frac{b+\sqrt{d}}{2}{)}^2=c_1\times (\frac{b+\sqrt{d}}{2})+c_2\text{\hspace{1em}(3)}$$
$$(\frac{b-\sqrt{d}}{2}{)}^2=c_1\times (\frac{b-\sqrt{d}}{2})+c_2\text{\hspace{1em}(4)}$$

化简得
$$\frac{b^2+2b\sqrt{d}+d}{4}=\frac{c_{1}b+c_1\sqrt{d}}{2}+c_2\text{\hspace{1em}(5)}$$
$$\frac{b^2-2b\sqrt{d}+d}{4}=\frac{c_{1}b-c_1\sqrt{d}}{2}+c_2\text{\hspace{1em}(6)}$$

$(3)+(4)$得
$$\frac{b^2+d}{2}=c_{1}b+2c_2\text{\hspace{1em}(7)}$$

$(3)-(4)$得$c_1=b$。带入$(7)$得到$c_2=\frac{d-b^2}{4}$。
于是我们成功推出了递推式：
$x_n=b\times x_{n-1}+\frac{d-b^2}{4}\times x_{n-2}$。
很容易得到边界$x_0=2,x_1=b$。
注意到$b\equiv 1mod2$，$d\equiv 1mod4$，所以$(\frac{b+\sqrt{d}}{2}{)}^n+(\frac{b-\sqrt{d}}{2}{)}^n$一定是整数。所以我们可以直接跑整数的矩阵乘法。
注意到$b^2\le d<(b+1{)}^2$，所以$(\frac{b-\sqrt{d}}{2}{)}^n$一定在$(-1,0]$内。$(\frac{b-\sqrt{d}}{2}{)}^n$对答案有影响，当且仅当$b\times b\ne d$且$n$为偶数。这时答案要减1。
注意模数很恶心，一加都爆$longlong$，不仅要用$O(1)$快速乘，还要用$unsignedlonglong$。
代码：

#include<cstdio>
#include<cstring>
#include<cmath>
#define int unsigned long long
const int mod=7528443412579576937;
int b,d,n,k,a[2][2],base[2][2],tmp[2][2];
int mul(int a,int b){
	long long c=a*(long double)b/mod+0.5;
	long long res=a*b-c*mod;
	if(res<0){
		res+=mod;
	}
	return res;
}
int fastpow(int a,int x){
	int res=1;
	while(x){
		if(x&1){
			res=mul(res,a);
		}
		x>>=1;
		a=mul(a,a);
	}
	return res;
}
void mul(int a[][2],int b[][2],int c[][2]){
	memset(tmp,0,sizeof(tmp));
	for(int i=0;i<2;i++){
		for(int j=0;j<2;j++){
			for(int k=0;k<2;k++){
				tmp[i][j]=(tmp[i][j]+mul(a[i][k],b[k][j]))%mod;
			}
		}
	}
	memcpy(c,tmp,sizeof(tmp));
}
signed main(){
	scanf("%lld%lld%lld",&b,&d,&n);
	if(!n){
		puts("1");
		return 0;
	}else if(n==1){
		printf("%lld\n",(int)(0.5*(b+sqrt(d))));
		return 0;
	}
	a[0][0]=2;
	a[0][1]=b;
	base[0][1]=(d-b*b)/4;
	base[1][0]=1;
	base[1][1]=b;
	int k=n;
	while(k){
		if(k&1){
			mul(a,base,a);
		}
		k>>=1;
		mul(base,base,base);
	}
	if(n%2==0&&b*b!=d){
		printf("%lld\n",a[0][0]-1);
	}else{
		printf("%lld\n",a[0][0]);
	}
	return 0;
}