"蔚来杯"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$ なる $i$ （ $1\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)$ とします。つまり、 $f$ は $A\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 A$ を $f(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$ について $i$ と $i$ の間に出現する数はすべて $i$ より大きい」という条件を満たすものを、 $n$ 次のStirling permutationと呼びます。 $B(n,k)$ とは、 $n$ 次のStirling permutation $a=(a_1,a_2,\dots,a_{2n})$ であって、「 $a_j\lt a_{j+1}$ を満たす $j$ （ $1\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$ を二つ挿入する方法を考えます。条件より $n$ と $n$ の間に出現する数はすべて $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=3$ 、 $n=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=7$ 、 $n=4$ より、ラベル $1,2,\dots,7$ を $4$ 個の集合に分割する方法の一つと対応しているはずです。

まず、 $/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,7$ を $4$ 個の集合に分割できました。どう表示したとしても見にくいので勘弁してください。

さて、分割を数列に戻しましょう。それぞれの集合のうち最小のラベルを仕切りに変換し、残りを倍にコピーして降順に並べると、 $/_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

kotatsugameの日記

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