"蔚来杯"2022牛客暑期多校训练营3 I Ice Drinking

"蔚来杯"2022牛客暑期多校训练营3_ACM/NOIP/NOI/CCPC/ICPC算法编程基础集训_牛客竞赛OJ

注(検索用):Multi-Universities Campとか、TwitterではNowCoderとか呼ばれている中国の合宿の問題です。

問題概要

長さnの順列P=(P_1,P_2,\dots,P_n)に対し、c(P)P_i=iなるi1\le i\le n)の個数として定めます。Pをランダムに決めたときのc(P)^kの期待値を\bmod 862118861で求めてください。

  • 1\le n\le 10^{18}
  • 0\le k\le n+5\times 10^3

解法

準備

まずk=0の場合は自明なので除きます。また素因数分解すると862118861=857\times 997\times 1009なので、素数p\approx 10^3について\bmod pで求めることができれば、中国剰余定理より答えが求まります。以降はこれらのことを前提にして解きます。

\displaystyle\sum_P c(P)^kを求める

解説スライドによれば、c(P)を固定して完全順列の式を使い式変形を頑張る方法もあるようですが、ここでは組み合わせ的解釈から導くことにします。

c(P){}_{c(P)}\mathrm{C}_1だと思うことにします。これはつまり、P_i=iであるようなiを一つ選ぶ場合の数です。今c(P)^kを求めているのですから、一つ取り出す操作をk回繰り返す方法を数え上げればよいです。このとき取り出す順番を考慮していることに注意します。

さて、k個取り出したiたちは当然大量に重複しますが、この重複を除いたものがm種類(ただし1\le m\le\min(n,k))だったとします。このmを固定してみます。

まず、m種類がどのように分布するかを考えます。k回の操作をそれぞれ区別される玉、m種類のiを区別しない箱だと思えば、求める分布の数はそれぞれの箱に1個以上の玉を入れる場合の数ですから、写像12相によればS(k,m)となります。ここでSは第2種スターリング数です。

次に、それを達成できる順列Pがどれくらいあるか考えます。1からnのうちm種類の箱として使うものを決めるのが、どの箱にどの数字を対応させるかまで含めて{}_n\mathrm{P}_m通りになります。先ほど箱を区別しなかった分、ここで区別する必要がありました。選んだm個の値は位置まで決定しており、残りは自由ですから(n-m)!通りで、掛け合わせると{}_n\mathrm{P}_m\times(n-m)!=n!になります。

順列としての重複を考える必要はありません。なぜなら今の考察は、最初「すべての順列について」「m種類の取り方すべてについて」の順で考えていたものを、逆順で表現しなおしただけだからです。

よって\sum_P c(P)^k=\sum_{m=1}^{\min(n,k)}S(k,m)n!であることがわかりました。

ベル数による書きなおし

今は期待値を求めているので、先ほど求めた式の両辺をn!で割ります。またk\lt mのときS(k,m)=0であるため、\min(n,k)の部分をnに書き換えることができます。これにより、求める値を\sum_{m=1}^n S(k,m)と表せました。

ところで、S(k,m)m=1\dots kで足し合わせたものはベル数B_kです。よって\sum_{m=1}^n S(k,m)=B_k-\sum_{m=n+1}^k S(k,m)と書けます。特にk\le nの場合は、B_kさえ求まればよいことになります。

Touchard's congruence

巨大なkに対してB_k、正確にはB_k\bmod pを求める必要があります。ここで素数pが十分小さいことを利用します。

Touchard's congruenceと呼ばれる合同式が存在します。非負整数n素数pについてB_{n+p}\equiv B_n+B_{n+1}\pmod pという式ですが、これについても組み合わせ的な解釈を説明します。

まず、ベル数B_nとは、n個の区別できる玉をいくつかの箱に分ける方法でした。ここでB_{n+p}を、B_nにおいて得られた分け方にさらにp個の玉を追加したものとして捉えます。追加の際に箱を増やしても構いません。

いったんn個の分け方を固定しましょう。ここにp個の玉、説明のため0,1,\dots,p-1と番号を付けますが、これを追加し、得られたn+p個の玉の分け方を集めた集合をAとします。最初に定めたn個の分け方が異なれば当然そこから得られるAの要素もすべて異なりますから、|A|\bmod pを考えてうまく足し合わせることでB_{n+p}\bmod pが得られます。

