zoukankan      html  css  js  c++  java
  • POJ3233:Matrix Power Series(矩阵快速幂+递推式)

    传送门

    题意

    给出n,m,k,求

    [sum_{i=1}^kA^i ]

    A是矩阵

    分析

    我们首先会想到等比公式,然后得到这样一个式子:

    [frac{A^{k+1}-E}{A-E} ]

    发现要用矩阵除法,可以用求矩阵逆来做,现在我们换一种做法,我们发现有这样一个性质:

    [left[ egin{matrix} A & E \ 0 & E \ end{matrix} ight]^n= left[ egin{matrix} A^{n} & sum_{i=0}^{n-1}A^i \ 0 & E \ end{matrix} ight] ]

    那么我们将原先矩阵扩大成四倍,对矩阵求k+1次幂,然后取右上角减去一个单位阵即可

    trick

    代码

    #include<cstdio>
    #include<cstring>
    #include<iostream>
    using namespace std;
    struct Matrix
    {
        int matrix[80][80];
    }ans,ret;
    int t,n,k,mod;
    Matrix multi(Matrix a,Matrix b)
    {
        Matrix tmp;
        for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
        {
            tmp.matrix[i][j]=0;
            for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
        {
            for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
         for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
        {
            tmp.matrix[i][j]=0;
            for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
        {
            for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
        {
            tmp.matrix[i][j]=0;
            for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
        {
            for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=n+1;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
        {
            tmp.matrix[i][j]=0;
            for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        for(int i=1+n;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
        {
            for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
        }
        return tmp;
    }
    void fast_mod(int p)
    {
        memset(ans.matrix,0,sizeof(ans.matrix));
        for(int i=1;i<=2*n;++i) ans.matrix[i][i]=1;
        for(;p;p>>=1,ret=multi(ret,ret)) if(p&1) ans=multi(ans,ret);
    }
    
    
    int main()
    {
        //freopen("data.in","w",stdout);
        while(scanf("%d %d %d",&n,&k,&mod)!=EOF)
        {
            for(int i=1;i<=n;++i)for(int j=1;j<=n;++j) scanf("%d",&ret.matrix[i][j]);
            for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
            for(int i=n+1;i<=2*n;++i)for(int j=1+n;j<=2*n;++j) if(i==j) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
            for(int i=n+1;i<=2*n;++i)for(int j=1;j<=n;++j)   ret.matrix[i][j]=0;
            //for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ret.matrix[i][j],j==2*n?'
    ':' ');
            fast_mod(k+1);
            //for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'
    ':' ');
            for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) (ans.matrix[i][j]+=(mod-1))%=mod;
            for(int i=1;i<=n;++i)for(int j=1+n;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'
    ':' ');
        }
        return 0;
    }
    
  • 相关阅读:
    VBS处理AD帐号密码到期提醒的脚本[zt]
    简单几步手工扩容LVM、缩小LVM及移除磁盘(笔记)
    python加入进度条:tqdm 和 progressbar
    python的map和reduce函数
    python的lambda表达式
    python的推导式 —— 列表推导式、集合和字典推导式
    Pyspark中遇到的 java.io.IOException: Not a file 和 pyspark.sql.utils.AnalysisException: 'Table or view not found
    pyecharts绘制map地图
    pyecharts绘制geo地图
    sklearn.feature_extraction.text 的TfidfVectorizer函数
  • 原文地址:https://www.cnblogs.com/chendl111/p/6690890.html
Copyright © 2011-2022 走看看