快速幂的思想以及应用

什么是快速幂？
如：3^4 = 3 * 3 * 3 * 3 在计算机中需要执行4次运算，优化算法如下：
3^4 = 3^2 * 3^2 = 9 * 9 比较以上方程，在计算机中只需要执行2次运算，足足少了一半的计算量。当次方越大时，优化效果最为明显。
上面例子为次方为偶数情况，当为奇数时，只需要先乘一个当前底数在继续缩半乘方就可。
如：3^5 = 3 * 3^2 * 3^2 = 3 * 9 * 9

由上面例子可知
快速幂思想：当需要 a^b 的值时， 判断 b是否为奇数，若为 ans = ans *a,b = b-1，若为偶数 a = a * a,同时b = b/2;重复上述，直到b为0，ans即为所求

快速幂代码如下：

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;

ll ksm(ll a,ll b)
{
	ll ans = 1;
	while(b)
	{
		if(b & 1) ans = ans * a;
		a *= a;
		b>>=1;
	}
	return ans;
}

我们来探究一下矩阵多次方的情况
如有矩阵 A ^5 = A * A * A * A * A，因为矩阵相乘运算满足结合律则
A^5 = A * (A * A) * (A * A) = A * A^2 * A^2，是否与数的多次幂的快速幂思想雷同？。
所以快速幂思想也可用在矩阵多次幂的运算下
矩阵快速幂代码如下：

#include<iostream>
#include<cstdio>
#include<vector>

using namespace std;
vector<vector<int>> mul_matrix(const vector<vector<int>> &a,const vector<vector<int>>& b)
{
    int ai = a.size(),aj = a[0].size();
    int bi = b.size(),bj = b[0].size();
    vector<vector<int>>c(ai,vector<int>(bj,0));
    if(aj != bi) return c;
    
    int ci = c.size(), cj = c[0].size();
    
    for(int k = 0; k < ci; k++){
        for(int i = 0 ; i < cj; i++){
           
            for(int j = 0; j <bi; j++){
                 c[k][i] = c[k][i] + a[k][j] * b[j][i];
            }
        }
    }
    
    return c;
}


vector<vector<int>> ksm_matrix(vector<vector<int>> a,int n)
{
    int ai = a.size(),aj = a[0].size();
    vector<vector<int>> c(ai,vector<int>(aj,0)); //单位矩阵E ， E*A = A
    for(int i = 0 ; i < ai; i++) c[i][i] = 1;
    
    while(n)
    {
        if(n&1) c =  mul_matrix(c,a);
        n>>=1;
        a = mul_matrix(a,a);
    }
    return c;
}

应用题目：
poj 1995

poj 3070

#include<iostream>
#include<cstdio>
#include<vector>

using namespace std;

const int MOD = 10000;
vector<vector<int>>a(2,vector<int>(2,0));
vector<vector<int>>f(2,vector<int>(1,0));




vector<vector<int>> mul_matrix(const vector<vector<int>> &a,const vector<vector<int>>& b)
{
    int ai = a.size(),aj = a[0].size();
    int bi = b.size(),bj = b[0].size();
    vector<vector<int>>c(ai,vector<int>(bj,0));
    if(aj != bi) return c;
    
    int ci = c.size(), cj = c[0].size();
    
    for(int k = 0; k < ci; k++){
        for(int i = 0 ; i < cj; i++){
           
            for(int j = 0; j <bi; j++){
                 c[k][i] = (c[k][i] + a[k][j] * b[j][i] %MOD)%MOD;
            }
        }
    }
    
    return c;
}


vector<vector<int>> ksm_matrix(vector<vector<int>> a,int n)
{
    vector<vector<int>> c(2,vector<int>(2,0));
    c[0][0] = c[1][1] = 1;
    
    while(n)
    {
        if(n&1) c =  mul_matrix(c,a);
        n>>=1;
        a = mul_matrix(a,a);
    }
    return c;
}

int main()
{
    while(true)
    {
        int n;
        cin >> n;
        if(n==-1) break;
        
        a[0][0] = 0;a[0][1] = a[1][0] = a[1][1] = 1;
        f[0][0] = 0;f[1][0] = 1;
    
        vector<vector<int>> c = ksm_matrix(a,n);
        c = mul_matrix(c,f);
        cout << c[0][0]<<endl;        
    }

}

poj 3734 dp + 矩阵快速幂

