POJ3233解题报告

最新推荐文章于 2018-01-31 16:43:00 发布

jhljx

最新推荐文章于 2018-01-31 16:43:00 发布

阅读量503

点赞数

分类专栏：矩阵文章标签： hdu acm algorithm 数学矩阵

本文链接：https://blog.csdn.net/u012294939/article/details/43924809

版权

矩阵专栏收录该内容

21 篇文章 0 订阅

订阅专栏

这篇博客介绍了一道ACM题目POJ3233的解题思路，重点讨论了矩阵快速幂在解决矩阵等比求和问题中的应用。文章强调了在矩阵快速幂运算中取模操作的时间消耗，并提供了两种优化方法：一是将矩阵视为单一元素，利用二分法和矩阵快速幂；二是利用分块矩阵，仅需一次矩阵快速幂运算，避免了二分步骤，提高了效率。

摘要由CSDN通过智能技术生成

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 16466		Accepted: 7012

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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

这道题是一道比较经典的矩阵快速幂运算的题目。在进行矩阵快速幂的过程中我们需要明确程序运行时间主要消耗在什么地方。其实在快速幂运算中最浪费时间的是取模运算。即在进行矩阵乘法时的取模运算很费时间。通过测试取模运算的处理会使得程序运行时间有10倍的差异。这一点是必须得注意的。为了防止超时，所以我们在做矩阵快速幂的题目时一定要注意在取模运算上多一些优化，防止超时。

下面说这道题的解法，这道题是矩阵的等比求和，当然我们可以类比等比数列求和，直接利用公式求解。但是对于矩阵来说，这样会涉及到逆矩阵的计算，还有矩阵是否可逆的判断，比较麻烦，所以我们避开这种做法。比较好的方法有两种，一种是多项式处理的常见方法，即二分法对多项式的和进行处理。构造出来的递推式是n与n/2之间的递推关系式。这也说明了矩阵具有优化递推式的作用。使用矩阵处理递推关系也是很重要的方法。

方法一：

我们可以把一个矩阵看成一个元素，利用二分法+矩阵快速幂求出该矩阵的等比之和。

参考代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<stack>
#include<queue>
#include<ctime>
#include<cstdlib>
#include<iomanip>
#include<utility>
#define pb push_back
#define mp make_pair
#define CLR(x) memset(x,0,sizeof(x))
#define _CLR(x) memset(x,-1,sizeof(x))
#define REP(i,n) for(int i=0;i<n;i++)
#define Debug(x) cout<<#x<<"="<<x<<" "<<endl
#define REP(i,l,r) for(int i=l;i<=r;i++)
#define rep(i,l,r) for(int i=l;i<r;i++)
#define RREP(i,l,r) for(int i=l;i>=r;i--)
#define rrep(i,l,r) for(int i=1;i>r;i--)
#define read(x) scanf("%d",&x)
#define put(x) printf("%d\n",x)
#define ll long long
#define lson l,m,rt<<1
#define rson m+1,r,rt<<11
using namespace std;
int n,k,m;
struct mat
{
    int d[30][30];
} A,E,M;

mat multi(mat &a,mat &b)
{
    mat ans;
    rep(i,0,n)
    {
        rep(j,0,n)
        {
            ans.d[i][j]=0;
            rep(k,0,n)
                ans.d[i][j]+=a.d[i][k]*b.d[k][j];
            ans.d[i][j]%=m;  //在循环外边取模这一步节省了10倍的时间，因为n最多是30,m又小于10^4,所以乘积不会超过10^8,所以乘积之和不会超过3*10^9,由于数据比较弱，所以这一层for循环求出的和可以保证在int范围内
        }
    }
    return ans;
}

mat add(mat &a,mat &b)
{
    mat ans;
    rep(i,0,n)
    {
        rep(j,0,n)
        {
            ans.d[i][j]=a.d[i][j]+b.d[i][j];
            ans.d[i][j]%=m;
        }
    }
    return ans;
}

mat quickmulti(int n)
{
    mat p=E,q=A;
    if(n==0) return p;
    if(n==1) return A;
    while(n!=1)
    {
        if(n&1)
        {
            n--;
            p=multi(p,q);
        }
        else
        {
            n>>=1;
            q=multi(q,q);
        }
    }
    return multi(p,q);
}

