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 long
using 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;
}