/*
    用一个三维数组来表示红黄蓝绿的格子个数
    一维表示绿，二维表示红，三维表示前i个
    
    f[i][0][0] 表示 前i个格子中，红，绿为偶数
    f[i][0][1] 表示 前i个格子中，红为偶数，绿为奇数
    f[i][1][0] 表示 前i个格子中，红为奇数，绿为偶数
    f[i][1][1] 表示 前i个格子中，红，绿均为奇数
    
    在f[i][0][0]这种情况下：
    当第i个格子为红色时 前i-1个格子有 f[i-1][1][0]
    当第i个格子为绿色时，前i-1个格子有 f[i-1][0][1]
    当第i个格子为黄或蓝，前i-1个格子有 f[i-1][0][0]*2 
    故f[i][0][0] = f[i-1][1][0] + f[i-1][0][1] + f[i-1][0][0] * 2
    
    在f[i][0][1]这种情况下:
    当第i个格子为红色时，前i-1个格子有f[i-1][1][1]
    当第i个格子为绿色时,前i-1个格子有 f[i-1][0][0]
    当第i个格子为黄或蓝时，前i-1个格子有f[i-1][0][1]*2
    故f[i][0][1] = f[i-1][1][1] + f[i-1][0][0] + f[i-1][0][1] *2
    
    在f[i][1][0]这种情况下:
    当第i个格子为红色时，前i-1个格子有f[i-1][0][0]
    当第i个格子为绿色时，前i-1个格子有f[i-1][1][1]
    当第i个格子为黄或者蓝时，前i-1个格子有f[i-1][1][0]*2
    故f[i][1][0] = f[i-1][0][0] + f[i-1][1][1] + f[i-1][1][0] * 2
    
    在f[i][1][1]这种情况下：
    当第i个格子为红时，前i-1个格子有f[i-1][0][1]
    当第i个格子为绿时，前i-1个格子有f[i-1][1][0]
    当第i个格子为黄或者蓝时，前i-1个格子有f[i-1][1][1]*2
    故f[i][1][1] = f[i-1][0][1] + f[i-1][1][0] + f[i-1][1][1]*2
      
    边界条件
    假设有m个格子
    i<=m中
    f[0][0][0] = 2;
    f[0][0][1] = f[0][1][0] = 1;
    f[0][1][1] = 0;
    
    
    f[i][0][0] =2*f[i-1][0][0] + 1*f[i-1][0][1]  + 1*f[i-1][1][0] + 0*f[i-1][1][1]
    f[i][0][1] =1*f[i-1][0][0] + 2*f[i-1][0][1]  + 0*f[i-1][1][0] + 1*f[i-1][1][1] 
    f[i][1][0] =1*f[i-1][0][0] + 0*f[i-1][0][1]  + 2*f[i-1][1][0] + 1*f[i-1][1][1] 
    f[i][1][1] =0*f[i-1][0][0] + 1*f[i-1][0][1]  + 1*f[i-1][1][0] + 2*f[i-1][1][1]
    
    
*/
#include<iostream>
#include<cstdio>
#include<vector>

using namespace std;

const int MOD = 10007;

vector<vector<int>> mul_matrix(const vector<vector<int>> &a,const vector<vector<int>>& b)
{
    int ai = a.size(),aj = a[0].size();
    int bi = b.size(),bj = b[0].size();
    vector<vector<int>>c(ai,vector<int>(bj,0));
    if(aj != bi) return c;
    
    int ci = c.size(), cj = c[0].size();
    
    for(int k = 0; k < ci; k++){
        for(int i = 0 ; i < cj; i++){
           
            for(int j = 0; j <bi; j++){
                 c[k][i] = (c[k][i] + a[k][j] * b[j][i] %MOD)%MOD;
            }
        }
    }
    
    return c;
}


vector<vector<int>> ksm_matrix(vector<vector<int>> a,int n)
{
    int ai = a.size(),aj = a[0].size();
    vector<vector<int>> c(ai,vector<int>(aj,0));
    for(int i = 0 ; i < ai; i++) c[i][i] = 1;
    
    while(n)
    {
        if(n&1) c =  mul_matrix(c,a);
        n>>=1;
        a = mul_matrix(a,a);
    }
    return c;
}

void show(const vector<vector<int>> & a)
{
    for(int i = 0 ; i <(int) a.size(); i++){
        for(int j = 0 ; j <(int) a[i].size(); j++){
            cout << a[i][j]<<" ";
        }
        cout << endl;
    }
}

int main()
{
    
    /*
        f[0][0][0] = 2;
        f[0][0][1] = f[0][1][0] = 1;
        f[0][1][1] = 0;
    */
    
    vector<vector<int>>matrix(4,vector<int>(4,0));
    for(int i = 0 ; i < 4 ; i++) matrix[i][i] = 2;
    matrix[0][1] = matrix[1][0] = 
    matrix[0][2] = matrix[2][0] = 
    matrix[1][3] = matrix[3][1] = 
    matrix[2][3] = matrix[3][2] = 1;
    vector<vector<int>>temp(1,vector<int>(4,0));
    temp[0][0] = 2;
    temp[0][1] = temp[0][2] = 1;
    temp[0][3] = 0;
    
    //show(matrix);
    int op = 0;
    scanf("%d",&op);
    
    while(op--)
    {
        int n;
        scanf("%d",&n);
        vector<vector<int>> ans = ksm_matrix(matrix,n-1);
        ans = mul_matrix(temp,ans);
        cout << ans[0][0]<<endl;
    }
    
    
    return 0;
}