POJ 3233 矩阵乘法

题意：求解A+A^2+...+A^k

题解：

1）利用通和公式，原式=(A^k+1 - A)(A - O)^-1 时间复杂度O(n^3lgk)

2）递归求解，A+A^2+...+A^k=(A+A^2+...+A^k/2)+A^k/2(A+A^2+...+A^k/2) 时间复杂度O(n^3lgk^2)

逆矩阵貌似繁琐，直接用第二种方法写的

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <utility>
#include <vector>
#include <queue>
#include <map>
#include <set>
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)>(y)?(y):(x))

using namespace std;

int n,MOD,k;

struct Matrix{
    int n,m;
    vector< vector<int> >a;
    Matrix(){};
    Matrix(const Matrix & T) : n(T.n),m(T.m)
    {
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
                a[i][j]=T.a[i][j];
        }
    }
    Matrix(int N, int M)
    {
        n=N;
        m=M;
        a.resize(N);
        for(int i=0; i<N; i++)
            a[i].resize(M);
    }
    Matrix & operator=(const Matrix &T)
    {
        n=T.n;
        m=T.m;
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
                a[i][j]=T.a[i][j];
        }
        return *this;
    }
    Matrix operator+(const Matrix &T) const
    {
        Matrix tmp(n,m);
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
            tmp.a[i][j]=(a[i][j]+T.a[i][j])%MOD;
        return tmp;
    }
    Matrix operator*(const Matrix &T) const
    {
        Matrix tmp(n,T.m);
        for(int i=0; i<n; i++)
            for(int j=0; j<T.m; j++)
                for(int k=0; k<m; k++)
                tmp.a[i][j]=(tmp.a[i][j]+a[i][k]*T.a[k][j])%MOD;
        return tmp;
    }
    void input(int N, int M)
    {
        n=N;
        m=M;
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
            scanf("%d",&a[i][j]);
        }
    }
    void output()
    {
        for(int i=0; i<n; i++)
        {
            for(int j=0; j<m; j++)
                printf("%d ",a[i][j]);
            printf("
");
        }
    }
    Matrix pow_m(int N)//矩阵满足n=m 矩阵快速幂
    {
        Matrix ret(n,n),tmp(*this);
        for(int i=0; i<n; i++)
            ret.a[i][i]=1;
        while(N)
        {
            if(N&1) ret=ret*tmp;
            tmp=tmp*tmp;
            N>>=1;
        }
        return ret;
    }
}ans,A;

void work(int k)
{
    if(k==1)
    {
        ans=A;
        return;
    }
    if(k==0)
    {
        ans=A.pow_m(0);
        return;
    }
    work(k/2);
    ans=ans*(A.pow_m(0)+A.pow_m(k/2));
    if(k&1) ans=ans+A.pow_m(k);
}

int main()
{
    scanf("%d%d%d",&n,&k,&MOD);
    A.input(n,n);
    work(k);
    ans.output();
    return 0;
}

太刁了

看到这种解法。

|A O|^k+1 =|A^k+1 　　　O|

|E E|　　　　|A^k+...+A^0 E|

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <utility>
#include <vector>
#include <queue>
#include <map>
#include <set>
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)>(y)?(y):(x))

using namespace std;

int n,MOD,k;

struct Matrix{
    int n,m;
    vector< vector<int> >a;
    Matrix(){};
    Matrix(const Matrix & T) : n(T.n),m(T.m)
    {
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
                a[i][j]=T.a[i][j];
        }
    }
    Matrix(int N, int M)
    {
        n=N;
        m=M;
        a.resize(N);
        for(int i=0; i<N; i++)
            a[i].resize(M);
    }
    Matrix & operator=(const Matrix &T)
    {
        n=T.n;
        m=T.m;
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
                a[i][j]=T.a[i][j];
        }
        return *this;
    }
    Matrix operator+(const Matrix &T) const
    {
        Matrix tmp(n,m);
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
            tmp.a[i][j]=(a[i][j]+T.a[i][j])%MOD;
        return tmp;
    }
    Matrix operator*(const Matrix &T) const
    {
        Matrix tmp(n,T.m);
        for(int i=0; i<n; i++)
            for(int j=0; j<T.m; j++)
                for(int k=0; k<m; k++)
                tmp.a[i][j]=(tmp.a[i][j]+a[i][k]*T.a[k][j])%MOD;
        return tmp;
    }
    void input(int N, int M)
    {
        n=N;
        m=M;
        a.resize(n);
        for(int i=0; i<n; i++)
        {
            a[i].resize(m);
            for(int j=0; j<m; j++)
            scanf("%d",&a[i][j]);
        }
    }
    void output()
    {
        for(int i=0; i<n; i++)
            for(int j=0; j<m; j++)
                if(j<m-1) 
                    printf("%d ",a[i][j]);
                else
                    printf("%d
",a[i][j]);
    }
    Matrix pow_m(int N)//矩阵满足n=m 矩阵快速幂
    {
        Matrix ret(n,n),tmp(*this);
        for(int i=0; i<n; i++)
            ret.a[i][i]=1;
        while(N)
        {
            if(N&1) ret=ret*tmp;
            tmp=tmp*tmp;
            N>>=1;
        }
        return ret;
    }
};

int main()
{
    scanf("%d%d%d",&n,&k,&MOD);
    Matrix A(2*n,2*n);
    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        scanf("%d",&A.a[i][j]);

    for(int i=0; i<2*n; i++)
        A.a[n+i%n][i]=1;
    A=A.pow_m(k+1);

    for(int i=0; i<n; i++)
        for(int j=0; j<n; j++)
        {
            if(i==j) A.a[n+i][j]=(A.a[n+i][j]+MOD-1)%MOD;
            if(j<n-1)
                printf("%d ",A.a[n+i][j]);
            else
                printf("%d
",A.a[n+i][j]);
        }
    return 0;
}