矩阵/整数快速幂/乘

参考来源： 整数快速乘法，整数快速幂，矩阵快速幂【参考博客】

参考题目 poj 3070

先说明一下基本的数学常识：

(a*b) mod c == ( (a mod c) * (b mod c) ) mod c //这最后一个mod c 是为了保证结果不超过c
对于2进制，2n可用1后接n个0来表示、对于8进制，可用公式 i+3j == n （其中 0<= i <=2 ），对于16进制，可用 i+4j==n（0 <= i <=3）来推算，表达形式为2i 后接 j 个0。

接下来让我们尽可能简单的描述快速乘法的思想：

a*b

快速乘法的基本思想 ，是二进制和乘法分配律的结合，（不由得想起浮点数不满足结合律，严重吐槽！！！╮(╯-╰)╭），比如说，13 ==（1101）2 ，413等于4（1101）2 ，用分配律展开得到413 == 4（1000+100+1）2，我们不难观察出，快速幂可以通过判断当前的位（bit）是1还是0，推断出是否需要做求和操作，每次移动到下一位（bit）时，就对ans进行*2操作，等待是否求和。由于除以2和位移操作是等效的，因此这也可以看作是二分思想的应用，这种算法将b进行二分从而减少了不必要的运算，时间复杂度是log（n）。

a^b

快速幂其实可以看作是快速乘法的特例，在快速幂中，我们不再对ans进行*2操作，因为在a^{b中b的意义已经从乘数变成了指数，但是我们可以仍然把b写成二进制，举例说明：此时，我们将4*13改为4}13，13=（1101）2 ，二进制13写开我们得到（1000+100+1），注意，这里的所有二进制是指数，指数的相加意味着底数相乘，因此有4^13 == 48 * 44 * 41。再注意到指数之间的2倍关系，我们就可以用很少的几个变量，完成这一算法。这样，我们就将原本用循环需要O(n)的算法，改进为O（logN）的算法。
整数的快速乘法， 和快速幂。

// 快速乘 
long long qmul(long long a, long long b){
	long long ans = 0;
	while(b){
		if(b & 1){
			ans += a;
		}
		a <<= 1;// kn 
		b >>= 1;
	}
	return ans;
}
// 快速幂 
long long qpow(long long a, long long b){
	long long ans = 1;
	while(b){
		if(b & 1){
			ans *= a;
		}
		a *= a;// n^k 
		b >>= 1;
	}
	return ans;
}

矩阵快速乘法，代码很简单。主要是利用了整数快速幂的思路。
下面是参考题目 poj 3070的AC代码：

#include<stdio.h>
#include<string.h>
// 矩阵快速乘法 poj 3070

struct matrix{
	long long num[2][2]; 
	matrix(){// 初始化 
		memset(num, 0, sizeof(num[0][0]) << 2); 
	} 
};
int N = 2; 
// 矩阵乘法
matrix matrix_mul(const matrix& a,const matrix& b) {
	matrix ans;
 
	for(int i = 0; i < N; i++){
		for(int j = 0; j < N; j++){
			for(int k = 0; k < N; k++){
				ans.num[i][j] += a.num[i][k] * b.num[k][j];
			}
		}
	}
	return ans;
}
// 矩阵乘法 取模 
matrix matrix_mul_mod(const matrix& a,const matrix& b, const long long& mod) {
	matrix ans;
	
	for(int i = 0; i < N; i++){
		for(int j = 0; j < N; j++){
			for(int k = 0; k < N; k++){
				ans.num[i][j] = (ans.num[i][j] % mod + (a.num[i][k] * b.num[k][j]) % mod) % mod;
			}
		}
	}
	return ans;
}
// 矩阵乘法 
matrix matrix_mul(const matrix& a, int aLen,const matrix& b, int bLen) {
	matrix ans;
 
	for(int i = 0; i < aLen; i++){
		for(int j = 0; j < bLen; j++){
			for(int k = 0; k < N; k++){
				ans.num[i][j] += a.num[i][k] * b.num[k][j];
			}
		}
	}
	return ans;
}
// 矩阵的快速幂 
matrix& q_matrix_pow(matrix ma, int len, int n){
 	
	matrix& new_ma = (*(new matrix));
	
	for(int i = 0; i < len; i++){ // 初始化为单位矩阵 
		new_ma.num[i][i] = 1;
	}
	while(n){
		if(n & 1){
			new_ma = matrix_mul(new_ma, ma);
		}
		ma = matrix_mul(ma, ma);
		n >>= 1;
	}
	return (new_ma); 
}
// 矩阵的快速幂 取模 
matrix& q_matrix_pow_mod(matrix ma, const int& len, int n, const long long& mod){
 
	matrix& new_ma = (*(new matrix));
	for(int i = 0; i < len; i++){ // 初始化为单位矩阵 
		new_ma.num[i][i] = 1;
	}
	while(n){
		if(n & 1){
			new_ma = matrix_mul_mod(new_ma, ma, mod);
		}
		ma = matrix_mul_mod(ma, ma, mod);
		n >>= 1;
	}
	return (new_ma); 
}
int main(){
  	int mod = 1e4;
	matrix ma; 
	ma.num[0][0] = 1;
	ma.num[0][1] = 1;
	ma.num[1][0] = 1;
	matrix t;
	matrix f;
	f.num[1][0] = 1;
	
	long long n = 0;
	while(scanf("%lld", &n) != -1 && n != -1) {
		t = q_matrix_pow_mod(ma, 2, n, mod);
		t = matrix_mul(t, 2, f, 1);
		printf("%lld\n", t.num[0][0]);
	}
	
	
	return 0;
}