POJ—3233—Matrix power series—【矩阵快速幂】【二分，递归，分治】

Matrix Power Series Time Limit: 3000MS   Memory Limit: 131072K Total Submissions: 15107   Accepted: 6474 Description Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak. Input The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order. Output Output the elements of S modulo m in the same way as A is given. Sample Input 2 2 4 0 1 1 1 Sample Output 1 2 2 3 Source POJ Monthly--2007.06.03, Huang, Jinsong

核心思想参考：

http://www.cnblogs.com/forever4444/archive/2009/05/12/1454736.html

思路:

设矩阵A ，设k=4

A+A^2+A^3+A^4+A^4

=(A+A^2)(1+A^2)

=A(1+A)(1+A^2)

意思就是A+……+A^k=(A+……A^k/2)(1+A^k/2)

如果k是奇数，先按上面的方法算A^k-1，最后加上A^k

【代码实现中有些不同，分治的，但是我不晓得怎么说明。。。o(╯□╰)o】

所以当k=4，要先求出(A+A^2)，

要求(A+A^2)，又要求出A，

这是一个递归的过程。。。

【好吧，承认自己说不清楚。。。对于递归也比较拎不清。。。】

直接上代码，模块，内容跟连接里的非常非常相似，但是我的跑的很慢，600+MS，链接里的100+MS。。。POJ上有大神0MS。。。膜拜啊。。。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <algorithm>
using namespace std;
#define LL __int64

typedef struct MATRIX
{
    int mat[35][35];
}MATRIX;
MATRIX A;
int n,m;

MATRIX mat_multiply (MATRIX a,MATRIX b)     //矩阵乘法
{
    MATRIX c;
    memset(c.mat,0,sizeof(c.mat));
    for(int k=0;k<n;k++)
        for(int i=0;i<n;i++)
            if(a.mat[i][k])
                for(int j=0;j<n;j++)
                    if(b.mat[k][j])
                    {
                        c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                        if(c.mat[i][j]>=m)
                            c.mat[i][j]%=m;
                    }
    return c;
}

MATRIX mat_add (MATRIX a,MATRIX b)          //矩阵加法
{
    MATRIX c;
    memset(c.mat,0,sizeof(c.mat));
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
        {
            if(a.mat[i][j] || b.mat[i][j])
            {
                c.mat[i][j]=a.mat[i][j] + b.mat[i][j];
                if(c.mat[i][j]>=m)
                    c.mat[i][j]%=m;
            }
        }
    return c;
}

MATRIX mat_quick_index (int N)           //矩阵求幂
{
    MATRIX E,a;
    a=A;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            E.mat[i][j]=(i==j);
    while(N)
    {
        if(N & 1)
            E=mat_multiply(E,a);
        N>>=1;
        a=mat_multiply(a,a);
    }
    return E;
}

MATRIX mat_sum (int k)
{
    if(k==1)
        return A;
    MATRIX temp,tnow;
    temp=mat_sum(k/2);
    if(k & 1)
    {
        tnow=mat_quick_index(k/2+1);
        temp=mat_add(temp,mat_multiply(temp,tnow));
        temp=mat_add(temp,tnow);
    }
    else
    {
        tnow=mat_quick_index(k/2);
        temp=mat_add(temp,mat_multiply(temp,tnow));
    }
    return temp;
}

int main()
{
    MATRIX B;
    int k;
    memset(B.mat,0,sizeof(B.mat));
    scanf("%d%d%d",&n,&k,&m);
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
        {
            scanf("%d",&A.mat[i][j]);
            A.mat[i][j]%=m;     //初始化的时候就先取余。。。
        }
    B=mat_sum(k);
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        {
            printf("%d",B.mat[i][j]);
            if(j!=n-1)
                printf(" ");
        }
        printf("\n");
    }
    return 0;
}

另外再贴一个方法吧。。。还没仔细看，但是感觉不错的。。。500+MS，但是代码量少了不少！！！

http://www.xuebuyuan.com/574359.html

这道题目是第三种矩阵乘法的应用。

然我们求 S =  A +  A ² + A ³ + … +
A^k.的结果矩阵

M67大牛的博客上面讲的是一种分治的方法。

设f[n]=A^1+A^2+....A^n;
当n是偶数,f[n]=f[n/2]+f[n/2]*A^(n/2);
但n是奇数,f[n]=f[n-1]+A^(n);

而看了DISCUSS之后发现一种编程复杂度更加小的方法。再次涨了姿势。。

Let B=[ [A I];
             [0 I] ]

B^(k+1) = [ [A^k   I+A+...+A^k];
                   [0                I        ] ]

只能说矩阵真是太神奇了。。。。。。。

我比较懒。所以。。。实现了第二种。。。。

#include<iostream>
#include<cstdio>

#define MAXN 100

using namespace std;

struct Matrix
{
    int size;
    long modulo;
    long element[MAXN][MAXN];
    void setSize(int);
    void setModulo(int);
    Matrix operator* (Matrix);
    Matrix power(int);

};

void Matrix::setSize(int a)
{
    for (int i=0; i<a; i++)
        for (int j=0; j<a; j++)
            element[i][j]=0;
    size = a;
}

void Matrix::setModulo(int a)
{
    modulo = a;
}

Matrix Matrix::operator* (Matrix param)
{
    Matrix product;
    product.setSize(size);
    product.setModulo(modulo);
    for (int i=0; i<size; i++)
        for (int j=0; j<size; j++)
            for (int k=0; k<size; k++)
            {
                product.element[i][j]+=element[i][k]*param.element[k][j];
                product.element[i][j]%=modulo;
            }
    return product;
}

Matrix Matrix::power(int exp)
{
    Matrix tmp = (*this) * (*this);
    if (exp==1)
        return *this;
    else if (exp & 1)
        return tmp.power(exp/2) * (*this);
    else
        return tmp.power(exp/2);
}

int n,m,k;
Matrix m1,ans;
int main()
{
    while(~scanf("%d%d%d",&n,&k,&m))
    {
        m1.setModulo(m);
        m1.setSize(n*2);
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                scanf("%d",&m1.element[i][j]);
                m1.element[i][j]%=m;
            }
            m1.element[i][i+n]=1;
            m1.element[n+i][n+i]=1;
        }
        m1.element[n][n]=1;
        ans=m1.power(k+1);
        for(int i=0;i<n;i++)
            for(int j=n;j<n*2;j++)
            {
                int out=ans.element[i][j];
                if(j==i+n)
                    out=out==0?m-1:out-1;
                printf("%d%s",out,j!=n*2-1?" ":"\n");
            }
        return 0;

    }
    return 0;
}