这题的dp式子应该挺好写的,我是不会告诉你我开始写错了的,设f[i][j]为操作前i次,取到j小球的期望个数(第一维这么大显然不可做),那么
f[i][j]=f[i-1][j](累加)+1*$frac{f[i-1][j-1]}{M}$ - 1* $frac{f[i-1][j]}{M}$(前i-1次拿到的j-1号球转化为j号球)以及(前i-1次拿到的j号球转化为j+1号球)注意1要特殊考虑。移项得
$f[i][j]=(1-1/m)*f[i-1][j]+(1/m)*f[i-1][j-1]$,开始并没有发现他和矩阵快速幂有啥关系,因为矩阵乘的式子是$f[i][j]=∑f[i][k]*f[k][j]$,但是想想矩阵快速幂是怎么优化肥不拉几数列的:
而上面的dp式子化一下和这个很像:第一维发现并没有什么卵用,而且k这么大显然是快速幂的指数啊,所以把它干掉,就成了f[j]=(1-1/m)*f[j]+(1/m)*f[j-1],以n=4为例初始矩阵和转移矩阵分别是这样的:
(为啥没有对齐啊)然后就可以用矩阵快速幂优化递推,复杂度$n^3logk$,然而n是一千啊,时间上不可过,而且即使不考虑时间,1000*1000的数组传参直接RE(亲测),
但是仔细看看转移矩阵有没有什么特点?它是一个循环矩阵!!!所以只需要存下第一行,矩阵乘的时候以$n^2$的复杂度重构两个数组,然后以$n^2$的复杂度计算ans第一行。
#include<iostream> #include<cstring> #include<cstdio> #define LL long long using namespace std; struct jz { double m[2][1010]; }cs; double tmpa[1010][1010],tmpb[1010][1010]; int n,m,k; jz operator * (const jz &a,const jz &b) { jz ans; memset(ans.m,0,sizeof(ans.m)); memset(tmpa,0,sizeof(tmpa)); memset(tmpb,0,sizeof(tmpb)); for(int j=1;j<=n;j++)tmpa[1][j]=a.m[1][j],tmpb[1][j]=b.m[1][j]; for(int i=2;i<=n;i++) for(int j=1;j<=n;j++) if(j==1)tmpa[i][j]=tmpa[i-1][n], tmpb[i][j]=tmpa[i-1][n]; else tmpa[i][j]=tmpa[i-1][j-1],tmpb[i][j]=tmpb[i-1][j-1]; for(int j=1;j<=n;j++) for(int k=1;k<=n;k++) ans.m[1][j]+=tmpa[1][k]*tmpb[k][j]; return ans; } jz operator ^ (jz &a,int &b) { jz ans=a,tem=a;b--; while(b) { if(b&1)ans=ans*tem; tem=tem*tem; b=b>>1; } return ans; } double a[1010],f[1010]; signed main() { // freopen("in.txt","r",stdin); cin>>n>>m>>k; for(int i=1;i<=n;i++) cin>>a[i]; cs.m[1][1]=(double)(1-1.0/m); cs.m[1][2]=(double)1.0/m; cs=cs^k; memset(tmpa,0,sizeof(tmpa)); for(int j=1;j<=n;j++)tmpa[1][j]=cs.m[1][j]; for(int i=2;i<=n;i++) for(int j=1;j<=n;j++) if(j==1)tmpa[i][j]=tmpa[i-1][n]; else tmpa[i][j]=tmpa[i-1][j-1]; for(int i=n;i;i--) { for(int j=1;j<=n;j++) f[i]+=a[j]*tmpa[n-i+1][n-j+1]; } for(int i=1;i<=n;i++) printf("%0.3lf ",f[i]); }