任意のa\in Aについて、aにおける玉0,1,\dots,p-2,p-1を一斉に1,2,\dots,p-1,0に置き換えたものもまたAの要素になっています。これをf(a)とします。つまり、fA\rightarrow A写像を与えます。また明らかにf^p=\mathrm{id}_Aです。

Aを、fを複数回適用することによって互いに移り変わる要素たちに分割します。巡回群\langle f\rangleによるAへの作用の軌道を考えているとも言えますが、以降特に群の性質は使っていません。分割後の各部分集合のサイズを見ましょう。

a\in Aであってf(a)\ne aであるようなものを取ります。f(a)\ne aと、p素数f^p(a)=aが成り立つということから、a,f(a),f^2(a),\dots,f^{p-1}(a)がすべて異なると言えます。よって分割後の部分集合でaが含まれるもののサイズはちょうどpだとわかりました。今は|A|\bmod pを考えているので、この部分集合はAの要素としては無視してよいです。

結局、f(a)=aなるa\in Aを数える問題に帰着されました。ではこのaはどんな構造をしているのでしょう。以下、a\in Af(a)=aを満たすものとして固定します。

iが入っていた箱に同時に入っていた玉0,1,\dots,p-1の数を数え、c_a(i)とします。例えば玉0に注目すると、c_a(0)個の玉がfによって置き換えられ、f(a)のほうで数えるとc_{f(a)}(1)個になることから、c_a(0)=c_{f(a)}(1)が成り立ちます。ここでf(a)=aよりc_a(0)=c_a(1)です。同様の議論でc_a(i)がすべて等しいことが言えました。つまりp個の玉が同じ数ずつに分けられてそれぞれ箱に入っていることになりますが、p素数ですから、その分け方は「1個ずつ」または「p個ずつ」です。

  • 1個ずつ分ける場合

当然ながら、すでに分けられていたn個のどれかと一緒に入ってしまうとf(a)=aが成り立たなくなるので、それぞれについて新しい箱を用意する1通りしかありません。これをn個の分け方であるB_n倍するのですから、B_{n+p}に対する寄与もB_nとなります。

  • p個ずつ分ける場合

どの箱に入れても常にf(a)=aとなります。ここでp個をまとめて1個の玉だと思うことで、n個の玉を分ける方法と合わせ、B_{n+p}にはB_{n+1}だけ寄与します。

以上よりB_{n+p}\equiv B_n+B_{n+1}\pmod pが言えました。ここまでの議論は参考文献[1]のLemma 3.1の証明に着想を得たものです。

また、より一般化した合同式として、非負整数mについてB_{n+p^m}\equiv mB_n+B_{n+1}\pmod pが成り立ちますが、これについては帰納法で簡単に示すことができます。以下pを法とします。

\displaystyle\begin{eqnarray}B_{n+p^{l+1}}&\equiv&lB_{n+p^{l+1}-p^l}+B_{n+p^{l+1}-p^l+1}\\
&\equiv&l^2B_{n+p^{l+1}-2\times p^l}+2lB_{n+p^{l+1}-2\times p^l+1}+B_{n+p^{l+1}-2\times p^l+2}\\
&\vdots&\\
&\equiv&\sum_{i=0}^p {}_p\mathrm{C}_i l^{p-i}B_{n+p^{l+1}-p\times p^l+i}\\
&\equiv&l^p B_n+B_{n+p}\equiv lB_n+B_{n+p}\equiv(l+1)B_n+B_{n+1}\end{eqnarray}

ただし最後の行への変形で、0\lt i\lt pに対して{}_p\mathrm{C}_i\equiv 0であること、フェルマーの小定理よりl^p\equiv lであることを用いました。

この等式を使うことで、巨大なkに対しても十分高速にB_k\bmod pを求めることができます。

S(k,m)を求める

さて、残るは\sum_{m=n+1}^k S(k,m)です。ここにk\le n+5\times 10^3という制約が効いてきます。つまり、改めてm\leftarrow k-mと置きなおして、0\le m\lt k-n\le 5\times 10^3に対してS(k,k-m)をそれぞれ高速に求められれば良いのです。m=0は自明なので取り除いておきます。

ところで、第2種オイラー数というものが存在します。これに関する日本語の文献はほとんどなく、英語で「Eulerian numbers of the second kind」などと検索する必要があります。三角かっこで2重にくくる表記が一般的ですが、ここでは第2種スターリング数(と参考文献[2])に合わせてB(n,k)で表記します。

