えびちゃんの日記

えびちゃん（競プロ）の日記です。

ある種の無限級数を DP で求める

というやり取りがありました。

$$ \sum_{k=1}^{\infty} \frac{k^{n+k-1}}{k!\, 2^k\, e^{k/2}} $$ という無限級数を求めたいです。ここでは $O(n^2)$ 時間の方法を述べますが、頭のいい人がちゃんと考察をすればもっと高速になるのかもしれません。

動機

たぶんここは読み飛ばしてもいいです。上記の無限級数がどうして求めたくなったかの話を書きます。

数え上げの講義資料（？） の 12 章の演習問題から抜粋です。 $N_n = \{0, 1, \dots, n-1\}$ の部分集合全体からなる集合 $2^{N_n}$ を考えます。

この部分集合 $\Delta\subseteq 2^{N_n}$ であって、$A, B \in\Delta \implies (A\cap B = \emptyset) \vee (A \subseteq B) \vee (A\supseteq B)$ を満たすもの全体の集合を $\mathcal{N}_n$ とおきます。$\Delta$ は、うまくネストされた部分集合たちという感じです。

たとえば、$\mathcal{N}_0 = \{ \{\}, \{\{\}\} \} = 2^{2^{N_0}}$、$\mathcal{N}_1 = \{ \{\}, \{\{\}\}, \{\{0\}\}, \{\{\}, \{0\}\} \} = 2^{2^{N_1}}$、$\mathcal{N}_2 = 2^{2^{N_2}}$ です。

$\mathcal{N}_3$ の要素としては、$\{\{0\}, \{0, 2\}\}$ や $\{\{1\}\}$ や $\{\{\}, \{0, 1\}, \{2\}, \{0, 1, 2\}\}$ や $\{\{\}, \{0\}, \{1\}, \{2\}, \{0, 2\}, \{0, 1, 2\}\}$ などがあります。一方 $2^{2^{N_3}} \smallsetminus \mathcal{N}_3$ の要素としては、$\{\{\}, \{0, 1\}, \{1, 2\}, \{0, 1, 2\}\}$ や $\{\{\}, \{0\}, \{1\}, \{0, 1\}, \{2\}, \{0, 2\}, \{1, 2\}, \{0, 1, 2\}\} = 2^{N_3}$ などがあります。

これに対してがちゃがちゃすると、 $$ \gdef\card#1{\#{#1}} \card{\mathcal{N}_n} = 4\sum_{k=1}^{\infty} \frac{k^{n+k-1}}{k!\, 2^k\, e^{k/2}} $$ が言えるらしいです。というわけでこれを求めたいです。

考察

$W$ 関数との関連

Lambert の $W$ 関数 というのがあります。これは、$w(x) = xe^x$ の逆関数で、$x = W(x) e^{W(x)}$ です。 たとえば、$w(-1/2) = -\frac{1}{2}\cdot e^{-1/2} = -\frac{1}{2\sqrt{e}}$ なので、$W{\left(-\frac{1}{2\sqrt{e}}\right)} = -1/2$ です。

これに対して、 $$ W(x) = \sum_{k=1}^{\infty} \frac{(-k)^{k-1}}{k!} x^k $$ という表現が知られています。 $$ \begin{aligned} -W(-x) &= -\sum_{k=1}^{\infty} \frac{(-k)^{k-1}}{k!} (-x)^k \\ &= \sum_{k=1}^{\infty} \frac{k^{k-1}}{k!} x^k \\ \end{aligned} $$ が成り立つので、$x = a \coloneqq \frac{1}{2\sqrt{e}}$ として $$ -W(-a) = \sum_{k=1}^{\infty} \frac{k^{k-1}}{k!\, 2^k\, e^{k/2}} = \frac{\card{\mathcal{N}_0}}{4} $$ となります。$W(-a) = -1/2$ より、$\card{\mathcal{N}_0} = 2$ とわかります。

さて、 $$ f_0(x) = -W(-x) = \sum_{k=1}^{\infty} \frac{k^{k-1}}{k!} x^k $$ とおき、$f_{i+1} = x\cdot f_i'(x)$ ($i = 0, 1, \dots$) と定めると、 $$ f_n(x) = \sum_{k=1}^{\infty} \frac{k^{n+k-1}}{k!} x^k $$ となります。

$W$ 関数の微分

$\gdef\dd{\mathrm{d}}$

逆関数の微分などを考えます。 $y = W(x)$ とすると $x = w(y) = ye^y$ で、$\dd x/\dd y = (y+1)e^y$ より、 $$ \frac{\dd y}{\dd x} = \frac{1}{\frac{\dd x}{\dd y}} = \frac{1}{(y+1)e^y} $$ です。$e^y = x/y$ に注意すると、 $$ W'(x) = \frac{W(x)}{x\cdot (1+W(x))} $$ となります。

$f_0(x) = -W(-x)$ であったので、 $$ \begin{aligned} f_1(x) &= x\cdot(f_0(x))' \\ &= x\cdot(-W(-x))' \\ &= x\cdot W'(-x) \\ &= -\frac{W(-x)}{1+W(-x)} \end{aligned} $$ とわかります。

ここで、$V(x) \coloneqq 1+W(-x)$ とおくと、$f_1(x) = -W(-x)/V(x)$ と書けます。

DP

手でがちゃがちゃやったりしていると、一般に $W(-x)^j/V(x)^{2i-1}$ の微分（に $x$ を掛けたもの）を求めたいな〜という気持ちになります。 計算過程は煩雑なので、結果だけ載せます。商の微分公式や合成関数の微分公式などを使いながらがんばるとよいです。

