"蔚来杯"2022牛客暑期多校训练营3_ACM/NOIP/NOI/CCPC/ICPC算法编程基础集训_牛客竞赛OJ
注(検索用):Multi-Universities Campとか、TwitterではNowCoderとか呼ばれている中国の合宿の問題です。
問題概要
長さの順列に対し、をなる()の個数として定めます。をランダムに決めたときのの期待値をで求めてください。
解法
準備
まずの場合は自明なので除きます。また素因数分解するとなので、素数についてで求めることができれば、中国剰余定理より答えが求まります。以降はこれらのことを前提にして解きます。
を求める
解説スライドによれば、を固定して完全順列の式を使い式変形を頑張る方法もあるようですが、ここでは組み合わせ的解釈から導くことにします。
をだと思うことにします。これはつまり、であるようなを一つ選ぶ場合の数です。今を求めているのですから、一つ取り出す操作を回繰り返す方法を数え上げればよいです。このとき取り出す順番を考慮していることに注意します。
さて、個取り出したたちは当然大量に重複しますが、この重複を除いたものが種類(ただし)だったとします。このを固定してみます。
まず、種類がどのように分布するかを考えます。回の操作をそれぞれ区別される玉、種類のを区別しない箱だと思えば、求める分布の数はそれぞれの箱に1個以上の玉を入れる場合の数ですから、写像12相によればとなります。ここでは第2種スターリング数です。
次に、それを達成できる順列がどれくらいあるか考えます。からのうち種類の箱として使うものを決めるのが、どの箱にどの数字を対応させるかまで含めて通りになります。先ほど箱を区別しなかった分、ここで区別する必要がありました。選んだ個の値は位置まで決定しており、残りは自由ですから通りで、掛け合わせるとになります。
順列としての重複を考える必要はありません。なぜなら今の考察は、最初「すべての順列について」「種類の取り方すべてについて」の順で考えていたものを、逆順で表現しなおしただけだからです。
よってであることがわかりました。
ベル数による書きなおし
今は期待値を求めているので、先ほど求めた式の両辺をで割ります。またのときであるため、の部分をに書き換えることができます。これにより、求める値をと表せました。
ところで、をで足し合わせたものはベル数です。よってと書けます。特にの場合は、さえ求まればよいことになります。
Touchard's congruence
巨大なに対して、正確にはを求める必要があります。ここで素数が十分小さいことを利用します。
Touchard's congruenceと呼ばれる合同式が存在します。非負整数、素数についてという式ですが、これについても組み合わせ的な解釈を説明します。
まず、ベル数とは、個の区別できる玉をいくつかの箱に分ける方法でした。ここでを、において得られた分け方にさらに個の玉を追加したものとして捉えます。追加の際に箱を増やしても構いません。
いったん個の分け方を固定しましょう。ここに個の玉、説明のためと番号を付けますが、これを追加し、得られた個の玉の分け方を集めた集合をとします。最初に定めた個の分け方が異なれば当然そこから得られるの要素もすべて異なりますから、を考えてうまく足し合わせることでが得られます。
任意のについて、における玉を一斉にに置き換えたものもまたの要素になっています。これをとします。つまり、はの写像を与えます。また明らかにです。
を、を複数回適用することによって互いに移り変わる要素たちに分割します。巡回群によるへの作用の軌道を考えているとも言えますが、以降特に群の性質は使っていません。分割後の各部分集合のサイズを見ましょう。
であってであるようなものを取ります。と、が素数でが成り立つということから、がすべて異なると言えます。よって分割後の部分集合でが含まれるもののサイズはちょうどだとわかりました。今はを考えているので、この部分集合はの要素としては無視してよいです。
結局、なるを数える問題に帰着されました。ではこのはどんな構造をしているのでしょう。以下、をを満たすものとして固定します。
玉が入っていた箱に同時に入っていた玉の数を数え、とします。例えば玉に注目すると、個の玉がによって置き換えられ、のほうで数えると個になることから、が成り立ちます。ここでよりです。同様の議論でがすべて等しいことが言えました。つまり個の玉が同じ数ずつに分けられてそれぞれ箱に入っていることになりますが、は素数ですから、その分け方は「個ずつ」または「個ずつ」です。
- 個ずつ分ける場合
当然ながら、すでに分けられていた個のどれかと一緒に入ってしまうとが成り立たなくなるので、それぞれについて新しい箱を用意する通りしかありません。これを個の分け方である倍するのですから、に対する寄与もとなります。
- 個ずつ分ける場合
どの箱に入れても常にとなります。ここで個をまとめて1個の玉だと思うことで、個の玉を分ける方法と合わせ、にはだけ寄与します。
以上よりが言えました。ここまでの議論は参考文献[1]のLemma 3.1の証明に着想を得たものです。
また、より一般化した合同式として、非負整数についてが成り立ちますが、これについては帰納法で簡単に示すことができます。以下を法とします。
ただし最後の行への変形で、に対してであること、フェルマーの小定理よりであることを用いました。
この等式を使うことで、巨大なに対しても十分高速にを求めることができます。
を求める
さて、残るはです。ここにという制約が効いてきます。つまり、改めてと置きなおして、に対してをそれぞれ高速に求められれば良いのです。は自明なので取り除いておきます。
ところで、第2種オイラー数というものが存在します。これに関する日本語の文献はほとんどなく、英語で「Eulerian numbers of the second kind」などと検索する必要があります。三角かっこで2重にくくる表記が一般的ですが、ここでは第2種スターリング数(と参考文献[2])に合わせてで表記します。
この値の組み合わせ的な解釈を述べます。まず、を並べ替えた数列であって、「任意のについてとの間に出現する数はすべてより大きい」という条件を満たすものを、次のStirling permutationと呼びます。とは、次のStirling permutation であって、「を満たす()がちょうど個ある」という条件を満たすものの数です。
具体例など詳しくはWikipediaを参照してください。だと思う定義も存在し、文献によってはが少しずれていました。また左右が逆のものもあります。
このような数を突如導入した理由は2点、まず単純な漸化式が存在し小さい範囲の値を高速に求められるから、次にに対しが成り立つからです。
前者によってを前計算しておけば、後者によってを高速に求めることができます。が小さい今回であればリュカの定理によりが達成でき、が大きければ二項係数を徐々に更新することでが達成できます。
参考までに、が小さいときは一般のスターリング数を高速に求めることができるそうです。大きなについてを高速に求めたくなる機会は今後訪れるのでしょうか。
第2種オイラー数の漸化式
というものです。ただしまたはの場合としておきます。の場合にとなることは、以下の漸化式の導出により確かめることができます。
この漸化式はdpだと思うことで簡単に導けます。最初にの場合、条件を満たすただ一つの数列を見てと定めることができます。以降とします。
次のStirling permutationにを二つ挿入する方法を考えます。条件よりとの間に出現する数はすべてより大きい必要があるので、特に今、二つのは隣接している必要があります。挿入できる位置は両端を含めて箇所あります。
まず、でカウントされる並べ方について考えます。を挿入することでなるを一つ増やす必要がありますから、すでにそうなっているの後ろに挿入してはいけません。そのような位置は箇所あります。さらに左端に入れるとその左の項が存在しないため、これもまた意図に沿いません。以上より挿入できる位置は箇所なので、場合の数としてに対しだけ寄与します。
次にについて考えます。今度は逆で、先ほど挿入できなかった箇所にしか挿入できません。なるは今回個ありますから、その箇所と左端を合わせてだけ寄与します。
二つを合わせることで先の漸化式が得られました。
について
探したらこれにも組み合わせ的な解釈がありました。参考文献[2]のTheorem 2.3を用いるものです。ちなみにこれが、先ほど言及した左右が逆になっているものです。以下、文献と左右が逆になった状態のまま進めます。
まず、を、次のStirling permutationに対し、その項の間または両端に合計本仕切りを入れたものの集合とします。ただしこの時、数列の「左端」「なるの直後」には必ず1本以上仕切りが入っているものとします。
例えば、とし、仕切りをで表せば、のようになります。
実はそのような仕切りの入った数列と個の玉を個の区別できない箱に入れる方法が一対一対応しており、したがってが成り立ちます。
先に、ここから上の式を導いておきましょう。となるので、次のStirling permutationに対し、条件を満たしつつ本の仕切りを入れる方法を数え上げればよいです。
「なる」が(ただしかつ)個あったとすると、そのような数列は通り存在し、左端と合わせすでに本の仕切りの位置が確定しています。
残り本の仕切りを入れる位置を箇所から選びます。これは重複組み合わせによってと書け、二項係数に直すととなります。ここまで変形すると、先ほどのの制限を外してよくなります。仕切りが足りない場合を考慮していましたが、そのような場合からとなるからです。
まとめるととなり、見事求める式が得られました。
について
参考文献[2]では一方をもう一方に変換する方法を事細かく説明していました。これを引き写してくるのはさすがに面倒なので、具体例を見た後大まかな説明を試みます。
基本的には仕切り一つが箱一つ+玉一つに対応します。玉が二つ以上入っている箱については数列の同じ値を持つ要素二つを対応させます。先ほど挙げた例を変換してみましょう。これは、より、ラベルを個の集合に分割する方法の一つと対応しているはずです。
まず、を、と分解していきます。同じ値を持つ要素二つとその間の仕切りすべてを削除するのを繰り返しますが、このとき削除した要素のすぐ左には必ず別の仕切りが存在します。このことは、の定義で設けた仕切りの位置の制約から言えます。
削除されたのと逆順にラベルを振っていきます。まず最初の仕切りが二つでとなります。右から左の順でラベルが付いているのは、参考文献と第2種オイラー数の定義が逆だからです。次に、に対応した仕切りのすぐ右の要素を削除しましたから、それをラベルとしてに追加します。です。
以下同様に、と更新されていき、無事を個の集合に分割できました。どう表示したとしても見にくいので勘弁してください。
さて、分割を数列に戻しましょう。それぞれの集合のうち最小のラベルを仕切りに変換し、残りを倍にコピーして降順に並べると、となります。これを適切にマージします。具体的には、仕切りがもともと対応していたラベルを昇順に見て、それより小さいラベルが出現する最も右の位置のすぐ左に塊ごと挿入します。となりました。番号を順序を保って振りなおせばが得られます。
これはいったい何だったんでしょうか。数列の要素との入れ子関係によってラベルの番号を管理しているのは明らかです。このとき、すでに決まった集合の構造が失われないように要素を二つずつ使っていたのでしょう。
つまり、すべての仕切りについて、一つ右にある要素を見るのと同じ要素の次の出現位置に飛ぶのを交互に繰り返せば、スタート地点の仕切りと同じ集合に入っていたラベルが列挙できるのです。特に要素のラベルは単調減少に並ぶようにしましたが、これは仕切りの位置の制約から、となる場合その間に別の仕切りが必要であることに対応しています。これによって一意性が担保されています。
今説明した両方向の変換が、それぞれ互いの逆変換になっていることも、この理由付けから納得できそうです。僕が納得した気になったので終わります。
実装例
ACLからmodint
を、手持ちのライブラリからの解を復元するのにgarnerのアルゴリズムを使っています。
#include<iostream> #include<atcoder/modint> using namespace std; #include"math/garner.cpp" template<int p> struct A{ using mint=atcoder::static_modint<p>; mint C[p][p],Bell[p],B[5000][5000]; A() { for(int i=0;i<p;i++)for(int j=0;j<p;j++)C[i][j]=0; C[0][0]=1; for(int i=1;i<p;i++) { C[i][0]=1; for(int j=1;j<=i;j++)C[i][j]=C[i-1][j-1]+C[i-1][j]; } Bell[0]=1; for(int i=1;i<p;i++) { Bell[i]=0; for(int j=1;j<=i;j++)Bell[i]+=C[i-1][j-1]*Bell[j-1]; } for(int i=0;i<5000;i++)for(int j=0;j<5000;j++)B[i][j]=0; B[0][0]=1; for(int i=1;i<5000;i++) { B[i][0]=B[i-1][0]; for(int j=1;j<i;j++) { B[i][j]=(2*i-j-1)*B[i-1][j-1]+(j+1)*B[i-1][j]; } } } mint bell(long long k) { int m=0; long long q=1; while(q<=k/p)m++,q*=p; vector<mint>now; now.push_back(1); for(;m>=1;m--,q/=p) { int t=k/q; k%=q; vector<mint>F(t+1); mint coef=1; for(int i=t;i>=0;i--) { F[i]=coef*C[t][i]; coef*=m; } vector<mint>nxt(now.size()+F.size()-1); for(int i=0;i<now.size();i++)for(int j=0;j<F.size();j++)nxt[i+j]+=now[i]*F[j]; while(nxt.size()>p) { mint tmp=nxt.back(); nxt.pop_back(); nxt[nxt.size()-p]+=tmp; nxt[nxt.size()-p+1]+=tmp; } now=nxt; } mint ret=0; for(int i=0;i<now.size();i++) { int id=i+k; if(id<p)ret+=Bell[id]*now[i]; else ret+=(Bell[id-p]+Bell[id-p+1])*now[i]; } return ret; } mint combination(long long n,long long k) { if(n<k)return 0; mint ret=1; while(k>0) { ret*=C[n%p][k%p]; n/=p; k/=p; } return ret; } mint S2(long long k,long long m) { m=k-m; if(m==0)return 1; //S2(k,k-m) mint ret=0; for(int i=0;i<m;i++)ret+=B[m][i]*combination(k+m-i-1,2*m); return ret; } int solve(long long n,long long k) { if(k==0)return 1; mint ans=bell(k); for(long long m=n+1;m<=k;m++)ans-=S2(k,m); return ans.val(); } }; int main() { long long n,k; cin>>n>>k; vector<long long>x(3),m(3); x[0]=A<857>().solve(n,k);m[0]=857; x[1]=A<997>().solve(n,k);m[1]=997; x[2]=A<1009>().solve(n,k);m[2]=1009; cout<<garner(x,m)<<endl; }
参考文献
[1]:Greg Hurst, Andrew Schultz、2009、[0906.0696] An elementary (number theory) proof of Touchard's congruence
[2]:Ira Gessel, Richard P Stanley、1976、Stirling polynomials - ScienceDirect