矩阵/整数快速幂/乘

朱雀 2022-03-10 03:44 296阅读 0赞

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

[参考题目 poj 3070][poj 3070]

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

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

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

## a\*b ##

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

## a^b ##

> 快速幂其实可以看作是快速乘法的特例，在快速幂中，我们不再对ans进行\*2操作，因为在ab中b的意义已经从乘数变成了指数，但是我们可以仍然把b写成二进制，举例说明：此时，我们将4\*13改为413，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][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;
    }

[Link 1]: https://www.cnblogs.com/luruiyuan/p/5570756.html
[poj 3070]: http://poj.org/problem?id=3070