题目大意
有\(M\)个球，一开始每个球均有一个初始标号，标号范围为\(1-N\)且为整数，标号为i的球有\(a_i\)个，并保证\(\sum a_i=M\)。
每次操作等概率取出一个球（即取出每个球的概率均为\(1/M\)），若这个球标号为\(k(k < N)\)，则将它重新标号为\(k + 1\)；若这个球标号为\(N\)，则将其重标号为\(1\)。（取出球后并不将其丢弃）
现在你需要求出，经过K次这样的操作后，每个标号的球的期望个数。
\(N ≤ 1000, M ≤ 100,000,000, K ≤ 2,147,483,647\)
题目分析
递推方程很好想，但是需要优化。
由于\(K\)很大，考虑使用矩阵乘法优化\(dp\)。但是\(n\)的范围太大，不能直接化出矩阵来相乘。通过观察矩阵，我们很容易发现矩阵的每一行之间是循环的，因此，我们可以只算一行即可，矩阵相乘的时间负责度瞬间降低为\(O(n^2)\)，矩阵快速幂优化后总时间复杂为\(O(n^2log_2 n)\)
#include <bits/stdc++.h>
using namespace std;
int n,m,k;
namespace Task1{
double p[2][1005];
void solve(){
for(int i=1;i<=n;i++)scanf("%lf",&p[0][i]);
int cur=0,pre;
for(int i=1;i<=k;i++){
pre=cur;cur^=1;
for(int j=1;j<=n;j++){
if(j==1)p[cur][j]=p[pre][j]+p[pre][n]/m-p[pre][j]/m;
else if(j==n)p[cur][j]=p[pre][j]+p[pre][j-1]/m-p[pre][j]/m;
else p[cur][j]=p[pre][j]+p[pre][j-1]/m-p[pre][j]/m;
}
}
for(int i=1;i<=n;i++)cout<<fixed<<setprecision(3)<<p[cur][i]<<"\n";
}
}
namespace Task2{
double tmp[1005][1005];
struct node{
double a[1005];
node(){memset(a,0,sizeof(a));}
double&operator[](int x){return a[x];}
node operator*(node &b){
node c;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++){
int pos=i+j-1;
if(pos>n)pos-=n;
tmp[j][pos]=b[i];
}
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
c[i]+=a[j]*tmp[j][i];
return c;
}
node operator^(int cnt){
node ret,mul=*this;
bool flag=0;
for(;cnt;cnt>>=1,mul=mul*mul)if(cnt&1){if(!flag)ret=mul,flag=1;else ret=ret*mul;}
return ret;
}
}p;
void solve(){
for(int i=1;i<=n;i++)scanf("%lf",&p[i]);
node mul;
mul[1]=(m-1)/double(m);
mul[2]=1.0/m;
mul=mul^k;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++){
int pos=i+j-1;
if(pos>n)pos-=n;
tmp[j][pos]=mul[i];
}
node ans;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
ans[i]+=p[j]*tmp[j][i];
for(int i=1;i<=n;i++)printf("%.3lf\n",ans[i]);
}
}
int main(){
scanf("%d%d%d",&n,&m,&k);
if(1ll*n*k<=5e7)Task1::solve();
else Task2::solve();
}
/*
[(m-1)/m,1/m,0,0,0]
[0,(m-1)/m,1/m,0,0]
[0,0,(m-1)/m,1/m,0]
[0,0,0,(m-1)/m,1/m]
[1/m,0,0,0,(m-1)/m]
[a1,a2,a3] [a1,a2,a3] [a1*a1+a2*a3+a3*a2,a1*a2+a2*a1+a3*a3,a1*a3+a2*a2+a3*a1]
[a3,a1,a2]*[a3,a1,a2]=[a3*a1+a1*a3+a2*a2,......
[a2,a3,a1] [a2,a3,a1] [......
[a1,a2,a3] [b1,b2,b3] [a1*b1+a2*b3+a3*b2,a1*b2+a2*b1+a3*b3,a1*b3+a2*b2+a3*b1]
[a3,a1,a2]*[b3,b1,b2]=[a3*b1+a1*b3+a2*b2,......
[a2,a3,a1] [b2,b3,b1] [......
[a1,a2,a3]*[b1,b2,b3]=>[a1*b1+a2*b3+a3*b2,a1*b2+a2*b1+a3*b3,a1*b3+a2*b2+a3*b1]
*/