無限和を計算するyosupo judgeの問題

スナネコです。yosupo judgeにある $\sum_{i=0}^{\infty}r^i i^d$ を解きました。
問題

d=0の時はただの等比数列の和です。
d=1の時を計算する問題もときどき見ますね。例えば「1の目が出るまでサイコロを振るときの回数の期待値は？」とかです。
…………冗談です。まさかこの問題でそんな計算しませんよね？

さて、dが大きなときでもやることはd=0,1の時と同じで、rを掛けて差を取っていきます。
例えばd=4のときは
$\begin{eqnarray} S &=& 1*r &+& 16*r^2 &+& 81*r^3 &+& 256*r^4 &+& 625*r^5 &+& \cdots \\ rS &=& \ &\ & 1*r^2 &+& 16*r^3 &+& 81*r^4 &+& 256*r^5 &+& \cdots \\ (1-r)S &= & 1*r &+& 15*r^2 &+& 65*r^3 &+& 175*r^4 &+& 369*r^5 &+& \cdots \\ r(1-r)S &=&\ &\ & 1*r^2 &+& 15*r^3 &+& 65*r^4 &+& 175*r^5 &+& \cdots \\ (1-r)^2 S &= & 1*r &+& 14*r^2 &+& 50*r^3 &+& 110*r^4 &+& 194*r^5 &+& \cdots \\ r(1-r)^2 S &=&\ &\ & 1*r^2 &+& 14*r^3 &+& 50*r^4 &+& 110*r^5 &+& \cdots \\ (1-r)^3 S &= & 1*r &+& 13*r^2 &+& 36*r^3 &+& 60*r^4 &+& 84*r^5 &+& \cdots \\ r(1-r)^3 S &=&\ &\ & 1*r^2 &+& 13*r^3 &+& 36*r^4 &+& 60*r^5 &+& \cdots \\ (1-r)^4 S &= & 1*r &+& 12*r^2 &+& 23*r^3 &+& 24*r^4 &+& 24*r^5 &+& \cdots \\ r(1-r)^4 S &=&\ &\ & 1*r^2 &+& 12*r^3 &+& 23*r^4 &+& 24*r^5 &+& \cdots \\ (1-r)^5 S &= & 1*r &+& 11*r^2 &+& 11*r^3 &+& 1*r^4 \end{eqnarray}$
という感じですね。
なんだか見たことのない数が残りました。こういうときはOEISです。
1,11,11,1 - OEIS
Eulerian Numberというらしいです。d=5の場合も計算すると一致したので多分あってるでしょう。
wikipediaを見てみます。
Eulerian number - Wikipedia
wikipediaに合わせて、Eulerian NumberをA(n,m)と表すことにしましょう。
スナネコは英語は読めませんが、数式を見れば大体の雰囲気はわかりますね。
A(n,m)を表す式は書いてありますが、nを固定したときに全部のmを高速に求める方法は書いてなさそうです。

このEulerian Numberを使って計算できればいいんですけど、うまくいきますかね。ちゃんと式に書いてみましょう。
$\displaystyle{ \sum_{m=0}^{d-1}A(d,m)r^{m+1}\\ =\sum_{m=0}^{d-1}\sum_{k=0}^{m} (-1)^k\binom{d+1}{k}(m+1-k)^d r^{m+1}\\ }$
愚直に計算するとO(d^2)です。困りました。
正解者の提出をちらっと見てみると、こういう式で計算してました。
$\displaystyle{ \sum_{i=0}^d (-r)^{d-i} \binom{d+1}{i+1} \sum_{j=0}^{i} r^j j^d }$
スナネコの考えてる式と似たような感じなので、今の方針であってそうです。
二項係数の部分と、(…)^dの部分で変数を分離して、独立に累積和っぽく計算する工夫は面白いですね。これが使えないか試してみます。m+1-kをxと変数変換してみます。
$\displaystyle{ \sum_{m=0}^{d-1}\sum_{k=0}^{m} (-1)^k\binom{d+1}{k}(m+1-k)^d r^{m+1}\\ =\sum_{k=0}^{d-1}\sum_{m=k}^{d-1} (-1)^k\binom{d+1}{k}(m+1-k)^d r^{m+1}\\ =\sum_{k=0}^{d-1}\sum_{x=1}^{d-k} (-1)^k\binom{d+1}{k}x^d r^{x+k}\\ =\sum_{k=0}^{d-1} (-1)^k\binom{d+1}{k} \sum_{x=1}^{d-k}x^d r^{x+k}\\ =\sum_{k=0}^{d-1} (-r)^k\binom{d+1}{k} \sum_{x=1}^{d-k}r^x x^d\\ }$
うまくできましたね。これで、xのシグマについてはkが大きい方から累積和のように計算していくことができるので、高速に計算できるようになりました。
この計算式の通りにやっても全体で $O(d\log mod)$ ですが、逆元テーブルの線形構築や細かい計算を配列にメモしておくことでいくらか定数倍高速化が効きます。

これでAC……ではありません。d=0の時はそもそも最初の $\sum_{m=0}^{d-1}A(d,m)r^{m+1}$ の時点で空和になって0になってしまいます。仕方がないのでこれだけ場合分けです。

ACは取れましたが、途中で紹介した他の正解者の計算式では場合分けをしてないみたいです。こっちも導出してみましょう。
最後の数式でd-kをyに変数変換します。
$\displaystyle{ \sum_{k=0}^{d-1} (-r)^k\binom{d+1}{k} \sum_{x=1}^{d-k}r^x x^d\\ =\sum_{y=1}^{d} (-r)^{d-y}\binom{d+1}{d-y} \sum_{x=1}^{y}r^x x^d\\ =\sum_{y=1}^{d} (-r)^{d-y}\binom{d+1}{y+1} \sum_{x=1}^{y}r^x x^d\\ }$
和を取る範囲が違いますね。ここが0から始まっているとd=0のときにも場合分けをしなくて済むようです。
wikipediaをよく見てみると、A(0,0)=1、n>0のときA(n,n)=0、と書いてありました。ということはそもそも一番最初の式の時点でm=dまで和を取ってよくて、そうすると整合性がありそうですね。
$\displaystyle{ \sum_{m=0}^{d}A(d,m)r^{m+1}\\ =\sum_{m=0}^{d}\sum_{k=0}^{m} (-1)^k\binom{d+1}{k}(m+1-k)^d r^{m+1}\\ =\sum_{k=0}^{d} (-r)^k\binom{d+1}{k} \sum_{x=0}^{d-k}r^x x^d\\ =\sum_{y=0}^{d} (-r)^{d-y}\binom{d+1}{y+1} \sum_{x=0}^{y}r^x x^d\\ }$
確かに同じ式が求まりました。