【BZOJ2510】弱题

题目大意

　　有\(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]
*/