数论之矩阵快速幂

数论之矩阵快速幂

快速幂是数论中最简单的几种算法之一,还记得我第一个学习的算法就是快速幂、快速乘。

快速幂顾名思义,就是快速算某个数的多少次幂。其时间复杂度为 O(log₂N), 与朴素的O(N)相比效率有了极大的提高。他的基本原理是二进制。(快速乘也使用了二进制)

我们都知道对于任何一个整数N,都能用二进制来表示。那么对于a^n , n一定可以用二进制表示。

比如a^156,而156(10)=10011100(2)

那么

我们就按照这个公式来求解a^156,原来要进行156-1=155次乘法运算,现在的差不多运算次数就是他 二进制的长度*二进制中1的个数=8*4=24次

对于一般的快速幂就不再多说什么了,下面直接给出普通快速幂的模板:

int fastpow(int base,int n,int mod){
	int ans=1;
	while(n){
		if(n&1) ans*=base%mod;
		base*=base;
		n>>=1;
	}
	return ans%mod;
}

相对而言,快速幂实在是简单不过来,所以很少有直接考快速幂的题。而对快速幂的应用大多是在其他比较复杂的算法中充当某个环节的组成,比如后面我们要讲的拓展欧几里得定理、欧拉公式等。


矩阵快速幂

PS:去年一开始接触矩阵快速幂的时候,我还完全不清楚矩阵的运算到底是怎么一回事儿,之前高中也没有接触过线代方面的知识。在刷了数篇博客后总算对矩阵的运算有了一些了解,但对关系矩阵到底这是怎么得出来的,还是一窍不通,所以那时候做题基本靠猜。

终于在这个假期集训的时候才总算搞懂了矩阵快速幂的原理,现在回头一看,卧槽,原来之前令我头皮发麻的题辣么简单!其实算法中很多知识点都一样,没搞懂之前一头雾水,懂了之后就会发现,哦,原来就这么回事儿而已。废话不多说,开始讲矩阵快速幂。



相对于一般的快速幂,矩阵快速幂仅仅是把他的底数和乘数换成了矩阵形式,而相应的乘法运算什么的也遵循矩阵的运算法则。对于矩阵的基本运算这里也不在多做赘述。

矩阵快速幂主要是用于求一个很复杂的递推式的某一项问题。

我觉得矩阵快速幂的难点主要是在关系矩阵的构造上,只要关系矩阵构造出来了,其他的就只是套模板运算而已。

下面我们用一个例题来了解矩阵快速幂的解法。

acm.hdu.edu.cn/showprob

题意:给你一个233矩阵,横纵坐标都是从0开始,第一行为0 233 2333 23333…… 依次往下延续,第一列除(0,0)坐标点之外,让你自由输入n个数,整个矩阵满足关系式:

a(i,j)= a(i-1,j)+a(i,j-1) ( i,j ≠ 0 且 a(i,j)表示矩阵中坐标为(i,j)的点的值)

求坐标为(n,m)的数的值。

这道题乍一看感觉很好解,你看矩阵都给已经给你了,套模板不就行了嘛,so easy!

幼稚!!!

我怎么可能把一道模板题当例题来讲!

这道题的图中的矩阵其实是迷惑你的障眼法,看似是关系矩阵,实质上是给出的递推式!(当初我就被坑的很惨,对着这个矩阵看了半天一点头绪没有,后来又回头看的时候恍然大悟,我想说,出这道题的人真腹黑~)

他的递推关系是这样给出的,对于一列的数我们可以经过某种变换产生下一列的数,而这中间的变换关系就是我们要求的关系矩阵。

我们要求这个矩阵中(n,m)位置的值,其实就是把第一列的n个数经过m次变换得到新的一列数,而答案就是得到这列数的第n个。

递推关系如下:

/*
	
	
	23		10 0 0 0 1				233				
	a1		10 1 0 0 1				a1+23*10+3
	a2	*	10 1 1 0 1	^m	=		a2+a1+23*10+3		^(m-1)
	a3		10 1 1 1 1				a3+a2+a1+23*10+3
	3   	        0  0 0 0 1				3
*/

你可能会说,哎,不对啊,第一列数不应该是,0 a1 a2 a3 吗?

为毛0变成23了,最后面还加了一个3?

其实这是为简化运算。

由条件可知,第一行为0 233 2333……为了简化运算,我们把第一个0改为23,这样一来,第一行就都满足 i=(i-1)*10+3 的关系。

那为什么要在最后加上3呢?还是与 i=(i-1)*10+3 这个表达式有关,为了得到下一行的第一个数,也就是2333,对于( 23,a1,a2,a3 )我们可以构造出( 10, 0, 0, 0 )这一行数据,但下一行还少个3,所以我们在第一行最后加一个3,在构造出的数据行最后加一个1,也就是( 10,0,0,0,1 ),这样一来,(23,a1,a2,a3,3)*(10,0,0,0,1)结果就为233了,中间三个0是因为233这个数与a1,a2,a3无关。

又因为 a(i,j)= a(i-1,j)+a(i,j-1) 所以233下面的数也就可以推出来了,为a1+233,我们可以将其拆成a1+23*10+3。

而这个数可以构造出第二行数据,即(10,1,0,0,1),剩下的数据也可以依次得出。

对于最后3这个数,为了使下一行与当前行拥有相同的性质,所以得出(0,0,0,1)这行数据。

至此,这道题的关系矩阵就已经全部构造出来了,剩下的就是套模板计算了。

