【xsy1103】随机数表（RanMat）矩阵快速幂

题目大意：你生成了一个随机数表，生成机制是这样子的：

$a[i]=A1a[i-1]+A2(2≤i≤m)$

$b[i]=B1b[i-1]+B2(2≤i≤m)$

$M[1][y]=a[y]%P，(1≤y≤m)$

$M[x][1]=b[x]%P，(2≤x≤n)$

$M[x][y]=(\sum\limits_{i=1}^{x-1}\sum\limits_{j=1}^{y-1} M[i][j])%P，(2≤x≤n，2≤y≤m)$

有$k$组询问，每次问你$M[x][y]$的值。

数据范围：$n≤50$，$m≤10^9$，$k≤10$，$P≤32768$

 

我们考虑$y=1$和$x=1$的情况，这两种情况直接等于$a$或$b$，直接矩阵快速幂就可以了。

对于非这两种的情况，我们考虑一个$1\times n+2$的矩阵

我们用这个矩阵的前$n$个数表示第i行的前缀和，第$n+1$个数为$M[1][j]$的值，第$n+2$个数恒为$1$，大概长这样：

$\begin{bmatrix}
s(1,i),s(2,i)\cdots s(n-1,i),s(n,i),a[i],1
\end{bmatrix}$

其中$s(x,y)=\sum\limits_{i=1}^{y} M[x][i]$

然后，我们考虑构造一个矩阵，使得上面这个矩阵乘上它后，可以变成

$\begin{bmatrix}
s(1,i+1),s(2,i+1)\cdots s(n-1,i+1),s(n,i+1),a[i+1],1
\end{bmatrix}$

不难推出这个矩阵是长这样的：

$\begin{bmatrix}1 , 1 , 1 ,\cdots 1 , 0 , 0\\0,1,1,\cdots 1,0,0\\0,0,1,\cdots 1,0,0\\ \vdots \ \ \ \ \ddots \ \ \ \ \ \ \vdots
\\ 0,0,0,\cdots 1,0,0\\
A1,0,0,\cdots A1,0\\
A2,0,0,\cdots A2,0\\ \end{bmatrix} $

假设我们需要求$M[x][y]$，我们可以通过矩阵快速幂，先求出

$\begin{bmatrix}
s(1,y-1),s(2,y-1)\cdots s(n-1,y-1),s(n,y-1),a[y-1],1
\end{bmatrix}$

然后$M[x][y]$显然等于$\sum\limits_{i=1}^{x-1} s(i,y-1)$。

然后就做完了。

完结撒花

#include<bits/stdc++.h>
#define M 55
using namespace std;
int MOD;
struct mat{
    int a[M][M],n,m; mat(){memset(a,0,sizeof(a));}
    mat(int nn,int mm){n=nn; m=mm; memset(a,0,sizeof(a));}
    int* operator [](int x) {return a[x];}
    friend mat operator *(mat a,mat b){
        mat c=mat(a.n,b.m);
        for(int i=1;i<=c.n;i++)
        for(int j=1;j<=c.m;j++)
        for(int k=1;k<=b.n;k++)
        c[i][j]=(c[i][j]+a[i][k]*b[k][j])%MOD;
        return c;
    }
    void dw(){memset(a,0,sizeof(a)); for(int i=1;i<=n;i++) a[i][i]=1;}
    friend mat operator ^(mat a,int b){
        mat ans=mat(a.n,a.m); ans.dw();
        while(b){
            if(b&1) ans=ans*a;
            a=a*a; b>>=1;
        }
        return ans;
    }
};

int n,m;
int a[M]={0},A1,A2,B1,B2;
int main(){
    cin>>n>>m>>MOD;
    cin>>a[1]>>B1>>B2;cin>>A1>>A2;
    for(int i=2;i<=n;i++) a[i]=(a[i-1]*A1+A2)%MOD;
    
    mat f=mat(n+2,n+2),g=mat(1,n+2);
    for(int i=1;i<=n;i++) g[1][i]=a[i]; 
    g[1][n+1]=a[1]; g[1][n+2]=1;
    
    for(int i=1;i<=n;i++) for(int j=1;j<=i;j++) f[j][i]=1;
    f[n+1][1]=f[n+1][n+1]=B1;
    f[n+2][1]=f[n+2][n+1]=B2;
    f[n+2][n+2]=1;
    
    int q; scanf("%d",&q);
    while(q--){
        int x,y; scanf("%d%d",&x,&y);
        if(y==1) {printf("%d\n",a[x]); continue;}
        mat F=f^(y-2);
        mat ans=g*F;
        if(x==1) {
            ans=ans*f;
            printf("%d\n",ans[1][n+1]); 
            continue;
        }
        int sum=0;
        for(int i=1;i<x;i++) sum=(sum+ans[1][i])%MOD;
        printf("%d\n",sum);
    }
}