みなさんは†AtCoder最難†と呼ばれる伝説の問題、「e」をご存じですか？
AtCoderのレーティング上位たった8人が集まり、究極の難しさの問題6問で競う、「World Tour Finals 2019」の最終問題です。
2020年3月13日現在もたった17人にしか解かれておらず、現在AtCoder史上最難といわれる問題です。
ここまで正解者数が少ない原因として、公式解説で詳細な説明がされておらず、非公式の解説ブログも非常に少ない（現在1件しかないと思われます）ことがあります。
今回は、ここでeの解説を執筆することで、eという問題に少しでも多くの人が挑戦できるようになることを願っています。
問題リンクはこちらです。

前提知識

• 数Ⅲの積分範囲（基礎のみで大丈夫です）
• テイラー展開・マクローリン展開の初歩的な知識

これだけで理解はできると思います。他になにか不足している前提知識があるようでしたら、コメント欄に質問してもらえると助かります。

問題概要

1$1$ 列に並んだ座席が M$M$ 個並んでいる。
ここに M$M$ 人の人が順に訪れ、次のような行動をする。

• 両隣に人がいない席があれば、そのうち1つをランダムに選んで座る。
• そのような席がなければ、どこにも座らないで立ち去る。

これがすべて終了した後、M$M$ 個の座席のうち連続する N$N$ 個を選んで写真を撮る。
このとき、この N$N$ 個の座席の人の有無が、与えられた文字列と同じようになる確率を考える。
この確率は M$M$ に依存するが、M$M$ が限りなく大きくなる時の確率の極限を求めよ。
ただし、この確率は、有理数 p,q,r$p,q,r$ を用いて

の形で表せることが証明できるので、この p,q,r$p,q,r$ を mod109+7$mod{10}^9+7$ で出力せよ。

入力例1

文字列の長さ 1$1$,文字列が X のとき、確率は

らしい。
もう、入力例の時点からどうしたらこれが出てくるのか謎ですね。
これが最難といわれる所以でもあります。
幸い、これに関する解説は公式の解説PDFにもあるので、理解はそこまで難しくないはずです。
ここでその証明をしてみましょう。

まず、座席の決め方を少し言い換えます。
0$0$ から 1$1$ までの任意の実数を、重複しないように、また、完全にランダムにすべての座席に割り当てます。i$i$ 番目の席の実数を Xi$X_i$ とします。
そして、座りに来た人はこの実数の小さい順に座席を見ていって、もし座れる席があればそこに座る、とします。
選ぶ実数のランダム性を考えれば、こう言い換えても結果は変わりません。
このとき、i$i$ 番目の座席に人が座る条件はなんでしょうか。
これは、i$i$ 番目から左右へ、 Xj$X_j$ が単調減少となる最長の部分列の長さがどちらも偶数となることです。
例えば、Xi−2>Xi−1<Xi$X_{i-2}>X_{i-1}<X_i$ のとき、Xi−1$X_{i-1}$ が先に埋まるため、Xi$X_i$ に人が座ることはありません。
一方、Xi−3>Xi−2<Xi−1<Xi$X_{i-3}>X_{i-2}<X_{i-1}<X_i$ のとき、Xi−2,Xi$X_{i-2},X_i$ が埋まります。
このように、単調な列では一つ飛ばしで席の埋まり方が決まっていくため、左右どちらからも Xi$X_i$ の隣が埋まらないのが、Xi$X_i$ に人が座る必要十分条件です。
もちろん逆に、座らない条件は左右どちらかが奇数となることです。
これを利用して、最終的にある席が埋まる確率を計算してみましょう。
P($P($不等式)$)$を、不等式が成り立つ確率とすると、Xi$X_i$ の左側の、最長の単調減少列の長さが偶数になる確率 E$E$ は、Xi=x$X_i=x$ のとき、包除原理を考えて、

e−x$e^{-x}$ のマクローリン展開を考えれば、この変形ができます。
これが左右で成立しないといけないので、左右合わせた確率は e−2x$e^{-2x}$ となります。
x$x$ の値は、閉区間 [0,1]$[0,1]$ の中からランダムに決まるので、答えは積分して求められます。

これで、入力例1の証明はできましたね。

解法の概要

入力例1での考察を発展させて、問題そのものを解くことを考えましょう。
左から順に Xi$X_i$ の値を決定していくことを考え、DPを行います。
ただし、もちろん文字列に状況が合致する確率は Xi$X_i$ に依存するので、Xi$X_i$ の関数として保持する必要があります。この方法は後に説明します。
ひとまず、DPの遷移式と初期値、答えの計算について考えましょう。

DPの定義

dp[i][j][k][l]$dp[i][j][k][l]$ を、「i$i$ 番目に割り当てた値から、文字列のうち閉区間 [0,i]$[0,i]$ が状況に合致し、かつ左側への最長の単調減少部分列の長さを 2$2$ で割った余りが j$j$ になる確率を求める関数」とします。
ただし、この時右側への最長の単調減少部分列の長さも当然ここには関わってくるので、k,l$k,l$ で、右側への最長単調減少部分列の長さが偶数の時成り立つか、また、奇数の時成り立つかも持っておきます。
実はこれで遷移式が立てられるのです。次の項で詳しく考えましょう。

