HDU 4965 Fast Matrix Calculation 矩阵快速幂

链接：http://acm.hdu.edu.cn/showproblem.php?pid=4965

题意：一个矩阵N*K的矩阵A，一个K*N的矩阵B，(4 <= N <= 1000,2 <=K <= 6)。先进行矩阵乘法C=A*B，然后进行D=C^（N*N）的幂运算，然后对D的每个数对6取模，最后求出D的所有位置数字之和。

思路：像之前那道矩阵乘法一样，特别大的矩阵直接进行乘法在没有小规律的帮助时是不可能直接过的（目前看即使是Strassen矩阵算法也不会加速到要求以内）题目中给的C矩阵是1000*1000的矩阵进行快速幂是一定超时的，所以我注意到了A矩阵的列数和B矩阵的行数K是不超过六的，而6*6的矩阵进行快速幂是没问题的。D=A*B*A*B*...*A*B*A*B，可以看成D=A*(B*A*...*B*A)*B。这样既不影响矩阵左乘和右乘的关系，又能防止行列数太大爆掉。（题里的小细节是挺厉害的，首先爆栈的问题需要解决，然后还有矩阵的赋值运算也是O(N)的复杂度）

代码：

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <ctype.h>
#include <iostream>
#include <map>
#include <queue>
#include <set>
#include <stack>
#include <string>
#include <vector>
#define eps 1e-8
#define INF 0x7fffffff
#define PI acos(-1.0)
#define seed 31//131,1313
typedef long long LL;
typedef unsigned long long ULL;
#define MOD 6
#define maxn 1005
#define maxm 1005
using namespace std;
struct Matrix
{
    int n,m;
    int a[maxn][maxm];
    void init()
    {
        n=m=0;
        memset(a,0,sizeof(a));
    }
    Matrix operator +(const Matrix &b) const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
                tmp.a[i][j]=a[i][j]+b.a[i][j];
        return tmp;
    }
    Matrix operator -(const Matrix &b) const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=m;
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
                tmp.a[i][j]=a[i][j]-b.a[i][j];
        return tmp;
    }
    Matrix operator *(const Matrix &b) const
    {
        Matrix tmp;
        tmp.n=n;
        tmp.m=b.m;
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<b.m;j++)
                tmp.a[i][j]=0;
        }
        for(int i=0; i<n; i++)
            for(int j=0; j<b.m; j++)
                for(int k=0; k<m; k++)
                {
                    tmp.a[i][j]+=a[i][k]*b.a[k][j];
                }
        for(int i=0; i<n; i++)
            for(int j=0; j<b.m; j++)
                tmp.a[i][j]%=MOD;
        return tmp;
    }
    void Copy(const Matrix &x)
    {
        n=x.n;
        m=x.m;
        for(int i=0;i<n;i++)
            for(int j=0;j<m;j++)
                a[i][j]=x.a[i][j];
    }
};//只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义
Matrix M_quick_pow(Matrix m,int k)
{
    Matrix tmp;
    tmp.n=m.n;
    tmp.m=m.m;//m=n才能做快速幂
    for(int i=0; i<tmp.n; i++)
    {
        for(int j=0; j<tmp.n; j++)
        {
            if(i==j)
                tmp.a[i][j]=1;
            else tmp.a[i][j]=0;
        }
    }
    while(k)
    {
        if(k&1)
            tmp.Copy(tmp*m);
        k>>=1;
        m.Copy(m*m);
    }
    return tmp;
}
Matrix a,b;
int main()
{
    int n,k;
    while(scanf("%d%d",&n,&k))
    {
        if(n==0&&k==0)
            break;
        a.n=n;
        a.m=k;
        for(int i=0; i<a.n; i++)
            for(int j=0; j<a.m; j++)
                scanf("%d",&a.a[i][j]);
        b.n=k;
        b.m=n;
        for(int i=0; i<b.n; i++)
            for(int j=0; j<b.m; j++)
                scanf("%d",&b.a[i][j]);
        Matrix e,c;
        c.Copy(b*a);
        c.Copy(M_quick_pow(c,n*n-1));
        e.Copy(a*c);
        e.Copy(e*b);
        LL ans=0;
        for(int i=0;i<e.n;i++)
        {
            for(int j=0;j<e.m;j++)
            {
                ans+=(LL)(e.a[i][j]);
            }
        }
        printf("%I64d\n",ans);
    }
    return 0;
}