この値の組み合わせ的な解釈を述べます。まず、1,1,2,2,\dots,n,nを並べ替えた数列であって、「任意のiについてiiの間に出現する数はすべてiより大きい」という条件を満たすものを、n次のStirling permutationと呼びます。B(n,k)とは、n次のStirling permutation a=(a_1,a_2,\dots,a_{2n})であって、「a_j\lt a_{j+1}を満たすj1\le j\lt 2n)がちょうどk個ある」という条件を満たすものの数です。

具体例など詳しくはWikipediaを参照してください。a_0=0だと思う定義も存在し、文献によってはkが少しずれていました。また左右が逆のものもあります。

このような数を突如導入した理由は2点、まず単純な漸化式が存在し小さい範囲の値を高速に求められるから、次にm\ge 1に対しS(k,k-m)=\sum_{i=0}^{m-1}B(m,i){}_{k+m-i-1}\mathrm{C}_{2m}が成り立つからです。

前者によってB(m,i)を前計算しておけば、後者によってS(k,k-m)\bmod pを高速に求めることができます。pが小さい今回であればリュカの定理によりO(m\log m)が達成でき、pが大きければ二項係数を徐々に更新することでO(m\log p)が達成できます。

参考までに、pが小さいときは一般のスターリング数を高速に求めることができるそうです。大きなpについてS(k,k-m)\bmod pを高速に求めたくなる機会は今後訪れるのでしょうか。

maspypy.com

第2種オイラー数の漸化式

B(n,k)=(2n-k-1)B(n-1,k-1)+(k+1)B(n-1,k)というものです。ただしk\lt 0またはn\le kの場合B(n,k)=0としておきます。n\le kの場合にB(n,k)となることは、以下の漸化式の導出により確かめることができます。

この漸化式はdpだと思うことで簡単に導けます。最初にn=0の場合、条件を満たすただ一つの数列(1,1)を見てB(1,0)=1と定めることができます。以降n\ge 2とします。

n-1次のStirling permutationにnを二つ挿入する方法を考えます。条件よりnnの間に出現する数はすべてnより大きい必要があるので、特に今、二つのnは隣接している必要があります。挿入できる位置は両端を含めて2(n-1)+1=2n-1箇所あります。

まず、B(n-1,k-1)でカウントされる並べ方について考えます。nを挿入することでa_j\lt a_{j+1}なるjを一つ増やす必要がありますから、すでにそうなっているjの後ろに挿入してはいけません。そのような位置はk-1箇所あります。さらに左端に入れるとその左の項が存在しないため、これもまた意図に沿いません。以上より挿入できる位置は2n-1-(k-1)-1=2n-k-1箇所なので、場合の数としてB(n,k)に対し(2n-k-1)B(n-1,k-1)だけ寄与します。

次にB(n-1,k)について考えます。今度は逆で、先ほど挿入できなかった箇所にしか挿入できません。a_j\lt a_{j+1}なるjは今回k個ありますから、そのk箇所と左端を合わせて(k+1)B(n-1,k)だけ寄与します。

二つを合わせることで先の漸化式が得られました。

\displaystyle S(k,k-m)=\sum_{i=0}^{m-1}B(m,i){}_{k+m-i-1}\mathrm{C}_{2m}について

探したらこれにも組み合わせ的な解釈がありました。参考文献[2]のTheorem 2.3を用いるものです。ちなみにこれが、先ほど言及した左右が逆になっているものです。以下、文献と左右が逆になった状態のまま進めます。

まず、P_{k,n}を、k次のStirling permutationに対し、その項の間または両端に合計n本仕切りを入れたものの集合とします。ただしこの時、数列の「左端」「a_j\lt a_{j+1}なるjの直後」には必ず1本以上仕切りが入っているものとします。

例えばk=3n=4とし、仕切りを/で表せば、/3/31/221/のようになります。

実はそのような仕切りの入った数列とn+k個の玉をn個の区別できない箱に入れる方法が一対一対応しており、したがって|P_{k,n}|=S(n+k,n)が成り立ちます。

先に、ここから上の式を導いておきましょう。S(k,k-m)=|P_{m,k-m}|となるので、m次のStirling permutationに対し、条件を満たしつつk-m本の仕切りを入れる方法を数え上げればよいです。

