POJ3233 Matrix Power Series（快速幂求等比矩阵和）

weixin_30294295

于 2019-02-15 19:09:00 发布

阅读量87

点赞数

原文链接：http://www.cnblogs.com/812-xiao-wen/p/10385473.html

版权

题面

$solution:$

首先，如果题目只要我们求$A^K$ 那这一题我们可以直接模版矩乘快速幂来做，但是它现在让我们求$\sum_{i=1}^{k}{(A^i)} $ 所以我们思考一下这两者是否有什么关系。仔细一想，不难发现几个东西：

一次矩阵乘法复杂度为$O(n^3)$，所以我们不能进行太多次矩阵乘法
快速幂的复杂度为$O(logk)$ 再乘一下矩阵乘法的复杂度，我们现在只能再接受$O(log)$级别的处理了
矩阵乘法满足交换律和结合律！！！！
若我们已经知道了$A^1+A^2+A^3+A^4$ 的值，我们需要$A^5+A^6+A^7+A^8$ 的值，我们可以直接将前者乘上一个$A^4$ 就可以了！

根据以上发现，我们不妨再设一个矩阵B，来帮助我们理解！

我们设$B_x=A^1+A^2+.....+A^x$ 根据上面第四个的原理我们可以得到：

$B_1=A^1$
$B_2=A^1+A^1*A^1=B_1+B_1*A^1$
$B_4=(A^1+A^2)+A^2*(A^1+A^2)=B_2+A^2*B_2$
$B_8=(A^1+A^2+A^3+A^4)+A^4*(A^1+A^2+A^3+A^4)=B_4+A^4*B_4$
$B_{16}=B_8+A^8*B_8$
$B_{2*x}==B_x+A_x*B_x$
在根据上面第四个规律可以得出：$B_{x*y}==B_x+A_x*B_y$

而我们要得到的最终结果就是$B_K$ 嘛。如果上述的B矩阵的底数能相加，那我们不就可以仿照二进制来得到$B_K$ 了吗？就像K等于19，如果我们可以直接$B_{19}=B_{1+2+16}=B_1+B_2+B_{16}$ （假设！）那我们就能log预处理所有$B_{1<<i}$ 然后log求出$B_K$ 了！

我们再思考一下，发现上述$B$ 矩阵的底数是可以相加的，但是不像上面那个式子那么直接相加，我们应该这么求：$B_{x+y}=B_x+A^x*B_y$ 这样底数就可以相加了！！！

然后我们再仿照二进制，就能log求出$B_K$ 了！

$code:$

#include<iostream>
#include<cstdio>
#include<iomanip>
#include<algorithm>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>

#define ll long long
#define db double
#define inf 0x7fffffff
#define rg register int

using namespace std;

int n,m,k;

inline int qr(){ char ch;//快读
    while((ch=getchar())<'0'||ch>'9');
    int res=ch^48;
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+(ch^48);
    return res;
}

struct su{
    int s[30][30];
    inline void read(){//读入一个矩阵
        for(rg i=0;i<n;++i)
            for(rg j=0;j<n;++j)
                s[i][j]=qr();
    }
    inline void write(){//输出一个矩阵
        for(rg i=0;i<n;++i){
            for(rg j=0;j<n;++j)
                printf("%d ",s[i][j]);
            puts("");
        }
    }
    inline su operator *(su x){//矩阵乘法
        su y; memset(y.s,0,sizeof(y.s));
        for(rg i=0;i<n;++i)
            for(rg j=0;j<n;++j)
                for(rg o=0;o<n;++o)
                    y.s[i][j]+=(ll)s[i][o]*x.s[o][j]%m,y.s[i][j]%=m;
        return y;
    }
    inline su operator +(su x){//矩阵加法
        for(rg i=0;i<n;++i)
            for(rg j=0;j<n;++j)
                x.s[i][j]+=s[i][j],x.s[i][j]%=m;
        return x;
    }
}ans,a[33],b[33];

int main(){
    //freopen("t1.in","r",stdin);
    //freopen("t1.out","w",stdout);
    n=qr();k=qr();m=qr();
    a[1].read(); b[1]=a[1];
    for(rg i=1;i<32;++i){//我们求出对应所有的a与b
        a[i+1]=a[i]*a[i]; //a数组表示A[1<<i-1]
        b[i+1]=b[i]*a[i]+b[i]; //b数组表示b[1<<i-1]
    }
    for(rg i=1;i<=32;++i){//根据二进制，不断累加，一直到b[k]
        if(k&1)ans=ans*a[i]+b[i]; //ans相当于b[ans]
        k>>=1; //k&1表示这一位上是一，不懂可以先学下快速幂的原理
    } ans.write();//累加完毕输出
    return 0;
}

我们发现上面main函数中的两个for循环可以换成一个，于是：

int main(){
    n=qr();k=qr();m=qr(); a.read(); b=a;
    while(k){
        if(k&1)ans=ans*a+b;
        b=b*a+b; a=a*a; k>>=1;
    } ans.write(); //不仅代码短，跑的还十分快！！！！！
    return 0;
}

转载于:https://www.cnblogs.com/812-xiao-wen/p/10385473.html