ARC104 E - Random LIS 解説 - Motsu_xe 競プロのhogeやfuga

ちょっとだけ面白い解き方をしたので解説として残しておく。
E - Random LIS

問題概要

長さ $N\leq6$ の数列 $X$ は、 $i$ 番目の要素が $1\leq X_i\leq A_i$ を満たす整数として一様分布で与えられる。最長増加部分列の長さの期待値を求めよ。

解法

LISは $X_i$ の大小関係で完全に決まるから、大小関係が一致しているような数列を数え上げられれば十分。

大小関係の型の取り方は、まず一致しているかどうかを無視して $N!$ 通りの順番で並べ、隣接する2つが一致しているか否かの $2^{N-1}$ 通りを考えば良い。この際被ったもの(例えば $1\lt 2=3,1\lt 3=2$ は同じ関係である)は適当に無視すれば良い。

あとは大小関係の型が一致する数列を数え上げれば良い。このとき一致している項に関しては $A_i$ の最小を上限として選べるので、長さ $L\leq N$ の数列 $Y$ が、 $i$ 番目の要素として $1\leq Y_i\leq B_i$ の範囲を全て動くとき、単調増加列の数を数え上げるという問題に言い換えられた(ここまでがわからなかったら公式解説とか読んで)。

「初項が $x$ 以上であるような長さ $L$ の数列 $Y$ が、 $i$ 番目の要素として $1\leq Y_i\leq B_i$ の範囲を全て動くとき、単調増加列の数」を $f_Y(x)$ と定義する。ただし与えられる $B$ は単調増加として良い。*1このとき求めるべきは $f_Y(1)$ となる。ここで $Y$ から初項を除いた数列 $Y'$ を考えると
$\begin{align} \qquad f_Y(x)=\sum_{z=x}^{B_0}f_{Y'}(z+1) \end{align}$
となることに注意する。ここで上昇階乗冪というものを導入する。
$\begin{align} \qquad x^{(n)}=\prod_{i=0}^{n-1}(x+i)=x(x+1)...(x+n-1) \end{align}$
このとき、重要な性質として
$\begin{align} \qquad \sum_{x=1}^{X} x^{(n)}=\frac1{n+1}X^{(n+1)} \end{align}$
が成立する(証明は簡単)。これを利用して実際に $f_Y$ を求める。 $f_Y$ が $l$ 次の上昇階乗冪多項式で表せることを帰納的に示す。

$l=0$ のとき、数列は空なので $1$ 通りであるから $f_Y(l)=1$ となり条件を満たす。

$l>0$ のとき、長さ $l-1$ で成立しているとして $f_{Y'}(x)=\sum_{i=0}^{l-1}a_ix^{(i)}$ と書けたとする。このとき
$\begin{align} \qquad f_Y(x)&=\sum_{z=x}^{B_0}f_{Y'}(z+1)\\ &=\sum_{z=x+1}^{B_0+1}f_{Y'}(z)\\ &=\sum_{z=x+1}^{B_0+1}\sum_{i=0}^{l-1}a_iz^{(i)}\\ &=\sum_{i=0}^{l-1}\sum_{z=x+1}^{B_0+1}a_iz^{(i)}\\ &=\sum_{i=0}^{l-1}a_i\left(\sum_{z=1}^{B_0+1}z^{(i)}-\sum_{z=1}^{x}z^{(i)}\right)\\ &=\sum_{i=0}^{l-1}a_i\left(\frac{1}{i+1}(B_0+1)^{(i+1)}-\frac{1}{i+1}x^{(i+1)}\right)\\ &=\sum_{i=0}^{l-1}a_i\left(\frac{1}{i+1}(B_0+1)^{(i+1)}-\frac{1}{i+1}x^{(i+1)}\right)\\ &=\sum_{i=1}^l\frac{a_{i-1}}{i}(B_0+1)^{(i)}-\sum_{i=1}^l\frac{a_{i-1}}{i}x^{(i)}\\ \end{align}$
となり、示された。

上昇階乗冪多項式で表せることを示したついでに、具体的な計算方法もわかったので、あとは計算するだけである。

最後に計算量を確認する。 $l$ 次上昇階乗冪多項式の積分操作は素朴にやって $O(l)$ *2、代入操作は素朴にやると $O(l^2)$ かかるが、上昇階乗冪を下位の次数から順番に求めれば $O(l)$ で実行できる。よって再帰1回は $O(l)$ で実行されているから、全体で $O(l^2)$ で計算される。Ordered Bell numberを $F_N$ とすれば、全体の計算量は $O(N!2^N\log F_N+N^2F_N)$ となり、 $F_N\leq N!2^{(N-1)}$ に注意して十分高速である。

実装例

#include<bits/stdc++.h>
#include<atcoder/all>
using namespace std;
using namespace atcoder;
using ll=long long;
template<class T,class U> inline bool chmin(T&x,U y){if(x>y){x=y;return true;}return false;}

template <class T> int LIS(vector<T>&a){
    int n=a.size();
    vector<T> l(n+1,numeric_limits<T>::max());
    for(int i{};i<(n);++i){
        *lower_bound(l.begin(),l.end(),a[i])=a[i];
    }
    for(int i{1};i<=(n);++i){
        if(l[i]==numeric_limits<T>::max()) return i;
    }
}

using mint=modint1000000007;

mint substitution(const vector<mint>&f,mint x){
    mint ans{},X{1};
    int n=f.size();
    for(int i{};i<(n);++i){
        ans+=f[i]*X;
        X*=x;
        ++x;
    }
    return ans;
}

vector<mint> rec(vector<int>&a){
    if(a.empty()) return {1};
    int A=a.back();
    a.pop_back();
    auto D=rec(a);
    int n=D.size();
    vector<mint> I(n+1);
    for(int i{1};i<=(n);++i){
        I[i]=-D[i-1]/i;
    }
    I[0]=-substitution(I,A+1);
    return I;
}

void solve(){
    unsigned n;
    cin>>n;
    vector<int> a(n);
    for(auto&i:a) cin>>i;
    vector<int> p(n);
    iota(p.begin(),p.end(),0);
    set<vector<vector<int>>> st;
    do{
        for(unsigned j{};j<(1u<<(n-1));++j){
            vector<vector<int>> v(1);
            unsigned b=j;
            for(auto&i:p){
                v.back().push_back(i);
                if(b&1u) v.emplace_back();
                b>>=1u;
            }
            for(auto&u:v) sort(u.begin(),u.end());
            st.insert(v);
        }
    }while(next_permutation(p.begin(),p.end()));
    mint ans{};
    for(auto&v:st){
        int N=v.size();
        vector<int> T(n),A(N,INT_MAX);
        for(int i{};i<(N);++i){
            for(auto&j:v[i]){
                T[j]=i;
                chmin(A[N-1-i],a[j]);
            }
        }
        for(int i{1};i<(N);++i){
            chmin(A[i],A[i-1]);
        }
        mint P=substitution(rec(A),1);
        int L=LIS(T);
        ans+=P*L;
    }
    for(int i{};i<(n);++i){
        ans/=a[i];
    }
    cout<<ans.val()<<'\n';
}

int main(){
    cin.tie(nullptr);
    ios::sync_with_stdio(false);
    solve();
}