簡単のため $B = 2 ^ K$ とする。また $\displaystyle f _ i(x) \coloneqq \prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i x + j) = (2 ^ i x + 1) (2 ^ i x + 2) \cdots (2 ^ i x + 2 ^ i - 1)$ と定める。
$f _ K(0), f _ K(1), \ldots, f _ K(\lfloor P/2 ^ K\rfloor)$ を計算することが目標である。
さて、$f _ i$ と $f _ {i + 1}$ の間には次の関係が成り立つ。
$$\begin{aligned}
f _ {i + 1} (x)
&= \prod _ {j = 1} ^ {2 ^ {i + 1} - 1} (2 ^ {i + 1} x + j) \\
&= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + j)\right)\cdot (2 ^ {i + 1} x + 2 ^ i) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ {i + 1} x + 2 ^ i + j)\right) \\
&= \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x) + j)\right)\cdot 2 ^ i (2x + 1) \cdot \left(\prod _ {j = 1} ^ {2 ^ i - 1} (2 ^ i (2x + 1) + j)\right) \\
&= f _ i(2x) \cdot f _ i(2x + 1) \cdot 2 ^ i (2x + 1).
\end{aligned}$$
$M _ t = \dfrac{M _ {t + 1} - (M _ {t + 1} \bmod A _ t)}{A _ t}$ より $R _ t \coloneqq M _ {t + 1} \bmod A _ t$ と定めると $M _ {t + 1} - x = A _ t M _ t + R _ t - x$ である。
$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ {t - 1} d _ {i, j} \dfrac{P _ {j + 1}}{P _ t} \leq M _ t - \left\lceil\dfrac{\left(\sum _ {i = 1} ^ {N + 1} d _ {i, t}\right) + x - R _ t}{A _ t}\right\rceil.$$
$\displaystyle s\coloneqq \sum _ {i = t + 1} ^ {N + 1} d _ {i, t}$ と定めるとさらに次のように表せる。
$$\sum _ {i = 1} ^ {N + 1} \sum _ {j = 0} ^ {t - 1} d _ {i, j} \dfrac{P _ {j + 1}}{P _ t} \leq M _ t - \left\lceil\dfrac{s + x - R _ t}{A _ t}\right\rceil.$$
従って、$\displaystyle \sum _ {i = t + 1} ^ {N + 1} d _ {i, t} = s,\ 0\leq d _ {i, t}\lt A _ t$ なる整数列 $(d _ {t + 1, t}, d _ {t + 2, t},\ldots,d _ {N + 1, t})$ の個数を $f(s)$ とすると、次が成り立つ。
$$\mathsf{dp}(t, x) = \sum _ {s = 0} ^ {(N-t+1)(A _ t - 1)} f(s) \cdot \mathsf{dp}(t - 1, \lceil (s + x - R _ t) / A _ t \rceil).$$
$0\leq x\leq N,\ 0\leq s\leq (N - t + 1)(A _ t - 1),\ 0\leq R _ t\leq A _ t - 1$ の下で、確かに再び $0\leq \lceil (s + x - R _ t) / A _ t \rceil\leq N$ が成り立っていることに注意する。
さて、$0\leq y\leq N$ なる整数 $y$ に対して $y = \lceil (s + x - R _ t) / A _ t \rceil$ を満たす $s$ の範囲は $(y A _ t + R _ t - x) - A _ t\lt s\leq (y A _ t + R _ t - x)$ である。
従って、$f(l,r)\coloneqq \sum _ {s = l + 1} ^ r f(s)$ と定めると次が成り立つ。
$$\mathsf{dp}(t, x) = \sum _ {y = 0} ^ {N} f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x) \cdot \mathsf{dp}(t - 1, y).$$
区間 $((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x\rbrack$ の長さはちょうど $A _ t$ なので、$f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x)$ は、$0$ 以上 $A _ t$ 未満の整数からなる長さ $N - t + 2$ の数列であって、総和が $y A _ t + R _ t - x$ であるようなものの個数に一致する。
この数え上げは有名問題で、包除原理を用いることで次を得る。
$$\begin{aligned}& f((y A _ t + R _ t - x) - A _ t, y A _ t + R _ t - x) \\
&{}=\sum _ {i = 0} ^ {N - t + 2} (-1) ^ i \binom{N - t + 2}{i} \binom{(y A _ t + R _ t - x) - i A _ t + (N - t + 1)}{N - t + 1}.\end{aligned}$$
式の簡単のため $u\coloneqq N - t + 1$ と定めると、$\mathsf{dp}(t, x)$ は次のように表せる。
$$\begin{aligned} \mathsf{dp}(t, x) &{}= \sum _ {y = 0} ^ {N} \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(y A _ t + R _ t - x) - i A _ t + u}{u} \cdot \mathsf{dp}(t - 1, y)\\
&{}=\sum _ {y = 0} ^ {N} \mathsf{dp}(t - 1, y) \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(R _ t - x + (y - i) A _ t) + u}{u}
.\end{aligned}$$
さて、$i\gt y$ のとき $R _ t + (y - i) A _ t \lt 0$ より $\displaystyle \binom{(R _ t - x + (y - i) A _ t) + u}{u} = 0$ である。また、$i\gt u + 1$ のとき $\displaystyle\binom{u + 1}{i} = 0$ である。従って、上式における内側の和で $i$ が動く範囲は $0$ 以上 $y$ 以下の範囲としてよい。
つまり、多項式 $a(z) = \displaystyle \sum _ {i = 0} ^ N (-1) ^ i \binom{u + 1}{i} z ^ i$ および $b(z) = \displaystyle \sum _ {i = 0} ^ N \binom{(R _ t - x + i A _ t) + u}{u} z ^ i$ について、$\lbrack z ^ y\rbrack (ab) = \displaystyle \sum _ {i = 0} ^ {u + 1} (-1) ^ i \binom{u + 1}{i} \binom{(R _ t - x + (y - i) A _ t) + u}{u}$ が成り立つ。
計算量を考える。$a(z), b(z)$ は各 $t,x$ に対してそれぞれ $O(N)$ 時間で計算できる。$b(z)$ の計算についてもう少し説明すると、これは $\displaystyle \binom{(R _ t - x + i A _ t) + u}{u} = \binom{(R _ t - (x - 1) + i A _ t) + u}{u}\cdot \dfrac{(R _ t - x + i A _ t) + u}{(R _ t - (x - 1) + i A _ t) + 1}$ を用いて SWAG で積をスライドさせながら求めることでならし $O(N)$ 時間が達成される。
// 1 <= M < 2^31uint32_t M = ...;
// 0 <= a,b < M の仮定下で a+b mod M を計算uint32_tadd(uint32_t x, uint32_t y) {
uint32_t c = a + b;
return c >= M ? c - M : c;
}
// 0 <= a,b < M の仮定下で a-b mod M を計算uint32_tsub(uint32_t x, uint32_t y) {
uint32_t c = a + (M - b);
return c >= M ? c - M : c;
}
$a\times b \bmod M $ を剰余演算無しで計算するのは容易ではないが、Barrett Reduction はこれを実現する方法の $1$ つとして知られている。12
