Scary Sum - メモ

他人の解説記事焼き直しシリーズ

min25さんの記事を読んでできることは知っていたし、確かオイラーでも1度見かけてad-hocに再帰して解いたような気がするんだけど、ちゃんと履修していなかったところ、sounansyaさんによるわかりやすい解説が生えていたので自分用に数式を打ち直したもの。
min25さんの記事より計算量は悪いものの、あちらは実装例なし＋2次元FFTとか言っているのでまずはこちらで入門……

他人が書いた数式って気持ちが微妙に違うから、こういう複雑な数式を追うときは自分で書き直したくなるんだよね。

Scary Sum - memo

floor sum アルゴリズムとその一般化 #競技プログラミング - Qiita

問題

$p,q \in \mathbb{Z}_{\geq 0},\ N, M \in\mathbb{Z}_{>0}, \ A, B \in\mathbb{Z}$ が与えられるので
$\displaystyle f_{p,q}(N,M,A,B):=\sum_{i=0}^{N-1}i^p\left\lfloor\frac{Ai+B}{M}\right\rfloor^q$ を高速に求めたい。

$0\leq A,B < M$ への帰着

$A_1, B_1 \in \mathbb{Z}, \ A_2,B_2 \in \{0,1,\ldots,M-1\}$ を、 $A=A_1M+A_2,\ B=B_1M+B_2$ となるよう取る。元の式に代入して
$\displaystyle{f_{p,q}(N,M,A,B)\\ =\sum_{i=0}^{N-1}i^p\left(\left\lfloor\frac{A_2i+B_2}{M}\right\rfloor+A_1i+B_1\right)^q\\ =\sum_{0\leq r+s\leq q}\binom{q}{r,s}A_1^sB_1^{q-r-s}\sum_{i=0}^{N-1}i^{p+s}\left\lfloor\frac{A_2i+B_2}{M}\right\rfloor^r}$
を得る。よって以下では $0\leq A,B < M$ であるとする。

総和公式

1変数多項式 $S_p$ を $S_p(N)=\sum_{n=0}^{N-1}n^p$ を満たすものとし、その係数を有理数 $s_{p,i}$ により $S_p(X)=\sum_{i=0}^{p+1}s_{p,i}X^i$ と表す。
q=0 のとき、答えは $S_p(N)$ である。以下では q > 0 であるとする。

主客転倒

floorの部分の値で主客転倒する。取りうる値は 0 以上 $K:=\displaystyle \left\lfloor\frac{A(N-1)+B}{M}\right\rfloor$ 以下である。
K=0のとき答えは明らかに0である。以下では K > 0 であるとする。
$d_0,\ldots,d_{K+1} \in \mathbb{Z}$ を、
$d_0=0$
$d_{K+1}=N$
$1\leq k \leq K$ のとき $d_k=\displaystyle \left\lceil\frac{Mk-B}{A}\right\rceil=\left\lfloor\frac{Mk-B+A-1}{A}\right\rfloor$
と定めると、
「 $0\leq i < N$ かつ $\displaystyle \left\lfloor\frac{Ai+B}{M}\right\rfloor=k$ 」 $\Longleftrightarrow i \in [d_k,d_{k+1})$ となることがわかる。
これはfloorとceilの処理に気をつけながら式変形するだけなので詳細は割愛する。

本計算

$\displaystyle{ f_{p,q}(N,M,A,B) \\ =\sum_{i=0}^{N-1}i^p\left\lfloor\frac{Ai+B}{M}\right\rfloor^q \\ =\sum_{k=1}^{K}k^q\sum_{i\in[d_k,d_{k+1})}i^p \\ =\sum_{k=1}^{K}k^q(S_p(d_{k+1})-S_p(d_k)) \\ }$
アーベルの級数変形法（これも主客転倒か？）により
$\displaystyle{ =K^qS(d_{k+1})-\sum_{k=0}^{K-1}((k+1)^q-k^q)S_p(d_{k+1})\\ =K^qS(N)-\sum_{k=0}^{K-1}\left(\sum_{i=0}^{q-1}\binom{q}{i}k^i\right)\left(\sum_{j=0}^{p+1}s_{p,j}d_{k+1}^j\right)\\ =K^qS_p(N)-\sum_{i,j}\binom{q}{i}s_{p,j}\sum_{k=0}^{K-1}k^id_{k+1}^j\\ =K^qS_p(N)-\sum_{i,j}\binom{q}{i}s_{p,j}\sum_{k=0}^{K-1}k^i\left\lfloor\frac{Mk+M-B+A-1}{A}\right\rfloor^j\\ =K^qS_p(N)-\sum_{i,j}\binom{q}{i}s_{p,j}f_{i,j}(K,A,M,M-B+A-1) }$
となり再帰的に解くことで求まる。

再帰のベースケースは

N=0 → 空和なので0
q=0 → $S_p(N)$
（q≠0 かつ K=0 → 0の和なので0）　←これはなくても良い。1回再帰するとN=0になるので

計算量

再帰呼び出しするときの引数はp,qによらずN,M,A,Bのみによるので、 $f_{p,q}(N,M,A,B)$ を計算する過程全体で呼び出される引数の種類数は $O(\log M)$ であり、p,qも合わせると評価する関数の値の個数は $O((p+q)^2\log M)$ 個であることがわかる。評価自体は総和が $O(pq)$ でできるので、全体の計算量は $O((p+q)^4\log M)$ になることが確かめられた。

実装

Pythonでの実装
 有理数クラスを使うなどいい加減な実装をしても (p,q)=(5,5), A,B,N,M≦1000 が 100ケースで4秒程度だった。