Day 070-Q4 — 多項式剰余木・多点評価（Polynomial Remainder Tree + Subproduct Tree）

2026-06-23 赤色 Master / Phase 8+ ★★★★★★★★★ 多点評価・NTT・Subproduct Tree・FPS

問題

次数 $N-1$ の多項式 $f(x) = \sum_{i=0}^{N-1} c_i x^i$ と $M$ 個の評価点 $x_1, \ldots, x_M$ が与えられる。$f(x_1), \ldots, f(x_M)$ をすべて $p = 998244353$ で割った余りを求めよ。

ナイーブな $O(NM)$ ではなく、Subproduct Tree（多項式剰余木）を使い $O((N+M)\log^2(N+M))$ で解け。

制約

パラメータ範囲
$N, M$$1 \le N, M \le 2 \times 10^5$
$c_i$$0 \le c_i < p$
$x_j$$0 \le x_j < p$
$p$$998244353$（NTT 対応素数）

入出力例

入力例 1

4 3
1 2 3 4
0 1 2

出力例 1

1
10
49

$f(x) = 4x^3+3x^2+2x+1$。$f(0)=1, f(1)=10, f(2)=49$。

入力例 2

3 4
5 0 3
1 2 3 4

出力例 2

8
17
32
53

$f(x) = 3x^2+5$。$f(1)=8, f(2)=17, f(3)=32, f(4)=53$。

概念図: Subproduct Tree と剰余伝播

Subproduct Tree: 評価点 x₁,x₂,x₃,x₄ (M=4) の例 葉: (x-x_j) | 内部ノード: 子2つの積 | 根: (x-x₁)(x-x₂)(x-x₃)(x-x₄) 根から f を剰余で分割し葉まで伝播 → 葉での剰余 = f(x_j) (x-x₁)(x-x₂)(x-x₃)(x-x₄) 根: 次数4の多項式 (x-x₁)(x-x₂) 次数2 (x-x₃)(x-x₄) 次数2 (x - x₁) 葉: f mod (x-x₁) (x - x₂) 葉: f mod (x-x₂) (x - x₃) 葉: f mod (x-x₃) (x - x₄) 葉: f mod (x-x₄) = f(x₁) mod p = f(x₂) mod p = f(x₃) mod p = f(x₄) mod p 重要な恒等式: f mod (x - x_j) = f(x_j)（剰余定理） 葉での剰余が 0 次多項式（定数）になり、それが評価値

ヒント（段階的開示）

ヒント1: 方向性
多点評価の高速化は Subproduct Tree（積多項式木）を使う。$M$ 個の評価点から積 $\prod_{j=1}^{M}(x - x_j)$ を分割統治で構築し、$f(x)$ を各リーフで剰余を取ることで評価値を得る。
ヒント2: アプローチ
  1. Build: 葉 = $(x - x_j)$、内部ノード = 子の積でツリー構築 $O(M \log^2 M)$
  2. Evaluate: 根から $f$ を分割。各ノードで親の多項式を子で割った剰余を子に渡す $O(M \log^2 M)$
  3. 葉での剰余 $= f \bmod (x - x_j) = f(x_j)$（剰余定理）
ヒント3: コード骨格 
# Subproduct Tree の構築
size = 1
while size < M: size <<= 1
tree = [None] * (2 * size)
for i in range(M):
    tree[size + i] = [(-xs[i]) % MOD, 1]  # x - x_j
for i in range(M, size):
    tree[size + i] = [1]  # padding
for i in range(size - 1, 0, -1):
    tree[i] = poly_mul(tree[2*i], tree[2*i+1])

# 根から剰余を伝播
tree[1] = poly_div_mod(f, tree[1])[1]
for i in range(1, size):
    _, rem_l = poly_div_mod(tree[i], tree[2*i])
    tree[2*i] = rem_l
    _, rem_r = poly_div_mod(tree[i], tree[2*i+1])
    tree[2*i+1] = rem_r

模範解答 (Python)

import sys
MOD = 998244353

def modinv(a, m=MOD): return pow(a, m-2, m)

