HDU2256&&HDU4565：给一个式子的求第n项的矩阵快速幂

解决的问题是：

HDU2256

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=2256

题意：求（sqrt(2)+sqrt(3)）^2n%1024是多少。

这个题算是hdu4565的一个常数版本了，所以我们先说这道题。对于这道题的做法我们可以计算((sqrt(2)+sqrt(3))^2)^n=(5+2*sqrt(6))^n，对于(5+2*sqrt(6))^n我们知道答案必定是以an+bn*sqrt(6)，而对于下一项我们只需要求(an+bn*sqrt(6))*(5+2*sqrt(6))=5*an+12*bn+2*an*sqrt(6)+5*bn*sqrt(6)，所以a(n+1)=5*an+12*bn; b(n+1)=2*an+5*bn。有了这个递推式我们就可以构造矩阵求an，bn。

这里还有一点对于(5+2*sqrt(6))^n=an+bn*sqrt(6); 同理(5-2*sqrt(6))^n=an-bn*sqrt(6)；两式相加(5+2*sqrt(6))^n+(5-2*sqrt(6))^n=2*an，当n趋于无穷的时候lim(5-2*sqrt(6))^n=0，因为5-2*sqrt(6)<1。

所以我们可以得到答案(5+2*sqrt(6))^n约等于2*an，且实际值是比2*an要小的且小于2*an-1要大的，所以由题目的意思我们向下取整，ans=2*an-1；具体看代码，其他都是矩阵快速幂的模板。

代码如下：

#include <bits/stdc++.h>using namespace std;const int MOD = 1024;struct Matrix{int n,m;int mat[3][3];void input(int n, int m){for(int i = 1;i <= n;i++){for(int j = 1;j <= m;j++) {scanf("%d",&mat[i][j]);}}}Matrix operator * (const Matrix &b) {Matrix ans;memset(ans.mat,0,sizeof(ans.mat));ans.n = n;ans.m = b.m;for(int i = 1;i <= n;i++){for(int j = 1;j <= b.m;j++){for(int k = 1;k <= m;k++){ans.mat[i][j] += mat[i][k] * b.mat[k][j];ans.mat[i][j] %= MOD;}}}return ans;}}T,p,One,q;void init(){One.n = One.m = 2;for(int i = 1;i <= One.n;i++){One.mat[i][i] = 1;}}Matrix q_pow(Matrix x,int k){Matrix res = One;if(k == 0) return One;while(k > 0){if(k & 1) res = res * x;x = x * x;k >>= 1;}return res;}int n;int main(){int kase;scanf("%d",&kase);init();p.n = 1;p.m = 2;p.mat[1][1] = 5;p.mat[1][2] = 2;while(kase--){scanf("%d",&n);T.n = T.m = 2;T.mat[1][1] = 5;T.mat[1][2] = 2;T.mat[2][1] = 12;T.mat[2][2] = 5;q = p * q_pow(T,n- 1);printf("%d\n",((2 * q.mat[1][1]) - 1) % MOD);}return 0;}

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4565

题意：

算是上面那道题的一个升级版本啦。现在是a和b不是固定的的常数了。和上面的做法一样。注意题目中给出a-sqrt(b)<1的条件，所以基本和上道题是一样的了，类比一下吧！很简单的。但是这里是向上取整，所以答案是2*an。具体看代码吧。由于上面那道题在常数情况下已经说得很明白了。这道题就不说了

代码如下：

#include <bits/stdc++.h>#define ll long longusing namespace std;ll MOD = 0;;struct Matrix{ll n,m;ll mat[3][3];void input(ll n, ll m){for(ll i = 1;i <= n;i++){for(ll j = 1;j <= m;j++) {scanf("%d",&mat[i][j]);}}}Matrix operator * (const Matrix &b) {Matrix ans;memset(ans.mat,0,sizeof(ans.mat));ans.n = n;ans.m = b.m;for(ll i = 1;i <= n;i++){for(ll j = 1;j <= b.m;j++){for(ll k = 1;k <= m;k++){ans.mat[i][j] += mat[i][k] * b.mat[k][j];ans.mat[i][j] %= MOD;}}}return ans;}}T,p,One,q;void init(){One.n = One.m = 2;for(ll i = 1;i <= One.n;i++){One.mat[i][i] = 1;}}Matrix q_pow(Matrix x,ll k){Matrix res = One;if(k == 0) return One;while(k > 0){if(k & 1) res = res * x;x = x * x;k >>= 1;}return res;}ll n,a,b;int main(){init();p.n = 1;p.m = 2;while(~scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&MOD)){p.mat[1][1] = a;p.mat[1][2] = 1;T.n = T.m = 2;T.mat[1][1] = a;T.mat[1][2] = 1;T.mat[2][1] = b;T.mat[2][2] = a;q = p * q_pow(T,n- 1);printf("%I64d\n",((2 * q.mat[1][1])) % MOD);}return 0;}