这题的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\n",f[i]);
}