矩阵快速幂的模板如下:

struct Matrix{
	ll m[15][15];
	
	Matrix(){
		memset(m,0,sizeof(m));
		for(int i=0;i<=n+1;i++){
			m[i][n+1]=1;
			if(i==n+1) continue;
			m[i][0]=10;
			for(int j=1;j<=i;j++)
				m[i][j]=1;
		}
	}
	
	void clear(){
		memset(m,0,sizeof(m));
		for(int i=0;i<=n+1;i++)
			m[i][i]=1;
	}
	
	void display(){
		cout<<"Matrix's begin:"<<endl;
		for(int i=0;i<=n+1;i++)
			for(int j=0;j<=n+1;j++)
				if(j<n+1) cout<<m[i][j]<<" ";
				else cout<<m[i][j]<<endl;
		cout<<"Matrix's end:"<<endl;
	}
	
	friend Matrix operator*(Matrix a,Matrix b){
		Matrix ans;
		for(int i=0;i<=n+1;i++)
			for(int j=0;j<=n+1;j++){
				ans.m[i][j]=0;
				for(int k=0;k<=n+1;k++)
					ans.m[i][j]+=a.m[i][k]*b.m[k][j]%mod;
			}
		return ans;
	}
	
	friend Matrix operator^(Matrix base,ll k){
		Matrix ans;
		ans.clear();
		while(k){
			if(k&1) ans=ans*base;
			base=base*base;
			k>>=1;
		}
		return ans;
	}
	
};

我全部用了一个结构体实现,这样看起来舒服一点~

AC代码:

#include<iostream>
#include<cstdio>
#include<cstring>
using namespace std;
typedef long long ll;
const int mod=10000007;
const int N=1e9+10;

ll n,M,arr[15],ans;

struct Matrix{
	ll m[15][15];
	
	Matrix(){
		memset(m,0,sizeof(m));
		for(int i=0;i<=n+1;i++){
			m[i][n+1]=1;
			if(i==n+1) continue;
			m[i][0]=10;
			for(int j=1;j<=i;j++)
				m[i][j]=1;
		}
	}
	
	void clear(){
		memset(m,0,sizeof(m));
		for(int i=0;i<=n+1;i++)
			m[i][i]=1;
	}
	
	void display(){
		cout<<"Matrix's begin:"<<endl;
		for(int i=0;i<=n+1;i++)
			for(int j=0;j<=n+1;j++)
				if(j<n+1) cout<<m[i][j]<<" ";
				else cout<<m[i][j]<<endl;
		cout<<"Matrix's end:"<<endl;
	}
	
	friend Matrix operator*(Matrix a,Matrix b){
		Matrix ans;
		for(int i=0;i<=n+1;i++)
			for(int j=0;j<=n+1;j++){
				ans.m[i][j]=0;
				for(int k=0;k<=n+1;k++)
					ans.m[i][j]+=a.m[i][k]*b.m[k][j]%mod;
			}
		return ans;
	}
	
	friend Matrix operator^(Matrix base,ll k){
		Matrix ans;
		ans.clear();
		while(k){
			if(k&1) ans=ans*base;
			base=base*base;
			k>>=1;
		}
		return ans;
	}
	
};

int main(){
//	cout<<"hello"<<endl;
	while(~scanf("%lld%lld",&n,&M)){
		Matrix base;ans=0;
		arr[0]=23;arr[n+1]=3;
		for(int i=1;i<=n;i++){
			scanf("%lld",&arr[i]);
		}
		base=base^M;
//		base.display();
		for(int i=0;i<=n+1;i++){
			ans+=arr[i]*base.m[n][i];
//			cout<<base.m[n][i]<<endl;
		}
		printf("%d\n",ans%mod);
	}
	return 0;
}
/*
	3 1
	1 2 3
	
	23		10 0 0 0 1				233				
	a1		10 1 0 0 1				a1+23*10+3
	a2	*	10 1 1 0 1	^m	=		a2+a1+23*10+3		^(m-1)
	a3		10 1 1 1 1				a3+a2+a1+23*10+3
	3   	        0  0 0 0 1				3
*/

其实刷的题多了,你就会发现,矩阵快速幂的难点在关系矩阵的构造上,而构造关系矩阵也就两个基本套路,要么是递推式、要么是线性变化,只要你理解了矩阵快速幂这个算法的本质,这种类型的题对你来说就so easy了~


应boss要求,加一个应用场景……这大概不算应用吧π_π,还是题……

哪个地方会直接用到矩阵快速幂呢?对,就是图论中求对应两个点长度为n的路径数,这是矩阵乘法在图论中的经典应用,其实就是用到了矩阵乘法的特殊性,因为矩阵乘法有三重循环,最里面的一层循环就是在枚举一个点ki -> k -> j,那么从i -> j的路径数长度为m的条数就等于i-> k路径长度为n的路径数乘上k -> j路径长度为m - n的的路径条数,枚举的过程中全加起来,可以看到,这个过程和矩阵乘法的过程是一样的。

矩阵An中的ai,j表示:图中点i到点j经过n条边的路径数。

例题:acm.hdu.edu.cn/showprob

题解下次补上π_π


到这里矩阵快速幂的学习就结束了,你可能认为仅仅是讲了一个用例,不够全面,也不够通透。但我认为最快的学习方法还是多刷几道题,自己多想想,就像上面说的,有些知识点懂了就是懂了,关键还是要靠自己。

一起努力吧~

怀挺!

编辑于 2018-08-22

文章被以下专栏收录