DPの初期値

ここはそこまで難しくないでしょう。

s[0]='X'の場合

添字については注意深く考えましょう。
左右ともに列の長さは偶数である必要があるので、

です。

s[0]='-'の場合

左右どちらかの列の長さが奇数である必要があるので、

です。

DPの遷移式

以下では煩雑さを避けるため、添え字が 0$0$ のときと 1$1$ のとき両方に対して適用されるような場合、[0/1]$[0/1]$ のように表記します。
場合分けして考えていきます。

s[i]='X' の場合

左右どちらの列の長さも偶数である必要があるので、遷移先は dp[i][0][1][0]$dp[i][0][1][0]$ のみです。
ここでさらに場合分けをします。
Xi−1<Xi$X_{i-1}<X_i$ のとき、左側が偶数であるには、その一つ前では左側が奇数である必要があります。
また、遷移前の右側は自動的に偶数になります。
この関数が区間 [0,Xi)$[0,X_i)$ でどのような値をとるかについて考えることになるので、遷移式は次のようになります。

Xi−1>Xi$X_{i-1}>X_i$ のとき、遷移後の左側は自動的に偶数になるので、遷移前の左側はどちらでもよいです。
遷移後の右側が偶数であることから、遷移前の右側は奇数でないといけないので、上と同様に、

s[i]='-' の場合

今度は逆に、左右どちらかの列の長さを奇数にしなければなりません。
このとき遷移先として考えられるのは、dp[i][0][0][1]$dp[i][0][0][1]$ と dp[i][1][1][1]$dp[i][1][1][1]$ の2つです。
前と同様な場合分けをします。
Xi−1<Xi$X_{i-1}<X_i$ のとき、さらに2パターンに分かれます。
遷移後に左が偶数なら、遷移前の左は奇数で、遷移後の右は奇数でなければいけません。
遷移前の右は自動的に偶数になります。

遷移後に左が奇数なら、右はどちらでもよいです。
遷移前の左は偶数であり、右は自動的に偶数になります。

Xi−1>Xi$X_{i-1}>X_i$ のときは、遷移後の左は自動的に偶数になり、遷移前の右は偶数、遷移後の右は奇数である必要があります。遷移前の左はどちらでもよいです。

これで遷移式の計算は終わりです。お疲れさまでした。下にすべての遷移式をまとめておきます。

答えの計算

さて、遷移式と初期値が分かったのでDPの計算はできますが、答えが直接求められるわけではありません。
DPの計算が終わった後、最終的な答えは、座席のさらに右側に続く列の状況がどうなるかも加味して考えなければいけません。
文字列の一番右に割り当てる実数を [0,1]$[0,1]$ の間からとり、それより右に続く最長の単調減少部分列の長さが、偶数でも奇数でも構わない場合は 1$1$ 倍、偶数でなければいけない場合は e−x$e^{-x}$ 倍、奇数でなければいけない場合は 1−e−x$1-e^{-x}$ 倍するので、最終的な答え A$A$ は、次のように計算できます。

あとは実装方針について考えていきましょう。

関数の保持

実は、DPの中で出てくる関数は、すべて N−1$N-1$ 次以下の有理数係数多項式 p,q$p,q$ と有理数定数 c$c$ を用いて、

と表せます。
これを積分してみればなぜこうなるかわかると思います。
まず、[0,x]$[0,x]$ の積分区間で積分してみましょう。

ふう。見ていただければわかるとおり、積分の前と同じ構造になり、中括弧の中の多項式を一つにまとめれば全く同じように処理できますね。
では、次、[x,1]$[x,1]$ の積分区間でいきましょう。

疲れましたね。というのも、TEXモードでこの数式を書き上げるのは結構疲れるのです.
この計算で、2種類の積分区間でどれだけ積分を繰り返しても関数はこの形となることがわかりました。
つまり、DPの際には、2 つの多項式と、ひとつの定数だけを保持すればよいのです。
また、1度の遷移で高々1回の積分をするだけなので、多項式の次数は最終的にも高々 N$N$ 次で済みます。
これで、多項式の積分に O(N)$O(N)$,DPの遷移に O(N)$O(N)$ で、全体で O(N2)$O(N^2)$ でDPの遷移までが完了しました。

あと一息

ここまでで終わった気になっていたそこのあなた。ここからももう1段階計算が要ります。
DP配列から答えを計算するには、次の積分が必要になります。

実は、多項式と e−x$e^{-x}$ の積の積分は結構難しいです。
もちろん部分積分法を用いて x$x$ の次数を落としていけば漸化式で答えが求められますが、初等的でもっと簡単な導出方法が、参考にした他の解説に記載されていたのでそれを使います。
0<|t|<1$0<|t|<1$ に対し、

とする。
積分してべき無限和（等比数列の無限和）に直すと、

一方で、テイラー展開を考えると、I$I$ は次のようにも表せる。

テイラー展開の一意性から、両辺を結んで、

