UVA11149 Power of Matrix（快速幂求等比矩阵和）

Posted 2021-02-09 812-xiao-wen

题面

---

---

(solution:)

首先这一题是(UVA11149)的题目，建议到(vjudge)上去测，没办法(UVA)太难注册了。然后其原题与本题不是完全一样的，毒瘤卡输出！但思想一模一样！

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

1. 一次矩阵乘法复杂度为(O(n^3))，所以我们不能进行太多次矩阵乘法
2. 快速幂的复杂度为(O(logk)) 再乘一下矩阵乘法的复杂度，我们现在只能再接受(O(log))级别的处理了
3. 矩阵乘法满足交换律和结合律！！！！
4. 若我们已经知道了(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) 根据上面第四个的原理我们可以得到：

1. (B^1=A^1)
2. (B^2=A^1+A^1*A^1=B^1+B^1*A^1)
3. (B^4=(A^1+A^2)+A^2*(A^1+A^2)=B^2+A^2*B^2)
4. (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)
5. (B^{16}==B^8+A^8*B^8)

而我们要得到的最终结果就是(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;
}