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

题面

\(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;
}

原文地址：https://www.cnblogs.com/812-xiao-wen/p/10385473.html