浅谈线性递推与矩阵加速

这个号是从同学那借来的，此后用这个号发些博客，这是我的第一篇文章。本人刚刚大一，写的很菜见谅，多多指教。

我的计算机科学概论娄老师发了个选做作业，正好以前有些印象，就勉强写写，

斐波那契数列指的是这样一个数列：

前两项为1，此后每一项为前两项之和。递推公式和通项公式为：

通项公式我也不会推导，据说是线性代数和母函数都可以推，详见百度百科，推导也不是我要讲的重点。先说说这个问题怎么用计算机语言解决，这里采用的是C++、C、Java和Python。觉得简单的跳过这一段：

#include
#define N 100001
#define int unsigned long long
#define RE register int
using namespace std;
int n,f[N]={0,1,1};
main(){
	cin>>n;
	for (RE i=3;i<=n;i++)
		f[i]=f[i-1]+f[i-2];
	cout< 
#include
#define N 100001
int n,f[N]={0,1,1};
int main(int argc,char*argv[]){
	int i;
	scanf("%d",&n);
	for (i=3;i<=n;i++)
		f[i]=f[i-1]+f[i-2];
	printf("%d",f[n]);
	return 0;
}
 
import java.util.Scanner;
public class Helloworld{
	public static void main(String[] args){
		Scanner s=new Scanner(System.in);
		int []f=new int [1001];
		f[1]=f[2]=1;
		int n=s.nextInt();
		for (int i=3;i<=n;i++)
			f[i]=f[i-1]+f[i-2];
		System.out.println(f[n]);
	}
}
 
f=[0]*1000
f[1]=f[2]=1
n=int(input())
for i in range(3,n+1):
    f[i]=f[i-1]+f[i-2]
print(f[n])
 
显然算法的复杂度是O(n)，这里有递归的写法，但递归比递推慢些我就不写了，其实是一样的。在一类离线分治算法（陈丹琦分治）中递归则表现出了其优越性。现在我步入重点，所谓的矩阵加速。
 
 
 
 
 嗯么么，真有意思。所以我们得到了
 
 
所以我们把这个O(n)的线性递推问题转化为了一个指数级算法问题。有人可能要问我，这两个有什么区别嘛，知道的也别剧透，下面我给大家说说快速幂，the theorem of quickpower。
 
对于a的b次方的一类问题中，只要b是整数，就可以把b转换成二进制来解决。复杂度从O(n)暴跌到O(logn)（以2为底），这个理解思路是李煜东提出的。具体如下
 
譬如2^101的101二进制为1100101，这个问题就可以转换成：
 
 
我们记录答案为answer=1，将上面的乘法从右往左做：
 
answer不断平方（想想为什么），当二进制下该位为1时，ans乘a即可，最后得到答案。有个技巧是二进制的左右移，这里我将会另开一章探讨基于二进制基的状态压缩规划。然后是代码
 
#include
using namespace std;
#define MOD 10000007
unsigned long long ans=1,a,n;
main(){
	cin>>a>>n;
	while(n){
		if (n&1)	ans*=a;
		ans%=MOD;
		a*=a;
		n>>=1;
	}cout<
 
这里防止溢出取模一个100000007，维护一个素域，Python不用担心爆long long。快速幂也可以用递归写，这样子更好理解噢，有空会补上的。
 
这里的矩阵乘法也满足结合律，快速幂也适用了，所以就可以用log2(n)的速度算了，比如我要算1024位，我只要算10次，实在6。
 
矩阵的加减、乘、求逆可以用C++的重载，代码暂时不放。然后就是矩阵加速的代码，我今天累了，以后再补吧。引用代码不用注明出处，你喜欢也可以。
 
#include
using namespace std;
#define RE register int
#define IL inline
#define N 1001
#define MOD 100000007
unsigned long long n;
struct matrix{
	long long a[3][3];
}A,B;
matrix operator*(matrix&S,matrix&T){
	matrix C;
	memset(C.a,0,sizeof(C.a));
	for (RE i=1;i<=3;i++)
		for (RE j=1;j<=3;j++)
			for (RE k=1;k<=3;k++)
				C.a[i][j]+=S.a[i][k]*T.a[k][j];
	return C;
}
int main(){
	cin>>n;
	A.a[1][1]=A.a[1][2]=1;
	B.a[1][2]=B.a[2][1]=B.a[2][2]=1;
	while(n){
		if (n&1)	A=A*B;
		n>>=1,B=B*B;
		cout<
 
23点熄灯，写到一半，下次再补。