矩阵快速幂

快速幂

思想

我们要让计算机很快地计算出$a^b$。
可以利用乘方的性质：$a^xa^y=a^{x+y}$。
我们要求$a^b$，可以将$b$拆成二进制数。比如$b=(11)_{10}=(1011)_2$，从左至右，每一个$1$，分别代表$8,2,1$，可以计算$a^{11}=a^8a^2a^1$ 。
可以得出，如果$b$在二进制上的第$n$位是$1$,我们就把答案乘上对应的$a^{2^n}$。
在快速幂中，我们通过$a$的不断自乘，得到$a^1,a^2,a^4...$，通过$b$的右移(>>)，自右向左取得$b$的二进制位，用按位与运算(&)判断每一位是否为$1$。

代码

#include<iostream>
using namespace std;
typedef long long ll;
ll quick_pow(ll a,ll b)
{
	ll res=1;
	while(b)
	{
		if(b&1) res*=a;
		b>>=1;
		a=a*a;
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	ll a,b;
	cin>>a>>b;
	cout<<quick_pow(a,b);
	return 0;
}

当然，一般会加上取余运算，计算$a^b mod $ $p$。

#include<iostream>
using namespace std;
typedef long long ll;
const ll p=1e9+7;
ll quick_pow_mod(ll a,ll b)
{
	ll res=1;
	while(b)
	{
		if(b&1)
		{
			res*=a;
			res%=p;
		}
		b>>=1;
		a=a*a;
		a%=p;
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	ll a,b;
	cin>>a>>b;
	cout<<quick_pow_mod(a,b);
	return 0;
}

矩阵快速幂

含义

设$A$是一个$n$阶方阵，我们需要快速地计算$A^b$。

类比

矩阵快速幂与数的快速幂类似，只需要将快速幂代码的数，换成矩阵即可。

实现

①利用结构体定义矩阵，成员为一个$n×n$的二维数组。
②重载‘*’，实现矩阵的乘法${(AB)_{ij}} = \sum\limits_{k = 1}^n {{a_{ik}}{b_{kj}}} $。
③套快速幂的模板。

代码实现

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=100;
ll n,b;
struct matrix
{
	ll m[N][N];
	matrix()//构造函数 
	{
		memset(m,0,sizeof(m));
	}
};
matrix operator *(const matrix &x,const matrix &y)
{
	matrix res;
	int i,j,k;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			for(k=1;k<=n;k++)
			{
				res.m[i][j]+=x.m[i][k]*y.m[k][j];
			}
		}
	}
	return res;
}
matrix quick_matrix_pow_mod(matrix a,ll b)
{
	matrix res;
	int i;
	for(i=1;i<=n;i++) res.m[i][i]=1;//单位矩阵 
	while(b)
	{
		if(b&1) res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	cin>>n>>b;
	int i,j;
	matrix a;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			cin>>a.m[i][j];
		}
	}
	matrix ans=quick_matrix_pow_mod(a,b);
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++) 
		{
			cout<<ans.m[i][j];
			if(j<n) cout<<' ';
		}
		cout<<endl;
	}
	return 0;
}

当然，一般会要求对矩阵的每一个元素取余，只需要在重载‘*’时对元素取余即可。

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=100;
const ll p=1e9+7;
ll n,b;
struct matrix
{
	ll m[N][N];
	matrix()//构造函数 
	{
		memset(m,0,sizeof(m));
	}
};
matrix operator *(const matrix &x,const matrix &y)
{
	matrix res;
	int i,j,k;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			for(k=1;k<=n;k++)
			{
				res.m[i][j]+=x.m[i][k]*y.m[k][j];
				res.m[i][j]%=p;
			}
		}
	}
	return res;
}
matrix quick_matrix_pow_mod(matrix a,ll b)
{
	matrix res;
	int i;
	for(i=1;i<=n;i++) res.m[i][i]=1;//单位矩阵 
	while(b)
	{
		if(b&1) res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	cin>>n>>b;
	int i,j;
	matrix a;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			cin>>a.m[i][j];
		}
	}
	matrix ans=quick_matrix_pow_mod(a,b);
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++) 
		{
			cout<<ans.m[i][j];
			if(j<n) cout<<' ';
		}
		cout<<endl;
	}
	return 0;
}

矩阵加速数列递推 (e.g. Fibonacci Sequence)

矩阵的递推思想

设$Fibonacci$ $Sequence$的第$n$项为$f_n$,$f_0=0,f_1=1$，那么有

$$f_n =f_{n - 1} ×1+ f_{n - 2}×1, f_{n-1}=f_{n-1}×1+f_{n-2}×0(n \geqslant 2)$$

可以写成矩阵的乘法

令,于是可得
而$A^n$可由矩阵快速幂得到, 最后利用矩阵乘法定义得到$f_n=A^n_{21}$，于是就这样在$O(logn)$的时间实现了数列递推。

代码

#include<iostream>
#include<cstring>
using namespace std;
typedef long long ll;
const int N=100;
const ll p=1e9+7;
const int n=2;
struct matrix
{
	ll m[N][N];
	matrix()//构造函数 
	{
		memset(m,0,sizeof(m));
	}
};
matrix operator *(const matrix &x,const matrix &y)
{
	matrix res;
	int i,j,k;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
		{
			for(k=1;k<=n;k++)
			{
				res.m[i][j]+=x.m[i][k]*y.m[k][j];
				res.m[i][j]%=p;
			}
		}
	}
	return res;
}
matrix quick_matrix_pow(matrix a,ll b)
{
	matrix res;
	int i;
	for(i=1;i<=n;i++) res.m[i][i]=1;//单位矩阵 
	while(b)
	{
		if(b&1) res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
int main()
{
	ios::sync_with_stdio(false);
	int k;
	cin>>k;//输出第k项的fib
	int i,j;
	matrix a;
	a.m[1][1]=a.m[1][2]=a.m[2][1]=1;
	a.m[2][2]=0;
	matrix ans=quick_matrix_pow(a,k);
	cout<<ans.m[2][1]%p;
	return 0;
}