競プロ典型90問5日目

# 競プロ典型90問5日目

https://raw.githubusercontent.com/E869120/kyopro_educational_90/main/problem/005.jpg

https://atcoder.jp/contests/typical90/tasks/typical90_e

# 定式化

$$
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)$ の定数項が解である。ここで漸化式
```math
\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 アルゴリズムより定数倍高速なのは[二種類の繰り返し二乗法](https://qiita.com/ryuhe1/items/da5acbcce4ac1911f47a#%E4%BA%8C%E7%A8%AE%E9%A1%9E%E3%81%AE%E7%B9%B0%E3%82%8A%E8%BF%94%E3%81%97%E4%BA%8C%E4%B9%97%E6%B3%95)と同様である。

## アルゴリズム 1
### LSB-first アルゴリズム

$N$が正の偶数のとき
```math
\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$ である。

```python
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$が正の偶数のとき
```math
\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$ の大きさによる。

```python
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$が正の偶数のとき
```math
\mathrm{prod}(N, D, Q)=
\mathrm{prod}\left(N/2, D^2\bmod B, Q(x) Q(x^{D}) \bmod x^B-1\right)
```

```python
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$が正の偶数のとき
```math
\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)
```

```python
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周期分だけ計算してべき乗する方法も考えられるが、最悪ケースでは計算量の改善は無さそうである。