【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)$。

然后就做完了。

完结撒花

 1 #include
 2 #define M 55
 3 using namespace std;
 4 int MOD;
 5 struct mat{
 6     int a[M][M],n,m; mat(){memset(a,0,sizeof(a));}
 7     mat(int nn,int mm){n=nn; m=mm; memset(a,0,sizeof(a));}
 8     int* operator [](int x) {return a[x];}
 9     friend mat operator *(mat a,mat b){
10         mat c=mat(a.n,b.m);
11         for(int i=1;i<=c.n;i++)
12         for(int j=1;j<=c.m;j++)
13         for(int k=1;k<=b.n;k++)
14         c[i][j]=(c[i][j]+a[i][k]*b[k][j])%MOD;
15         return c;
16     }
17     void dw(){memset(a,0,sizeof(a)); for(int i=1;i<=n;i++) a[i][i]=1;}
18     friend mat operator ^(mat a,int b){
19         mat ans=mat(a.n,a.m); ans.dw();
20         while(b){
21             if(b&1) ans=ans*a;
22             a=a*a; b>>=1;
23         }
24         return ans;
25     }
26 };
27 
28 int n,m;
29 int a[M]={0},A1,A2,B1,B2;
30 int main(){
31     cin>>n>>m>>MOD;
32     cin>>a[1]>>B1>>B2;cin>>A1>>A2;
33     for(int i=2;i<=n;i++) a[i]=(a[i-1]*A1+A2)%MOD;
34     
35     mat f=mat(n+2,n+2),g=mat(1,n+2);
36     for(int i=1;i<=n;i++) g[1][i]=a[i]; 
37     g[1][n+1]=a[1]; g[1][n+2]=1;
38     
39     for(int i=1;i<=n;i++) for(int j=1;j<=i;j++) f[j][i]=1;
40     f[n+1][1]=f[n+1][n+1]=B1;
41     f[n+2][1]=f[n+2][n+1]=B2;
42     f[n+2][n+2]=1;
43     
44     int q; scanf("%d",&q);
45     while(q--){
46         int x,y; scanf("%d%d",&x,&y);
47         if(y==1) {printf("%d\n",a[x]); continue;}
48         mat F=f^(y-2);
49         mat ans=g*F;
50         if(x==1) {
51             ans=ans*f;
52             printf("%d\n",ans[1][n+1]); 
53             continue;
54         }
55         int sum=0;
56         for(int i=1;i1][i])%MOD;
57         printf("%d\n",sum);
58     }
59 }