a_j\lt a_{j+1}なるj」がi(ただし0\le i\lt mかつi\lt k-m)個あったとすると、そのような数列はB(m,i)通り存在し、左端と合わせすでにi+1本の仕切りの位置が確定しています。

残りk-m-i-1本の仕切りを入れる位置を2m+1箇所から選びます。これは重複組み合わせによって{}_{2m+1}\mathrm{H}_{k-m-i-1}と書け、二項係数に直すと{}_{(2m+1)+(k-m-i-1)-1}\mathrm{C}_{k-m-i-1}={}_{k+m-i-1}\mathrm{C}_{2m}となります。ここまで変形すると、先ほどのi\lt k-mの制限を外してよくなります。仕切りが足りない場合を考慮していましたが、そのような場合k+m-i-1\lt 2mから{}_{k+m-i-1}\mathrm{C}_{2m}=0となるからです。

まとめるとS(k,k-m)=\sum_{i=0}^{m-1}B(m,i){}_{k+m-i-1}\mathrm{C}_{2m}となり、見事求める式が得られました。

\displaystyle|P_{k,n}|=S(n+k,n)について

参考文献[2]では一方をもう一方に変換する方法を事細かく説明していました。これを引き写してくるのはさすがに面倒なので、具体例を見た後大まかな説明を試みます。

基本的には仕切り一つが箱一つ+玉一つに対応します。玉が二つ以上入っている箱については数列の同じ値を持つ要素二つを対応させます。先ほど挙げた例/3/31/221/を変換してみましょう。これはn+k=7n=4より、ラベル1,2,\dots,74個の集合に分割する方法の一つと対応しているはずです。

まず、/3/31/221/を、/3/31/221/\rightarrow/1/221/\rightarrow/1/1/\rightarrow//と分解していきます。同じ値を持つ要素二つとその間の仕切りすべてを削除するのを繰り返しますが、このとき削除した要素のすぐ左には必ず別の仕切りが存在します。このことは、P_{k,n}の定義で設けた仕切りの位置の制約から言えます。

削除されたのと逆順にラベルを振っていきます。まず最初の仕切りが二つで/_2/_1\leftrightarrow\{1\},\{2\}となります。右から左の順でラベルが付いているのは、参考文献と第2種オイラー数の定義が逆だからです。次に、2に対応した仕切りのすぐ右の要素を削除しましたから、それをラベル3として\{2\}に追加します。/_21_3/_41/_1\leftrightarrow\{1\},\{2,3\},\{4\}です。

以下同様に/_21_3/_42_521/_1\leftrightarrow\{1\},\{2,3\},\{4,5\}/_23_6/_731_3/_42_521/_1\leftrightarrow\{1\},\{2,3,6\},\{4,5\},\{7\}と更新されていき、無事1,2,\dots,74個の集合に分割できました。どう表示したとしても見にくいので勘弁してください。

さて、分割を数列に戻しましょう。それぞれの集合のうち最小のラベルを仕切りに変換し、残りを倍にコピーして降順に並べると、/_1,/_26633,/_455,/_7となります。これを適切にマージします。具体的には、仕切りがもともと対応していたラベルを昇順に見て、それより1小さいラベルが出現する最も右の位置のすぐ左に塊ごと挿入します。/_1\rightarrow/_26633/_1\rightarrow/_2663/_4553/_1\rightarrow/_26/_763/_4553/_1となりました。番号を順序を保って振りなおせば/3/31/221/が得られます。

これはいったい何だったんでしょうか。数列の要素との入れ子関係によってラベルの番号を管理しているのは明らかです。このとき、すでに決まった集合の構造が失われないように要素を二つずつ使っていたのでしょう。

つまり、すべての仕切りについて、一つ右にある要素を見るのと同じ要素の次の出現位置に飛ぶのを交互に繰り返せば、スタート地点の仕切りと同じ集合に入っていたラベルが列挙できるのです。特に要素のラベルは単調減少に並ぶようにしましたが、これは仕切りの位置の制約から、a_j\lt a_{j+1}となる場合その間に別の仕切りが必要であることに対応しています。これによって一意性が担保されています。

今説明した両方向の変換が、それぞれ互いの逆変換になっていることも、この理由付けから納得できそうです。僕が納得した気になったので終わります。

実装例

ACLからmodintを、手持ちのライブラリから\bmod 862118861の解を復元するのに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