题意：给a, b, n, m
求 $left lceil ( a+ sqrt b )^n ight ceil$ % m
看到 $( a+ sqrt b )^n$ 虽然很好联想到共轭 但是推出矩阵还是比较难的
$(a+sqrt b)^n + (a-sqrt b)^n$
= $(C^0_n a^n + C^1_n a^{n-1} sqrt{b} + ... + C^{n-1}_n a sqrt{b}^{n-1} + C^n_n (sqrt b)^n)$
+ $(C^0_n a^n + (-1)^1 C^1_n a^{n-1} sqrt{b} + ... + (-1)^{n-1} C^{n-1}_n a sqrt{b}^{n-1} + (-1)^n C^n_n sqrt b^{ n})$
对于 n+1 项 $C^i_n$ 其中所有 i 为奇数的都消掉了 只剩下偶数项
= $2 sumlimits_{i=0,i为偶数}^n C^i_n a {sqrt b}^{ i}$
偶数 所以可以把根号开掉 得到
= $2 sumlimits_{i=0}^{frac{n}{2}} C^i_n a {b}^{ i}$
很明显这是一个整数
数据范围为：$(a-1)^2 < b < a^2$
那么 $a-1 < sqrt b < a$
那么 $0 < a-sqrt b < 1$
$(a+sqrt b)^n + (a-sqrt b)^n$ 是个整数 同时 $a-sqrt b$ 大于0 小于1
因此 $left lceil ( a+ sqrt b )^n ight ceil = (a+sqrt b)^n + (a-sqrt b)^n$ (画上这个等号真是不容易!)
求 $(a+sqrt b)^n + (a-sqrt b)^n$ 就变得很简单了
通过二次特征方程$x^2-2ax+(a^2-b)=0$
可以得到$S_n = 2aS_{n-1}+(b-a^2)S_{n-2}$
写成矩阵形式：
$$egin{pmatrix} S_n\S_{n-1}end{pmatrix} = egin{pmatrix} a & {b-a^2}\ 1 & 0end{pmatrix} egin{pmatrix} S_{n-1}\S_{n-2}end{pmatrix}$$
$$= {egin{pmatrix} a & {b-a^2}\ 1 & 0end{pmatrix}}^{n-1} egin{pmatrix} 2a \ 2 end{pmatrix}$$
$b-a^2$ 很容易为负 n-1次方之后 更加负
为了防止答案为负 记得多加几个mod哟~
1 #include <bits/stdc++.h> 2 using namespace std; 3 typedef long long LL; 4 LL mod; 5 struct Mat 6 { 7 LL t[2][2]; 8 }; 9 Mat mul(Mat a, Mat b) 10 { 11 Mat c; 12 memset(c.t, 0, sizeof(c.t)); 13 for(int i=0;i<2;i++) 14 for(int k=0;k<2;k++) 15 if(a.t[i][k]) 16 for(int j=0;j<2;j++) 17 c.t[i][j]=(c.t[i][j]+a.t[i][k]*b.t[k][j])%mod; 18 return c; 19 } 20 Mat expo(Mat p, LL k) 21 { 22 if(k==1) 23 return p; 24 Mat e; 25 memset(e.t, 0, sizeof(e.t)); 26 for(int i=0;i<2;i++) 27 e.t[i][i]=1; 28 if(k==0) 29 return e; 30 while(k) 31 { 32 if(k & 1) 33 e=mul(p, e); 34 p=mul(p, p); 35 k>>=1; 36 } 37 return e; 38 } 39 LL MOD(LL x) 40 { 41 while(x<0) 42 x+=mod; 43 return x%mod; 44 } 45 int main() 46 { 47 //freopen("in.txt", "r", stdin); 48 //freopen("out.txt", "w", stdout); 49 LL a, b, n; 50 while(cin>>a>>b>>n>>mod) 51 { 52 Mat m; 53 m.t[0][0]=2*a, m.t[0][1]=b-a*a; 54 m.t[1][0]=1, m.t[1][1]=0; 55 Mat ans=expo(m, n-1); 56 cout<<MOD(ans.t[0][0]*2*a+2*ans.t[0][1])<<endl; 57 } 58 return 0; 59 }