mat binarymulti(int k)
{
    if(k==1) return A;
    mat p=binarymulti(k/2),tmp,q;
    if(k&1)          //k为奇数时sum(k)=(1+A^(k/2+1))*sum(k/2)+A^(k/2+1);
    {
        tmp=quickmulti(k/2+1);
        q=multi(tmp,p);
        p=add(p,q);
        p=add(p,tmp);
    }
    else            //k为偶数时sum(k)=(1+A^(k/2))*sum(k/2)
    {
        tmp=quickmulti(k/2);
        q=multi(tmp,p);
        p=add(p,q);
    }
    return p;
}

int main()
{
    scanf("%d%d%d",&n,&k,&m);
    rep(i,0,n)
        rep(j,0,n)
        {
            read(A.d[i][j]);
            A.d[i][j]%=m;   //这是比较关键的一步。题目中m的范围很小，所以提前取模，在矩阵乘法时取模次数就可以减少，防止溢出
            E.d[i][j]=(i==j);
        }
    M=binarymulti(k);
    rep(i,0,n)
    {
        rep(j,0,n)
            printf("%d ",M.d[i][j]);
        printf("\n");
    }
}

方法二：

利用好分块矩阵的性质，构造矩阵。只用一次矩阵快速幂。省掉了方法一中的二分多项式的步骤，使得运算速度有了很大提高。

这个方法中需要用到分块矩阵，这样会有些麻烦。矩阵中套矩阵。不过速度比第一种方法要快一些。

参考代码：

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<stack>
#include<queue>
#include<ctime>
#include<cstdlib>
#include<iomanip>
#include<utility>
#define pb push_back
#define mp make_pair
#define CLR(x) memset(x,0,sizeof(x))
#define _CLR(x) memset(x,-1,sizeof(x))
#define REP(i,n) for(int i=0;i<n;i++)
#define Debug(x) cout<<#x<<"="<<x<<" "<<endl
#define REP(i,l,r) for(int i=l;i<=r;i++)
#define rep(i,l,r) for(int i=l;i<r;i++)
#define RREP(i,l,r) for(int i=l;i>=r;i--)
#define rrep(i,l,r) for(int i=1;i>r;i--)
#define read(x) scanf("%d",&x)
#define put(x) printf("%d\n",x)
#define ll long long
#define lson l,m,rt<<1
#define rson m+1,r,rt<<11
using namespace std;
int n,k,M;
struct mat
{
    int d[30][30];
};
struct Mat
{
    mat m[2][2];
} A,E;

mat add(mat &a,mat &b)
{
    mat ans;
    rep(i,0,n)
    {
        rep(j,0,n)
            ans.d[i][j]=(a.d[i][j]+b.d[i][j])%M;
    }
    return ans;
}

mat multi(mat &a,mat &b)
{
    mat ans;
    rep(i,0,n)
    {
        rep(j,0,n)
        {
            ans.d[i][j]=0;
            rep(k,0,n)
            {
                if(a.d[i][k]&&b.d[k][j])
                    ans.d[i][j]+=a.d[i][k]*b.d[k][j];
            }
            ans.d[i][j]%=M;
        }
    }
    return ans;
}

Mat Multi(Mat &a,Mat &b)
{
    Mat p;
    CLR(p.m);
    rep(i,0,2)
    {
        rep(j,0,2)
        {
            rep(k,0,2)
            {
                mat q=multi(a.m[i][k],b.m[k][j]);
                p.m[i][j]=add(p.m[i][j],q);
            }
        }
    }
    return p;
}

Mat quickmulti(int n)
{
    Mat p=E,q=A;
    if(n==0) return p;
    if(n==1) return A;
    while(n)
    {
        if(n&1)
        {
            n--;
            p=Multi(p,q);
        }
        else
        {
            n>>=1;
            q=Multi(q,q);
        }
    }
    return p;
}

int main()
{
    CLR(E.m);
    rep(i,0,2)
        rep(j,0,30)
    E.m[i][i].d[j][j]=1;
    scanf("%d%d%d",&n,&k,&M);
    rep(i,0,30)
        A.m[0][0].d[i][i]=1;
    rep(i,0,n)
    {
        rep(j,0,n)
        {
            read(A.m[0][1].d[i][j]);
            A.m[0][1].d[i][j]%=M;
            A.m[1][1].d[i][j]=A.m[0][1].d[i][j];
        }
    }
    Mat m1=quickmulti(k);
    rep(i,0,n)
    {
        rep(j,0,n)
            printf("%d ",m1.m[0][1].d[i][j]);
        printf("\n");
    }
}

jhljx

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
POJ3233解题报告

Matrix Power SeriesTime Limit: 3000MS Memory Limit: 131072KTotal Submissions: 16466 Accepted: 7012DescriptionGiven a n × n matrix A and a positive integer k,
复制链接

扫一扫

专栏目录