競プロ典型90問5日目

定式化

$$
Q(x) := \sum_{i=1}^K x^{c_i}
$$

とおくと、

$$
G(x) := Q(x) Q(x^{10}) Q(x^{100}) \dotsm Q(x^{10^{N-1}}) \bmod x^B-1
$$

の定数項が解である。

アルゴリズム

概要

$$
\mathrm{prod}(N, D, Q) := Q(x) Q(x^{D}) Q(x^{D^2}) \dotsm Q(x^{D^{N-1}}) \bmod x^B-1
$$

と定義すると、$\mathrm{prod}(N, 10, Q)$ の定数項が解である。ここで漸化式

\mathrm{prod}(N, D, Q)=
\begin{cases}
1,& \text{if }N = 0\\
Q(x)\cdot\mathrm{prod}(N-1, D, Q(x^D)\bmod x^B-1)\bmod x^B-1,& \text{if $N$ is odd}\\
\cdots,& \text{otherwise}
\end{cases}

が成り立つので、後は$N$が正の偶数の場合の漸化式を考えればよい。

以下では本質的に二つの異なるアルゴリズムを紹介する。この二つのアルゴリズムの算術計算量は等しい。また、そのそれぞれに LSB-first と MSB-first のアルゴリズムがあるので合計で四つのアルゴリズムを紹介する。MSB-first のアルゴリズムは $Q$ が疎であるときに、LSB-first のアルゴリズムよりも高速である。MSB-first アルゴリズムが特殊ケースで LSB-first アルゴリズムより定数倍高速なのは二種類の繰り返し二乗法と同様である。

アルゴリズム 1

LSB-first アルゴリズム

$N$が正の偶数のとき

\mathrm{prod}(N, D, Q)=
\mathrm{prod}\left(N/2, D, Q(x) Q(x^{D^{N/2}})\bmod x^B-1\right)

二つの高々 $B-1$ 次の多項式を掛けて $\bmod x^B-1$ をとるための計算量は定数 $c$ について $\mathsf{M}(B-1)+cB$ で抑えられる。
よって上記の漸化式に基づくアルゴリズムの算術計算量 $T(N)$ は

$$
T(N) \le 2(\mathsf{M}(B-1) + cB) + T(\lfloor N/2\rfloor)
$$

を満たすので、 $T(N)\le 2(\mathsf{M}(B-1)+cB)\lfloor\log N\rfloor$ である。

N, B, K = map(int, input().split())
C = list(map(int, input().split()))

mod = 1000000007
F = GF(mod)
R.<x> = F[]

# P(x^b) mod x^B-1
def reduce(P, b = 1):
  l = P.list()
  r = [F(0)] * B
  for i in range(len(l)):
    r[i*b%B] += l[i]
  return R(r)

def prod(N, Q):
  A = R(1);
  while N:
    if N % 2 == 1:
      A = reduce(A * Q)
      Q = reduce(Q, 10)
    b = pow(10, N//2, B)
    Q = reduce(Q * reduce(Q, b))
    N //= 2
  return A[0]

P = R(0)
for c in C:
  P += x^(c % B)

print(prod(N, P))

MSB-first アルゴリズム

$N$が正の偶数のとき

\mathrm{prod}(N, D, Q)= A(x)A(x^{D^{N/2}}) \bmod x^B-1\quad \text{for }
A(x):=\mathrm{prod}\left(N/2, D, Q(x)\right)

算術計算量は LSB-first アルゴリズムと同様に $T(N)\le 2(\mathsf{M}(B-1)+cB)\lfloor \log N\rfloor$ となるが、$Q(x)$ が疎であることを利用すると、$N$ が奇数のときに発生する乗算の算術計算量を $KB$ で抑えられる。よって $T(N)\le (\mathsf{M}(B-1)+cB+KB)\lfloor \log N\rfloor$ となる。ここで $\mathsf{M}(B-1)+cB$ と $KB$ のどちらが小さいかは $K$ と $B$ の大きさによる。

N, B, K = map(int, input().split())
C = list(map(int, input().split()))

mod = 1000000007
F = GF(mod)
R.<x> = F[]

# P(x^b) mod x^B-1
def reduce(P, b = 1):
  l = P.list()
  r = [F(0)] * B
  for i in range(len(l)):
    r[(i*b)%B] += l[i]
  return R(r)

def prod(N, Q):
  if N == 1:
    return Q
  if N % 2 == 1:
    return reduce(Q * prod(N-1, reduce(Q, 10)))
  A = prod(N//2, Q)
  b = pow(10, N//2, B)
  A = reduce(A * reduce(A, b))
  return A

P = R(0)
for c in C:
  P += x^(c % B)

print(prod(N, P)[0])

アルゴリズム2

LSB-first アルゴリズム

$N$が正の偶数のとき

\mathrm{prod}(N, D, Q)=
\mathrm{prod}\left(N/2, D^2\bmod B, Q(x) Q(x^{D}) \bmod x^B-1\right)

N, B, K = map(int, input().split())
C = list(map(int, input().split()))

mod = 1000000007
F = GF(mod)
R.<x> = F[]

# P(x^b) mod x^B-1
def reduce(P, b = 1):
  l = P.list()
  r = [F(0)] * B
  for i in range(len(l)):
    r[i*b%B] += l[i]
  return R(r)

def prod(N, D, Q):
  A = R(1);
  while N:
    if N % 2 == 1:
      A = reduce(A * Q)
      Q = reduce(Q, D)
    Q = reduce(Q * reduce(Q, D))
    D = D * D % B
    N //= 2
  return A[0]

P = R(0)
for c in C:
  P += x^(c % B)

print(prod(N, 10, P))

MSB-first アルゴリズム

$N$が正の偶数のとき

\mathrm{prod}(N, D, Q)= A(x)A(x^D) \bmod x^B-1 \quad\text{ for }
A(x):=\mathrm{prod}\left(N/2, D^2\bmod B, Q(x)\right)

N, B, K = map(int, input().split())
C = list(map(int, input().split()))

mod = 1000000007
F = GF(mod)
R.<x> = F[]

# P(x^b) mod x^B-1
def reduce(P, b = 1):
  l = P.list()
  r = [F(0)] * B
  for i in range(len(l)):
    r[i*b%B] += l[i]
  return R(r)

def prod(N, D, Q):
  if N == 1:
    return Q
  if N % 2 == 1:
    return reduce(Q * prod(N-1, D, reduce(Q, D)))
  A = prod(N//2, D * D % B, Q)
  return reduce(A * reduce(A, D))

P = R(0)
for c in C:
  P += x^(c % B)

print(prod(N, 10, P)[0])

その他の手法

等比数列 $1,10,100,1000,\dotsc$ は $\bmod B$ で周期的になるので(途中で0になるケースには注意)、1周期分だけ計算してべき乗する方法も考えられるが、最悪ケースでは計算量の改善は無さそうである。