矩阵经典操作（2）

关于矩阵经典操作（2）只要是关于矩阵求幂或者求前K次方的和

这些题目一般题目会给定矩阵，这属于模板题，比较简单

关于整数快速幂，在数论里有一个比较快速的算法，矩阵的快速幂和这个差不多，有异曲同工的意思，都是二分的思想

整数快速幂：

long long power(long long a,long long b)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
        {
            ans=ans*a%m;
        }
        b>>=1;
        a=a*a%m;
    }
    return ans;
}

这是针对结果%mod的情况，对于不需要%的数，直接输出就可以了，根据情况来选择

这种算法的思想就是A^K=A^(K>>1)*A^(K>>1)或者A^K=A^(K>>1)*A^(K>>1）*A

由于乘法具有结合律，因此A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2 。我们可以得A^(n/2) * A  （其中n/2 取整）有这样的结论：当n为偶数时，A^n = A^(n/2) * A^(n/2) ；当n 为奇数时，A^n = A^(n/2) * 。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3 的值即可。我们可以在计算过程中不断取模，避免高精度运算

关于矩阵的快速幂，思想差不多，只是在初始化的时候ans要注意变成单位矩阵

矩阵快速幂{

matrix power(matrix a,int b)
{
    matrix ans;
    for(int i=0; i<n; i++)//构造单位矩阵，n为矩阵大小，而且单位矩阵一定为n*n
    {
        for(int j=0; j<n; j++)
        {
            ans.a[i][j]=(i==j);
        }
    }

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

在需要%mod时，可以重载*操作，在乘法的过程中%mod

关于 重载*

matrix operator +(const matrix &b)const
    {
        matrix c;
        c.clear();
        for(int i=0; i<2; i++)
        {
            for(int j=0; j<2; j++)
            {
                c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
            }
        }
        return c;
    }

    matrix operator *(const matrix &b)const
    {
        matrix c;
        c.clear();
        for(int i=0; i<2; i++)
        {
            for(int j=0; j<2; j++)
            {
                for(int k=0; k<2; k++)
                {
                    c.a[i][j]+=a[i][k]*b.a[k][j];
                }
                c.a[i][j]%=mod;
            }
        }
        return c;
    }

然后如果题目求的是矩阵的前K次方之和，比如说给定一个矩阵A，求A^1+A^2+……A^K

思想：两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k 进行二分。比如，当k=6 时，有：

    A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3)  + A^3*(A + A^2 + A^3)  
应用这个式子后，规模k 减小了一半。我们二分求即可得到原问题的答案。 对于不能整除2的k只需要再乘以一个A^K就可以了

Matrix MatrixSum(int k)//求A+A^2+A^3+....A^k
{
    if(k==1)
        return a;
    Matrix temp,b;
    temp=MatrixSum(k/2);
    if(k&1)
    {
        b=power(k/2+1);
        temp=add(temp,multi(temp,b));
        temp=add(temp,b);
    }
    else
    {
        b=power(k/2);
        temp=add(temp,multi(temp,b));
    }
    return temp;
}

或者这个代码也可以的，比较简洁

matrix pow_sum(matrix a,int k)
{
    matrix ans,temp;
    if(k==1)return a;
    temp=pow_sum(a,k>>1);
    ans=temp+temp*pow(a,k>>1);
    if(k&1)
    {
        ans=ans+pow(a,k);
    }
    return ans;
}

POJ3233

题目连接：http://poj.org/problem?id=3233

#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <cstring>
using namespace std;

#define maxn 101

struct Matrix
{
    int m[maxn][maxn];
} a,per;
int n,mm;

void init()//初始化
{
    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        {
            scanf("%d", &a.m[i][j]);
            a.m[i][j] %=mm;
            per.m[i][j]=(i==j);//单位矩阵
        }
}

Matrix multi(Matrix a, Matrix b)//矩阵乘法
{
    Matrix c;
    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        {
            c.m[i][j]=0;
            for(int k=0; k<n; k++)
            {
                c.m[i][j]+=a.m[i][k]*b.m[k][j];
            }
            c.m[i][j]%=mm;
        }
    return c;
}

Matrix add(Matrix a, Matrix b)//矩阵加法
{
    Matrix c;
    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        {
            c.m[i][j]=(a.m[i][j]+b.m[i][j])%mm;
        }
    return c;
}

Matrix power(int k)//矩阵求幂
{
    Matrix p=a,ans=per;
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p);
            k--;
        }
        else
        {
            k/=2;
            p=multi(p,p);
        }
    }
    return ans;
}

Matrix MatrixSum(int k)//求A+A^2+A^3+....A^k
{
    if(k==1)
        return a;
    Matrix temp,b;
    temp=MatrixSum(k/2);
    if(k&1)
    {
        b=power(k/2+1);
        temp=add(temp,multi(temp,b));
        temp=add(temp,b);
    }
    else
    {
        b=power(k/2);
        temp=add(temp,multi(temp,b));
    }
    return temp;
}

int main()
{
    int k,i,j;
    while(cin>>n>>k>>mm)
    {
        init();
        Matrix ans=MatrixSum(k);
        for(i=0; i<n; i++)
        {
            for(j=0; j<n-1; j++) //输出
                printf("%d ",ans.m[i][j]);
            printf("%d\n",ans.m[i][j]);
        }
    }
    return 520;
}

HDU1575

题目连接：http://acm.hdu.edu.cn/showproblem.php?pid=1575

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
using namespace std;
const int maxn=11;

struct matrax
{
    int m[maxn][maxn];
} a,per;
int n,k;

void init()
{
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
        {
            scanf("%d",&a.m[i][j]);
            per.m[i][j]=(i==j);
        }
    }
}

matrax mul(matrax a,matrax b)
{
    matrax c;
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
        {
            c.m[i][j]=0;
            for(int k=0; k<n; k++)
            {
                c.m[i][j]+=a.m[i][k]*b.m[k][j];
            }
            c.m[i][j]%=9973;
        }
    }
    return c;
}

matrax power(int k)
{
    matrax p=a,ans=per;
    while(k)
    {
        if(k&1)
        {
            ans=mul(ans,p);
            k--;
        }
        else
        {
            k=k/2;
            p=mul(p,p);
        }
    }
    return ans;
}

int main()
{
    //freopen("C:\\Users\\Administrator\\Desktop\\in.txt" , "r" , stdin);
    int t;
    int sum;
    scanf("%d",&t);
    while(t--)
    {
        sum=0;
        scanf("%d%d",&n,&k);
        init();
        matrax ans=power(k);
        /*
        for(int i=0; i<n; i++)
        {
            for(int j=0; j<n-1; j++)
                printf("%d ",ans.m[i][j]);
            printf("%d\n",ans.m[i][n-1]);
        }*/
        for(int i=0; i<n; i++)
        {
            sum+=ans.m[i][i];
            sum%=9973;
        }
        printf("%d\n",sum);
    }
    return 0;
}