ココロコネクト・ヒトランダムの問題
こんな問題
とりあえず思いつく解法としては、それぞれの状態になる確率を記憶しておいての動的計画法。
// 順列を一意な整数に変換 int perm2id(vector<int>v) { int N=(int)v.size(); int r=0; for (int i=0;i<N;i++){ int t=v[i]; for(int j=0;j<i;j++)if(v[j]<v[i])t--; r=r*(N-i)+t;} return r; } // 一意な整数から順列に変換 vector<int>id2perm(int N,int r) { vector<int>v(N); for(int i=N-1;i>=0;i--)v[i]=r%(N-i),r/=N-i; vector<bool>f(N); for(int i=0;i<N;i++){ int t=0; while(f[t])t++; for(int j=0;j<v[i];j++){t++;while(f[t])t++;} v[i]=t;f[t]=true;} return v; } double solve( int N, int k ) { int L = 1; for ( int i=1; i<=N; i++ ) L *= i; vector<double> T(L); T[0] = 1; for ( int i=0; i<k; i++ ) { vector<double> P = T; T = vector<double>(L); for ( int j=0; j<L; j++ ) { vector<int> v = id2perm(N,j); for ( int a=0; a<N; a++ ) for ( int b=0; b<a; b++ ) { swap(v[a],v[b]); T[perm2id(v)] += P[j]/(N*(N-1)/2); swap(v[a],v[b]); } } } return T[0]; }
計算量(kN3N!)で、N=10くらいが限界。
N k prob time[s] --- --- -------------- ------ 2 10 1.000000000000 0.00 3 10 0.333333333333 0.00 4 10 0.083346034649 0.00 5 10 0.016927126000 0.00 6 10 0.003201527728 0.02 7 10 0.000662616042 0.27 8 10 0.000159123443 3.28 9 10 0.000044105344 44.39 10 10 0.000013859489 640.31
N=10で640秒だから、N=11なんて計算すると、おねえさん死んじゃう！やめて！
指数的であることに変わりはないけど、もう少し速い方法。部員の置換を相異なる巡回置換に分解したときの、巡回置換の大きさがそれぞれ等しければ、その状態になる確率は等しい。例えば、
1→2→3→1 4→4 5→6→5
という入れ替わりと、
1→1 2→6→3→2 4→6→4
という入れ替わりが発生する確率は等しい。整数Nを順番の違いを除いて和で表す方法の総数（分割数）は階乗に比べると小さく、N=42でも37338にしかならない。これを使うとN=40くらいまで計算できる。
// Nの分割を一意な整数に変換・逆変換 class Part { vector<vector<int> >P; void f(int N,int k,vector<int>&v){ if(N==0)P.push_back(v); for(int i=k;i<=N;i++){ v.push_back(i); f(N-i,i,v); v.pop_back();}} public: Part(int N){ vector<int>v; f(N,1,v);} int size()const{ return(int)P.size();} int part2id(vector<int>v)const{ int r=0,b=1; while(b<(int)P.size())b*=2; for(;b>0;b/=2) if(r+b<(int)P.size()&&P[r+b]<=v) r+=b; return r;} vector<int>id2part(int r)const{ return P[r];} }; double solve( int N, int k ) { Part Pt(N); vector<double> T(Pt.size()); T[0] = 1; for ( int i=0; i<k; i++ ) { vector<double> P = T; T = vector<double>(Pt.size()); for ( int j=0; j<Pt.size(); j++ ) { vector<int> p = Pt.id2part(j); // 異なる巡回置換に属する要素を交換する場合 for ( int a=0; a<(int)p.size(); a++ ) for ( int b=0; b<a; b++ ) { vector<int> q = p; q.erase(q.begin()+a); q.erase(q.begin()+b); q.push_back(p[a]+p[b]); sort(q.begin(),q.end()); T[Pt.part2id(q)] += P[j]*p[a]*p[b]/(N*(N-1)/2); } // 同じ巡回置換に属する要素を入れ替える場合 for ( int a=0; a<(int)p.size(); a++ ) for ( int k=1; k<p[a]; k++ ) { vector<int> q = p; q.erase(q.begin()+a); q.push_back(k); q.push_back(p[a]-k); sort(q.begin(),q.end()); T[Pt.part2id(q)] += P[j]*p[a]/2/(N*(N-1)/2); } } } return T[0]; }
実行結果。
N k prob time[s] --- --- -------------- ------ 2 10 1.000000000000 0.00 3 10 0.333333333333 0.00 4 10 0.083346034649 0.00 5 10 0.016927126000 0.00 6 10 0.003201527728 0.00 7 10 0.000662616042 0.00 8 10 0.000159123443 0.00 9 10 0.000044105344 0.00 10 10 0.000013859489 0.00 11 10 0.000004847868 0.00 12 10 0.000001857739 0.01 13 10 0.000000769656 0.01 14 10 0.000000341006 0.02 15 10 0.000000160139 0.02 16 10 0.000000079123 0.04 17 10 0.000000040881 0.05 18 10 0.000000021974 0.07 19 10 0.000000012235 0.10 20 10 0.000000007031 0.14 21 10 0.000000004157 0.21 22 10 0.000000002522 0.27 23 10 0.000000001566 0.37 24 10 0.000000000993 0.53 25 10 0.000000000642 0.71 26 10 0.000000000423 0.91 27 10 0.000000000283 1.23 28 10 0.000000000193 1.83 29 10 0.000000000133 2.13 30 10 0.000000000093 2.79 31 10 0.000000000066 4.81 32 10 0.000000000047 5.67 33 10 0.000000000034 6.18 34 10 0.000000000025 7.87 35 10 0.000000000018 10.29 36 10 0.000000000014 12.98 37 10 0.000000000010 16.43 38 10 0.000000000008 30.46 39 10 0.000000000006 26.77 40 10 0.000000000005 34.07
ちなみに、kをもう少し大きくするとこんな感じで、2/N!になる。
N k prob time[s] --- --- -------------- ------ 2 100 1.000000000000 0.00 3 100 0.333333333333 0.00 4 100 0.083333333333 0.00 5 100 0.016666666667 0.00 6 100 0.002777777778 0.00 7 100 0.000396825397 0.01 8 100 0.000049603175 0.01 9 100 0.000005511464 0.02 10 100 0.000000551146 0.03
kが奇数なら、0。偶置換の個数がN!/2で、それぞれが等確率で発生するということなのだと思う。