$$ \begin{aligned} x\cdot \left(\frac{W(-x)^j}{V(x)^{2i-1}}\right)' &= \frac{W(-x)^j}{V(x)^{2i+1}}\cdot(j+(j-(2i-1))\cdot W(-x)) \\ &= j\cdot\frac{W(-x)^j}{V(x)^{2i+1}} +(j-(2i-1))\cdot \frac{W(-x)^{j+1}}{V(x)^{2i+1}}. \\ \end{aligned} $$

$\gdef\DP{\mathrm{dp}}$ ここで、$f_i$ における $W(-x)^j/V(x)^{2i-1}$ の係数を $\DP[i][j]$ で定めます。 $f_1(x) = -W(-x)/V(x) = (-1)\cdot W(-x)^1/V(x)^{2\cdot 1-1}$ だったので、$\DP[1][1] = -1$ です。

上記の微分より、 $$ \gdef\xgets#1{\xleftarrow{#1}} \begin{aligned} \DP[i+1][j] &\xgets{+} j\cdot\DP[i][j]; \\ \DP[i+1][j+1] &\xgets{+} (j-(2i-1))\cdot\DP[i][j] \end{aligned} $$ の更新を $(i, j)$ の昇順に行えばよいです。

さて、$W(-a) = -1/2$ および $V(a) = 1 + W(-a) = 1/2$ より、 $$ \begin{aligned} \left.\frac{W(-x)^j}{V(x)^{2i-1}}\right|_{x=a} &= \frac{(-1/2)^j}{(1/2)^{2i-1}} \\ % &= \frac{2^{2i-1}}{(-2)^j} \\ % &= -\frac{(-2)^{2i-1}}{(-2)^j} \\ &= -(-2)^{2i-j-1} \end{aligned} $$ となるので、$n \ge 1$ に対して $$ \begin{aligned} \card{\mathcal{N}_n} &= 4\sum_{j=1}^n \DP[n][j] \cdot (-(-2)^{2n-j-1}) \\ &= -\sum_{j=1}^n \DP[n][j] \cdot (-2)^{2n-j+1} \end{aligned} $$ とわかります。

あるいは、最後に $(-2)^{*}$ を掛けるのではなく、遷移しながら $-2$ を掛けていくような別の方針の DP もできるかもしれません。 そちらの方針についてはあまり考えていません。

実装

実装自体は、（遷移が与えられていれば）DP 初級編くらいの内容な気がします。すぐにめちゃ大きくなるので、998244353 を法として求めます。

use proconio::input;

type Mi = ModInt998244353;

fn main() {
    input! {
        n: usize,
    }
    let mut dp = vec![vec![Mi::new(0); n + 2]; n + 1];
    dp[1][1] = Mi::new(-1);
    for i in 0..n {
        for j in 1..=i {
            let x = dp[i][j];
            dp[i + 1][j] += Mi::new(j) * x;
            dp[i + 1][j + 1] += (Mi::new(j) - Mi::new(2 * i - 1)) * x;
        }
    }

    let powm2 = {
        let m = 2 * n;
        let mut pow = vec![Mi::new(1); m + 1];
        for i in 1..=m {
            pow[i] = pow[i - 1] * Mi::new(-2);
        }
        pow
    };

    let sum: Mi = (1..=n).map(|j| dp[n][j] * powm2[2 * n - j + 1]).sum();
    println!("{}", -sum);
}

$(\card{\mathcal{N}_n})_{n=0}^{\infty} = (2, 4, 16, 128, 1664, 30208, \dots)$ となりそうです。

所感

自分語りのコーナーです。 唐突ですが、$e^{e^x}$ の $i$ 階微分は、$e^{e^x + jx}$ ($1\le j\le i$) の線形和で表せます。 たとえば、 $$ \frac{\dd^4}{\dd x^4} e^{e^x} = e^{e^x + x} + 7 e^{e^x + 2 x} + 6 e^{e^x + 3 x} + e^{e^x + 4 x} $$ といった具合です。

高校生くらいのえびちゃんは競プロや DP という用語を知りませんでしたが、これの係数を手計算で求めようとして DP をしていた記憶があります。

$$ \begin{aligned} \frac{\dd}{\dd x} e^{e^x + jx} &= (e^{e^x})'\cdot e^{jx} + e^{e^x}\cdot (e^{jx})' \\ % &= (e^{e^x+x})\cdot e^{jx} + j\cdot e^{e^x}\cdot (e^{jx}) \\ &= (e^{e^x+x})\cdot e^{jx} + j\cdot e^{e^x + jx} \\ &= e^{e^x+(j+1)x} + j\cdot e^{e^x + jx} \\ \end{aligned} $$ となることから、$\DP[i+1][j] \xgets{+} j\cdot \DP[i][j]$ と $\DP[i+1][j+1] \xgets{+} \DP[i][j]$ で遷移できるわけですね。 係数に $j$ が含まれるのが大変そうな気もしますが、形自体は単純なので、なんかいい感じの手法で高速化できるのでしょうか。

ともあれ、今回出てきた DP と発想がほぼ同じで、昔のことを思い出してなつかしくなったという話でした。

おまけ

WolframAlpha すごい。

www.wolframalpha.com

www.wolframalpha.com

もう少し大きい $n$ に対しては、無限級数として与える方では、誤差が出たり TLE したりしました。

おわり

にゃんこ。