結局、

となる。
これを使って、多項式と e−x$e^{-x}$ の積の積分は当然 O(N2)$O(N^2)$ で求められるので、DPパートも合わせて全体で O(N2)$O(N^2)$ の計算量が達成できます。
定数倍は重めですが、実行時間制限は長いので余裕をもって通ると考えられます。
お疲れさまでした。

実装

さて、この問題は実装も結構重いです。
かなりコードは汚いかもしれませんが、自分の実装例を示しておきます。

マクロなど

使ったマクロなどを置いておきます。

#define rep(i,n) for(int i=0;i<(n);i++)
#define REP(i,n) for(int i=1;i<=(n);i++)
#define ll long long

剰余類クラス

これは他の問題でもよく使いますが、剰余類を扱うためのクラスを持っておきましょう。もちろん、有理数も法を定義した上で問題なく扱えます。
ただ、法が素数でない場合除算などができないので、注意しましょう。

constexpr int mod = 100000007;
class modInt {
    ll value;
public:
    modInt() : value(0) {}
    template<typename T>
    modInt(T value = 0) : value(value) {
        if (value < 0)value = -(-value % mod) + mod;
        this->value = value % mod;
    }
    inline operator int()const { return value; }
    inline modInt& operator+=(const modInt& x) {
        value += x.value;
        if (value >= mod)value -= mod;
        return *this;
    }
    inline modInt& operator++() {
        if (value == mod - 1)value = 0;
        else value++;
        return *this;
    }
    inline modInt operator-()const {
        return modInt(0) -= *this;
    }
    inline modInt& operator-=(const modInt& x) {
        value -= x.value;
        if (value < 0)value += mod;
        return *this;
    }
    inline modInt& operator--() {
        if (value == 0)value = mod - 1;
        else value--;
        return *this;
    }
    inline modInt& operator*=(const modInt& x) {
        value = value * x.value % mod;
        return *this;
    }
    inline modInt& operator/=(modInt rhs) {
        int exp = mod - 2;
        while (exp) {
            if (exp & 1)*this *= rhs;
            rhs *= rhs;
            exp >>= 1;
        }
        return *this;
    }
    template<typename T> modInt operator+(T x)const { return modInt(*this) += x; }
    template<typename T> modInt& operator+=(T x) { return operator+=(modInt(x)); }
    template<typename T> modInt operator-(T x)const { return modInt(*this) -= x; }
    template<typename T> modInt& operator-=(T x) { return operator-=(modInt(x)); }
    template<typename T> modInt operator*(T x)const { return modInt(*this) *= x; }
    template<typename T> modInt& operator*=(T x) { return operator*=(modInt(x)); }
    template<typename T> modInt operator/(T x)const { return modInt(*this) /= x; }
    template<typename T> modInt& operator/=(T x) { return operator/=(modInt(x)); }
};
istream& operator>>(istream& ist, modInt& x) {
    ll a;
    ist >> a;
    x = a;
    return ist;
}
modInt modpow(modInt a, int b) {
    if (!b)return modInt(1);
    if (b & 1)return modpow(a, b - 1) * a;
    modInt memo = modpow(a, b / 2);
    return memo * memo;
}

多項式クラス

まずは有理数係数多項式クラスを用意しましょう。積分をする機能もつけておくと便利です。
外から項をいじりやすくするためにメンバは全部publicにしています（本当はメソッドを通すなどするべき）
足し算の機能だけつけておきます。

class Polynomial {
public:
    vector<modInt> vec;
    Polynomial() {}
    Polynomial(const vector<modInt>& vec) :vec(vec) {}
    modInt get(const modInt& x) {
        modInt res = 0;
        rep(i, vec.size()) {
            res += modpow(x, i) * vec[i];
        }
        return res;
    }
    Polynomial integrate()const {
        vector<modInt> res = vec;
        res.emplace_back();
        for (int i = (int)res.size() - 2; i >= 0; i--)res[i + 1] = res[i] / (i + 1);
        res[0] = 0;
        return Polynomial(res);
    }
    Polynomial& operator+=(const Polynomial& rhs) {
        vec.resize(max(vec.size(), rhs.vec.size()));
        rep(i, rhs.vec.size())vec[i] += rhs.vec[i];
        return *this;
    }
};

関数クラス

DPに載せる関数もクラス化します。
Polynomial2つと、有理数の定数ひとつを持っておけばよいです。
ここにも足し算の機能はつけておきます。

class func {
public:
    Polynomial first, second;
    modInt third;
    func() {}
    func(const Polynomial& first, const Polynomial& second, const modInt& third) :first(first), second(second), third(third) {}
    void set(const Polynomial& first, const Polynomial& second, const modInt& third) {
        this->first = first;
        this->second = second;
        this->third = third;
    }
    func& operator+=(const func& f) {
        first += f.first;
        second += f.second;
        third += f.third;
        return *this;
    }
};

定積分メソッド

DPの遷移で使う定積分をメソッド化しておきます。
積分区間は高々2種類なので、bool値で管理します。
前に手計算した積分（下の式）をコードに写せばよいです。