矩阵幂求和

题目大意

给定矩阵 $A$，求 $A+A^2+A^3+...+A^k$ 中每一个数字对 $p$ 取模的结果。

$n \le 30,k \le 10^9,p \le 10^4$。

思路

先把矩阵乘法、加法、快速幂写出来（不做讲解）。

然后写一个函数 $f(x)$，返回值为 $A+A^2+A^3+...+A^x$。

先用 $f(\lfloor x/2 \rfloor)$ 把 $A+A^2+A^3+...+A^{\lfloor x/2 \rfloor}$ 求出来。

假设 $res = f(\lfloor x/2 \rfloor)$，那么

$$\begin{aligned} res \times A^{\lfloor x/2 \rfloor} &= (A+A^2+A^3+...+A^{\lfloor x/2 \rfloor}) \times A^{\lfloor x/2 \rfloor}\\ &= A \times A^{\lfloor x/2 \rfloor}+A^2 \times A^{\lfloor x/2 \rfloor}+A^3 \times A^{\lfloor x/2 \rfloor}+...+A^{\lfloor x/2 \rfloor} \times A^{\lfloor x/2 \rfloor}\\ &= A^{ {\lfloor x/2 \rfloor}+1}+A^{ {\lfloor x/2 \rfloor}+2}+A^{ {\lfloor x/2 \rfloor}+3}+...+A^{ {\lfloor x/2 \rfloor}+{\lfloor x/2 \rfloor} } \end{aligned}$$

这样就可以得出：

$$\begin{aligned} res+res \times A^{\lfloor x/2 \rfloor} &=(A+A^2+A^3+...+A^{\lfloor x/2 \rfloor})+(A^{ {\lfloor x/2 \rfloor}+1}+A^{ {\lfloor x/2 \rfloor}+2}+A^{ {\lfloor x/2 \rfloor}+3}+...+A^{ {\lfloor x/2 \rfloor}+{\lfloor x/2 \rfloor} })\\ &=A+A^2+A^3+...+A^{ {\lfloor x/2 \rfloor\times 2} } \end{aligned}$$

$\, \, \, \, \, \, \, \, res+res \times A^{\lfloor x/2 \rfloor}$

$=(A+A^2+A^3+...+A^{\lfloor x/2 \rfloor})+(A^{ {\lfloor x/2 \rfloor}+1}+A^{ {\lfloor x/2 \rfloor}+2}+A^{ {\lfloor x/2 \rfloor}+3}+...+A^{ {\lfloor x/2 \rfloor}+{\lfloor x/2 \rfloor} })$

$=A+A^2+A^3+...+A^{ {\lfloor x/2 \rfloor\times 2} }$

那么，当 $x$ 为偶数，$\lfloor x/2 \rfloor\times 2=x$，所以 $res+res \times A^{\lfloor x/2 \rfloor}$ 即为答案。

当 $x$ 为奇数，$\lfloor x/2 \rfloor\times 2=x-1$，所以 $res+res \times A^{\lfloor x/2 \rfloor}=A+A^2+A^3+...+A^{x-1}$，所以这个结果加上 $A^x$ 即为答案。

因为每次递归只需要求 $f(\lfloor x/2 \rfloor)$，所以时间复杂度为 $\mathcal{O}(\log k)$。

摆上代码

matrix a;//输入的矩阵

matrix f(int x)
{
    if(x==1) return a;//当x==1，答案就是a.
    matrix res=f(x/2);//求出f(⌊x/2⌋)
    if(x%2==0) return res+(qpow(a,x/2)*res);//当x为偶数，答案为res+res×A^⌊x/2⌋
    else return res+(res*qpow(a,x/2))+qpow(a,x);//否则答案就是 res+res×A^⌊x/2⌋+A^x
}

代码

#include<bits/stdc++.h>
using namespace std;
int n,k,p;
struct matrix
{
    int a[35][35];
    matrix operator *(matrix b)//矩阵乘法
    {
        matrix res;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
            {
                res.a[i][j]=0;
                for(int k=1;k<=n;k++)
                {
                    res.a[i][j]=(res.a[i][j]+a[i][k]*b.a[k][j])%p;
                }
            }
        return res;
    }
    
    matrix operator +(matrix b)//矩阵加法
    {
        matrix res;
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                res.a[i][j]=(a[i][j]+b.a[i][j])%p;
        return res;
    }
};

matrix qpow(matrix a,int k)//快速幂
{
    matrix res;
       for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) res.a[i][j]=0;//清空
    for(int i=1;i<=n;i++) res.a[i][i]=1;//把res变成单位矩阵
    //单位矩阵性质：一个矩阵乘以单位矩阵=这个矩阵。
    while(k)
    {
        if(k&1) res=res*a;
        a=a*a;
        k>>=1;
    }
    return res;
}

matrix a;

matrix f(int x)
{
    if(x==1) return a;//当x==1，答案就是a.
    matrix res=f(x/2);//求出f(⌊x/2⌋)
    if(x%2==0) return res+(qpow(a,x/2)*res);//当x为偶数，答案为res+res×A^⌊x/2⌋
    else return res+(res*qpow(a,x/2))+qpow(a,x);
}

matrix ans;

signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin>>n>>k>>p;
    for(int i=1;i<=n;i++)for(int j=1;j<=n;j++) cin>>a.a[i][j],a.a[i][j]%=p;//读入矩阵
        
    ans=f(k);//求出这个矩阵
    
    for(int i=1;i<=n;i++)//直接输出矩阵就可以了
    {
        for(int j=1;j<=n;j++) cout<<ans.a[i][j]%p<<" ";
        cout<<endl;
    }
    return 0;
}

---