def ntt(a, invert):
    n = len(a)
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j & bit: j ^= bit; bit >>= 1
        j ^= bit
        if i < j: a[i], a[j] = a[j], a[i]
    length = 2
    while length <= n:
        w = 3 if not invert else modinv(3)
        w = pow(w, (MOD-1)//length, MOD)
        for i in range(0, n, length):
            wn = 1
            for k in range(length//2):
                u = a[i+k]; v = a[i+k+length//2]*wn%MOD
                a[i+k] = (u+v)%MOD
                a[i+k+length//2] = (u-v)%MOD
                wn = wn*w%MOD
        length <<= 1
    if invert:
        inv_n = modinv(n)
        for i in range(n): a[i] = a[i]*inv_n%MOD

def poly_mul(a, b):
    rl = len(a)+len(b)-1; n = 1
    while n < rl: n <<= 1
    fa = a+[0]*(n-len(a)); fb = b+[0]*(n-len(b))
    ntt(fa, False); ntt(fb, False)
    for i in range(n): fa[i] = fa[i]*fb[i]%MOD
    ntt(fa, True); return fa[:rl]

def poly_inv(a, deg):
    r = [modinv(a[0])]; length = 1
    while length < deg:
        length <<= 1
        ar = poly_mul(a[:length], r)[:length]
        rr = poly_mul(r, r)[:length]
        arr = poly_mul(a[:length], rr)[:length]
        new_r = [(2*r[i] if i < len(r) else 0)%MOD for i in range(length)]
        for i in range(len(arr)): new_r[i] = (new_r[i]-arr[i])%MOD
        r = new_r[:length]
    return r[:deg]

def poly_div_mod(a, b):
    if len(a) < len(b): return [0], a[:]
    n = len(a)-len(b)+1
    a_rev = a[::-1][:n]; b_rev = b[::-1][:n]
    q_rev = poly_mul(a_rev, poly_inv(b_rev, n))[:n]
    q = q_rev[::-1]
    qb = poly_mul(q, b)
    r = [(a[i]-(qb[i] if i= len(tree[1]) else f[:]
    for i in range(1, size):
        if tree[i] is None: continue
        _, r = poly_div_mod(tree[i], tree[2*i]); tree[2*i] = r
        _, r = poly_div_mod(tree[i], tree[2*i+1]); tree[2*i+1] = r
    return [(tree[size+i][0]%MOD if tree[size+i] else 0) for i in range(M)]

def solve():
    data = sys.stdin.read().split()
    idx = 0
    N = int(data[idx]); idx += 1
    M = int(data[idx]); idx += 1
    f = [int(data[idx+i])%MOD for i in range(N)]; idx += N
    xs = [int(data[idx+i])%MOD for i in range(M)]; idx += M
    print(*multipoint_eval(f, xs), sep='\n')

solve()

Step-by-Step 解説

Step 1: 剰余定理の応用

任意の多項式 $f(x)$ に対し $f(x) = (x - a) \cdot q(x) + f(a)$ が成立（剰余定理）。よって $f \bmod (x-a) = f(a)$。

Step 2: Subproduct Tree の構築

葉から根に向かって積を計算。深さ $O(\log M)$、各レベルで多項式サイズの和が $O(M)$。NTT を使った多項式乗算 $O(M \log M)$ × $O(\log M)$ レベル = $O(M \log^2 M)$。

Step 3: 剰余の根から葉への伝播

$f$ を根の多項式で割った余り $r$ を計算し、$r$ を左右の子で割った余りを子に渡す。葉まで伝播すると $f(x_j)$ が得られる。

Step 4: 多項式割り算の実装

次数 $n$ の多項式の逆元を Newton 法で計算（$O(n \log n)$）し、割り算に使う。全体 $O(M \log^2 M)$。

Step 5: 計算量比較

手法計算量
ナイーブ$O(NM)$
Subproduct Tree$O((N+M)\log^2(N+M))$

よくあるミス

ミス原因正しい書き方
葉の多項式を $(x_j - x)$ にする符号ミス[(-xs[i]) % MOD, 1] = $(x - x_j)$
葉の剰余が空リストになる次数0の多項式処理r[0] if r else 0
NTT の mod 対応を確認しない$998244353$ 以外では素朴 NTT 不可$p = 998244353$ を使用
多項式逆元の Newton 法バグ実装が複雑テスト用に手計算で検証

次のステップ

  • 発展問題: Lagrange 補間の高速化（多点評価 + 剰余木で $O(N \log^2 N)$）
  • FPS exp/log との組み合わせ
  • Bostan-Mori アルゴリズムとの比較
  • 多項式の多点評価 + 内挿の一体実装

自己評価