結果
| 問題 | No.3619 Compositional Power with Schröder Coordinate |
| コンテスト | |
| ユーザー |
 lif4635
|
| 提出日時 | 2026-08-09 15:39:23 |
| 言語 | PyPy3 (7.3.17) |
| 結果 |
AC
|
| 実行時間 | 6,405 ms / 10,000 ms |
| + 515µs | |
| コード長 | 34,720 bytes |
| 記録 | |
| コンパイル時間 | 231 ms |
| コンパイル使用メモリ | 95,972 KB |
| 実行使用メモリ | 340,716 KB |
| 最終ジャッジ日時 | 2026-08-10 20:54:07 |
| 合計ジャッジ時間 | 34,638 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / judge1_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 6 |
ソースコード
# harurun's library（依存を含む貼り付け用コード）
# 依存: convolution/NTT998.py
"""998244353 固定の高速NTTと係数畳み込み。
`multiply(first, second)` は2つの昇べき順係数列を畳み込み、長さ
`len(first) + len(second) - 1` の係数列を返す。汎用mod判定、原始根探索、
CRTを通らず、固定したradix-4の変換表だけを使う。
"""
MOD = 998244353
PRIMITIVE_ROOT = 3
MAX_LOG = 23
_convolution_ntt998_IMAG = 911660635
_convolution_ntt998_IIMAG = 86583718
_convolution_ntt998_RATE2 = [0, 911660635, 509520358, 369330050, 332049552, 983190778, 123842337, 238493703, 975955924, 603855026, 856644456, 131300601, 842657263, 730768835, 942482514, 806263778, 151565301, 510815449, 503497456, 743006876, 741047443, 56250497, 867605899, 0]
_convolution_ntt998_IRATE2 = [0, 86583718, 372528824, 373294451, 645684063, 112220581, 692852209, 155456985, 797128860, 90816748, 860285882, 927414960, 354738543, 109331171, 293255632, 535113200, 308540755, 121186627, 608385704, 438932459, 359477183, 824071951, 103369235, 0]
_convolution_ntt998_RATE3 = [0, 372528824, 337190230, 454590761, 816400692, 578227951, 180142363, 83780245, 6597683, 70046822, 623238099, 183021267, 402682409, 631680428, 344509872, 689220186, 365017329, 774342554, 729444058, 102986190, 128751033, 395565204, 0]
_convolution_ntt998_IRATE3 = [0, 509520358, 929031873, 170256584, 839780419, 282974284, 395914482, 444904435, 72135471, 638914820, 66769500, 771127074, 985925487, 262319669, 262341272, 625870173, 768022760, 859816005, 914661783, 430819711, 272774365, 530924681, 0]
_convolution_ntt998_INVERSE_SIZE = {1: 1}
def _convolution_ntt998_check_length(size):
if size < 1 or size & size - 1:
raise ValueError('NTT length must be a positive power of two')
if size > 1 << MAX_LOG:
raise ValueError('NTT length exceeds 2^23')
def _convolution_ntt998_butterfly(values):
"""In-place forward radix-4 NTT without coefficient normalization."""
size = len(values)
_convolution_ntt998_check_length(size)
if size == 1:
values[0] %= MOD
return values
height = (size - 1).bit_length()
level = 0
mod = MOD
imag = _convolution_ntt998_IMAG
rate2 = _convolution_ntt998_RATE2
rate3 = _convolution_ntt998_RATE3
while level < height:
if height - level == 1:
width = 1 << height - level - 1
rotation = 1
for block in range(1 << level):
offset = block << height - level
for index in range(width):
left = values[offset + index]
right = values[offset + index + width] * rotation
values[offset + index] = (left + right) % mod
values[offset + index + width] = (left - right) % mod
rotation = rotation * rate2[(~block & -~block).bit_length()] % mod
level += 1
else:
width = 1 << height - level - 2
rotation = 1
for block in range(1 << level):
rotation2 = rotation * rotation % mod
rotation3 = rotation2 * rotation % mod
offset = block << height - level
for index in range(width):
value0 = values[offset + index]
value1 = values[offset + index + width] * rotation
value2 = values[offset + index + 2 * width] * rotation2
value3 = values[offset + index + 3 * width] * rotation3
difference = (value1 - value3) % mod * imag
values[offset + index] = (value0 + value2 + value1 + value3) % mod
values[offset + index + width] = (value0 + value2 - value1 - value3) % mod
values[offset + index + 2 * width] = (value0 - value2 + difference) % mod
values[offset + index + 3 * width] = (value0 - value2 - difference) % mod
rotation = rotation * rate3[(~block & -~block).bit_length()] % mod
level += 2
return values
def _convolution_ntt998_butterfly_inv(values):
"""In-place inverse radix-4 transform without division by the length."""
size = len(values)
_convolution_ntt998_check_length(size)
if size == 1:
values[0] %= MOD
return values
height = (size - 1).bit_length()
level = height
mod = MOD
inverse_imag = _convolution_ntt998_IIMAG
irate2 = _convolution_ntt998_IRATE2
irate3 = _convolution_ntt998_IRATE3
while level:
if level == 1:
width = 1 << height - level
rotation = 1
for block in range(1 << level - 1):
offset = block << height - level + 1
for index in range(width):
left = values[offset + index]
right = values[offset + index + width]
values[offset + index] = (left + right) % mod
values[offset + index + width] = (left - right) * rotation % mod
rotation = rotation * irate2[(~block & -~block).bit_length()] % mod
level -= 1
else:
width = 1 << height - level
rotation = 1
for block in range(1 << level - 2):
rotation2 = rotation * rotation % mod
rotation3 = rotation2 * rotation % mod
offset = block << height - level + 2
for index in range(width):
value0 = values[offset + index]
value1 = values[offset + index + width]
value2 = values[offset + index + 2 * width]
value3 = values[offset + index + 3 * width]
difference = (value2 - value3) * inverse_imag % mod
values[offset + index] = (value0 + value1 + value2 + value3) % mod
values[offset + index + width] = (value0 - value1 + difference) * rotation % mod
values[offset + index + 2 * width] = (value0 + value1 - value2 - value3) * rotation2 % mod
values[offset + index + 3 * width] = (value0 - value1 - difference) * rotation3 % mod
rotation = rotation * irate3[(~block & -~block).bit_length()] % mod
level -= 2
return values
def ntt(values):
"""`values`を破壊的に順変換し、同じlistを返す。O(N log N)。"""
return _convolution_ntt998_butterfly(values)
def _convolution_ntt998__intt(values):
size = len(values)
_convolution_ntt998_butterfly_inv(values)
inverse_size = _convolution_ntt998_INVERSE_SIZE.get(size)
if inverse_size is None:
inverse_size = pow(size, MOD - 2, MOD)
_convolution_ntt998_INVERSE_SIZE[size] = inverse_size
for index in range(size):
values[index] = values[index] * inverse_size % MOD
return values
def intt(values):
"""`values`を破壊的に正規化済み逆変換し、同じlistを返す。O(N log N)。"""
return _convolution_ntt998__intt(values)
def _convolution_ntt998_multiply_naive(first, second):
first_size = len(first)
second_size = len(second)
if first_size == 0 or second_size == 0:
return []
if first_size < second_size:
(first, second) = (second, first)
(first_size, second_size) = (second_size, first_size)
result = [0] * (first_size + second_size - 1)
mod = MOD
for (index, left) in enumerate(first):
left %= mod
if left:
for (offset, right) in enumerate(second):
result[index + offset] += left * right
if index & 7 == 7:
start = index
stop = min(index + second_size, len(result))
for position in range(start, stop):
result[position] %= mod
return [value % mod for value in result]
def multiply(first, second):
"""2つの係数列の積を長さ`len(first)+len(second)-1`で返す。O(N log N)。"""
first_size = len(first)
second_size = len(second)
if first_size == 0 or second_size == 0:
return []
if min(first_size, second_size) <= 60:
return _convolution_ntt998_multiply_naive(first, second)
output_size = first_size + second_size - 1
size = 1 << (output_size - 1).bit_length()
_convolution_ntt998_check_length(size)
left = [value % MOD for value in first]
left.extend([0] * (size - first_size))
_convolution_ntt998_butterfly(left)
if first is second:
for index in range(size):
left[index] = left[index] * left[index] % MOD
else:
right = [value % MOD for value in second]
right.extend([0] * (size - second_size))
_convolution_ntt998_butterfly(right)
for index in range(size):
left[index] = left[index] * right[index] % MOD
_convolution_ntt998_butterfly_inv(left)
inverse_size = _convolution_ntt998_INVERSE_SIZE.get(size)
if inverse_size is None:
inverse_size = pow(size, MOD - 2, MOD)
_convolution_ntt998_INVERSE_SIZE[size] = inverse_size
for index in range(output_size):
left[index] = left[index] * inverse_size % MOD
del left[output_size:]
return left
def square(series):
"""係数列の二乗を長さ`2*len(series)-1`で返す。O(N log N)。"""
size = len(series)
if size == 0:
return []
if size <= 60:
result = [0] * (2 * size - 1)
for (index, left) in enumerate(series):
left %= MOD
result[index << 1] += left * left
for offset in range(index + 1, size):
result[index + offset] += 2 * left * series[offset]
return [value % MOD for value in result]
return multiply(series, series)
# 依存: fps998/FPS.py
"""998244353上の形式的冪級数を昇べき順の係数listで計算する。
`a[i]`は$x^i$の係数を表す。inv・log・exp・pow・sqrt、微分・積分、
多項式除算、Taylor shift、一括積を固定modのradix-4 NTTで計算する。
入力listは変更せず、新しい係数listを返す。
"""
from heapq import heapify, heappop, heappush
_fps998_fps_INVERSES = [0, 1]
def _fps998_fps_mod_sqrt(value):
value %= MOD
if value < 2:
return value
if pow(value, MOD - 1 >> 1, MOD) != 1:
return -1
odd = MOD - 1
exponent = 0
while odd & 1 == 0:
odd >>= 1
exponent += 1
nonresidue = 3
root = pow(value, odd + 1 >> 1, MOD)
remainder = pow(value, odd, MOD)
generator = pow(nonresidue, odd, MOD)
level = exponent
while remainder != 1:
position = 1
squared = remainder * remainder % MOD
while position < level and squared != 1:
squared = squared * squared % MOD
position += 1
if position == level:
return -1
adjustment = pow(generator, 1 << level - position - 1, MOD)
root = root * adjustment % MOD
adjustment = adjustment * adjustment % MOD
remainder = remainder * adjustment % MOD
generator = adjustment
level = position
return min(root, MOD - root)
def _fps998_fps_degree(degree, default):
if degree is None:
return default
if degree < 0:
raise ValueError('degree must be nonnegative')
return degree
def _fps998_fps_inverses(size):
if size >= MOD:
raise ValueError('formal integration requires degree < 998244353')
values = _fps998_fps_INVERSES
for index in range(len(values), size + 1):
values.append(-values[MOD % index] * (MOD // index) % MOD)
return values
def shrink(series):
"""係数をmodで正規化し、末尾の0を除いた新しいlistを返す。O(N)。"""
result = [value % MOD for value in series]
while result and result[-1] == 0:
result.pop()
return result
def fps_add(first, second):
"""2つのFPSを係数ごとに加え、長い方と同じ長さのlistを返す。O(N)。"""
size = max(len(first), len(second))
result = [0] * size
common = min(len(first), len(second))
mod = MOD
for index in range(common):
result[index] = (first[index] + second[index]) % mod
for index in range(common, len(first)):
result[index] = first[index] % mod
for index in range(common, len(second)):
result[index] = second[index] % mod
return result
def fps_sub(first, second):
"""`first-second`の係数列を長い方と同じ長さで返す。O(N)。"""
size = max(len(first), len(second))
result = [0] * size
common = min(len(first), len(second))
mod = MOD
for index in range(common):
result[index] = (first[index] - second[index]) % mod
for index in range(common, len(first)):
result[index] = first[index] % mod
for index in range(common, len(second)):
result[index] = -second[index] % mod
return result
def fps_neg(series):
"""各係数の加法逆元を同じ長さのlistで返す。O(N)。"""
return [-value % MOD for value in series]
def fps_diff(series):
"""形式微分の係数を昇べき順で返す。O(N)。"""
mod = MOD
return [index * series[index] % mod for index in range(1, len(series))]
def fps_integral(series):
"""定数項を0とした形式積分の係数を昇べき順で返す。O(N)。"""
inverse = _fps998_fps_inverses(len(series))
mod = MOD
result = [0] * (len(series) + 1)
for (index, value) in enumerate(series, 1):
result[index] = value * inverse[index] % mod
return result
def fps_eval(series, value):
"""FPSを多項式とみなし`value`へ代入した値を返す。O(N)。"""
result = 0
value %= MOD
mod = MOD
for coefficient in reversed(series):
result = (result * value + coefficient) % mod
return result
def _fps998_fps_inverse_step(series, result, current, target):
size = current << 1
mod = MOD
left = [value % mod for value in series[:target]]
left.extend([0] * (size - len(left)))
right = result + [0] * (size - current)
_convolution_ntt998_butterfly(left)
_convolution_ntt998_butterfly(right)
for index in range(size):
left[index] = left[index] * right[index] % mod
_convolution_ntt998__intt(left)
for index in range(current):
left[index] = 0
for index in range(current, target):
left[index] = -left[index] % mod
for index in range(target, size):
left[index] = 0
_convolution_ntt998_butterfly(left)
for index in range(size):
left[index] = left[index] * right[index] % mod
_convolution_ntt998__intt(left)
result.extend(left[current:target])
def fps_inv(series, degree=None):
"""$1/f(x)\\bmod x^{degree}$の係数を`degree`個返す。O(N log N)。"""
degree = _fps998_fps_degree(degree, len(series))
if degree == 0:
return []
if not series or series[0] % MOD == 0:
raise ZeroDivisionError('fps inverse requires nonzero constant coefficient')
result = [pow(series[0] % MOD, MOD - 2, MOD)]
current = 1
while current < degree:
target = min(current << 1, degree)
_convolution_ntt998_check_length(current << 1)
_fps998_fps_inverse_step(series, result, current, target)
current = target
return result
def fps_div(numerator, denominator, degree=None):
"""Return ``numerator / denominator mod x^degree``. O(N log N)."""
degree = _fps998_fps_degree(degree, len(numerator))
if degree == 0:
return []
inverse = fps_inv(denominator, degree)
result = multiply(numerator[:degree], inverse)[:degree]
result.extend([0] * (degree - len(result)))
return result
def fps_log(series, degree=None):
"""$\\log f(x)\\bmod x^{degree}$を返す。`f[0]`は1。O(N log N)。"""
degree = _fps998_fps_degree(degree, len(series))
if degree == 0:
return []
if not series or series[0] % MOD != 1:
raise ValueError('fps logarithm requires constant coefficient 1')
product = multiply(fps_diff(series), fps_inv(series, degree))
result = fps_integral(product[:degree - 1])
result.extend([0] * (degree - len(result)))
return result
def _fps998_fps_fps_exp_ntt(series, degree):
mod = MOD
b = [1, series[1] % mod if len(series) > 1 else 0]
c = [1]
z2 = [1, 1]
inverse = [0, 1]
size = 2
while size < degree:
doubled = size << 1
y = b + [0] * size
_convolution_ntt998_butterfly(y)
z1 = z2
z = [y[index] * z1[index] % mod for index in range(size)]
_convolution_ntt998__intt(z)
for index in range(size >> 1):
z[index] = 0
_convolution_ntt998_butterfly(z)
for index in range(size):
z[index] = -z[index] * z1[index] % mod
_convolution_ntt998__intt(z)
c.extend(z[size >> 1:])
z2 = c + [0] * size
_convolution_ntt998_butterfly(z2)
source_size = min(len(series), size)
x = [series[index] % mod for index in range(source_size)]
x.extend([0] * (size - source_size))
x = fps_diff(x)
x.append(0)
_convolution_ntt998_butterfly(x)
for index in range(size):
x[index] = x[index] * y[index] % mod
_convolution_ntt998__intt(x)
for index in range(1, len(b)):
x[index - 1] = (x[index - 1] - index * b[index]) % mod
x.extend([0] * size)
for index in range(size - 1):
(x[size + index], x[index]) = (x[index], 0)
_convolution_ntt998_butterfly(x)
for index in range(doubled):
x[index] = x[index] * z2[index] % mod
_convolution_ntt998__intt(x)
x.pop()
for index in range(len(inverse), len(x) + 1):
inverse.append(-inverse[mod % index] * (mod // index) % mod)
x = [0] + [value * inverse[index + 1] % mod for (index, value) in enumerate(x)]
for index in range(size):
x[index] = 0
for index in range(size, min(len(series), doubled)):
x[index] = (x[index] + series[index]) % mod
_convolution_ntt998_butterfly(x)
for index in range(doubled):
x[index] = x[index] * y[index] % mod
_convolution_ntt998__intt(x)
b.extend(x[size:])
size = doubled
return b[:degree]
def fps_exp(series, degree=None):
"""$\\exp f(x)\\bmod x^{degree}$を返す。`f[0]`は0。O(N log N)。"""
degree = _fps998_fps_degree(degree, len(series))
if degree == 0:
return []
if series and series[0] % MOD:
raise ValueError('fps exponential requires constant coefficient 0')
_convolution_ntt998_check_length(1 << (degree - 1).bit_length())
if degree == 1:
return [1]
return _fps998_fps_fps_exp_ntt(series, degree)
def fps_pow(series, exponent, degree=None):
"""$f(x)^{exponent}\\bmod x^{degree}$を係数`degree`個で返す。O(N log N)。"""
degree = _fps998_fps_degree(degree, len(series))
if degree == 0:
return []
if exponent == 0:
return [1] + [0] * (degree - 1)
leading = 0
while leading < len(series) and series[leading] % MOD == 0:
leading += 1
if leading == len(series):
if exponent < 0:
raise ZeroDivisionError('a zero series cannot have negative exponent')
return [0] * degree
if exponent < 0 and leading:
raise ValueError('negative power requires an invertible series')
shift = leading * exponent
if shift >= degree:
return [0] * degree
coefficient = series[leading] % MOD
inverse_coefficient = pow(coefficient, MOD - 2, MOD)
needed = degree - shift
normalized = [value * inverse_coefficient % MOD for value in series[leading:]]
logarithm = fps_log(normalized, needed)
for index in range(needed):
logarithm[index] = logarithm[index] * exponent % MOD
result = fps_exp(logarithm, needed)
scale = pow(coefficient, exponent, MOD)
return [0] * shift + [value * scale % MOD for value in result]
def fps_sqrt(series, degree=None):
"""$g(x)^2=f(x)\\bmod x^{degree}$となる係数列を返し、なければ`None`。O(N log N)。"""
degree = _fps998_fps_degree(degree, len(series))
if degree == 0:
return []
leading = 0
limit = min(len(series), degree)
while leading < limit and series[leading] % MOD == 0:
leading += 1
if leading == limit:
return [0] * degree
if leading & 1:
return None
shift = leading >> 1
needed = degree - shift
source = [value % MOD for value in series[leading:]]
root = _fps998_fps_mod_sqrt(source[0])
if root == -1:
return None
inverse_two = MOD + 1 >> 1
result = [root]
current = 1
while current < needed:
target = min(current << 1, needed)
quotient = multiply(source[:target], fps_inv(result, target))[:target]
result.extend([0] * (target - len(result)))
for index in range(target):
value = quotient[index] if index < len(quotient) else 0
result[index] = (result[index] + value) * inverse_two % MOD
current = target
return [0] * shift + result[:needed]
def taylor_shift(series, shift):
"""$f(x+shift)$の係数を`f`と同じ長さのlistで返す。O(N log N)。"""
size = len(series)
if size == 0:
return []
factorial = [1] * size
for index in range(1, size):
factorial[index] = factorial[index - 1] * index % MOD
inverse_factorial = [0] * size
inverse_factorial[-1] = pow(factorial[-1], MOD - 2, MOD)
for index in range(size - 1, 0, -1):
inverse_factorial[index - 1] = inverse_factorial[index] * index % MOD
left = [series[index] * factorial[index] % MOD for index in range(size)]
left.reverse()
right = [0] * size
power = 1
shift %= MOD
for index in range(size):
right[index] = power * inverse_factorial[index] % MOD
power = power * shift % MOD
product = multiply(left, right)
return [product[size - 1 - index] * inverse_factorial[index] % MOD for index in range(size)]
def fps_product(polynomials):
"""複数の多項式をすべて掛けた係数列を返す。O(S log S log K)。"""
heap = []
serial = 0
for polynomial in polynomials:
values = [value % MOD for value in polynomial]
if not values:
return []
heap.append((len(values), serial, values))
serial += 1
if not heap:
return [1]
heapify(heap)
while len(heap) > 1:
(_, _, first) = heappop(heap)
(_, _, second) = heappop(heap)
product = multiply(first, second)
heappush(heap, (len(product), serial, product))
serial += 1
return heap[0][2]
# 依存: fps998/PowerProjection.py
"""998244353上で多項式の冪に関する係数をまとめて列挙する。"""
def _fps998_power_projection_power_projection_zero_constant(polynomial, weights, count):
if not polynomial or not weights:
return [0] * count
original_size = len(weights)
size = 1
while size < original_size:
size <<= 1
product_numerator = [0] * (size << 1)
denominator = [0] * (size << 1)
padded_weights = list(weights) + [0] * (size - original_size)
padded_weights.reverse()
product_numerator[:size] = [value % MOD for value in padded_weights]
limit = min(len(polynomial), original_size)
for index in range(1, limit):
denominator[index] = -polynomial[index] % MOD
block_count = 1
while size > 1:
reflected = denominator[:]
for index in range(1, len(reflected), 2):
reflected[index] = -reflected[index] % MOD
next_numerator = multiply(product_numerator, reflected)
next_denominator = multiply(denominator, reflected)
expanded = size * block_count << 2
next_numerator.extend([0] * (expanded - len(next_numerator)))
next_denominator.extend([0] * (expanded - len(next_denominator)))
offset = size * block_count << 1
source_length = size * block_count << 1
for index in range(source_length):
next_numerator[offset + index] = (next_numerator[offset + index] + product_numerator[index]) % MOD
next_denominator[offset + index] = (next_denominator[offset + index] + denominator[index] + reflected[index]) % MOD
new_length = size * block_count << 1
new_numerator = [0] * new_length
new_denominator = [0] * new_length
half = size >> 1
for block in range(block_count << 1):
source = block * size * 2
destination = block * size
for index in range(half):
new_numerator[destination + index] = next_numerator[source + (index << 1) + 1]
new_denominator[destination + index] = next_denominator[source + (index << 1)]
product_numerator = new_numerator
denominator = new_denominator
size >>= 1
block_count <<= 1
result = [product_numerator[index << 1] for index in range(block_count)]
result.reverse()
result = result[:count]
result.extend([0] * (count - len(result)))
return result
def power_projection(polynomial, weights, count):
"""`result[i]=sum_j weights[j][x^j]polynomial(x)^i`を`0<=i<count`で返す。O(N log^2 N)。"""
if count < 0:
raise ValueError('count must be nonnegative')
if count == 0:
return []
if not polynomial or not weights:
return [0] * count
constant = polynomial[0] % MOD
shifted = list(polynomial)
shifted[0] = 0
result = _fps998_power_projection_power_projection_zero_constant(shifted, weights, count)
if constant == 0:
return result
factorial = [1] * count
inverse_factorial = [1] * count
for index in range(1, count):
factorial[index] = factorial[index - 1] * index % MOD
inverse_factorial[-1] = pow(factorial[-1], MOD - 2, MOD)
for index in range(count - 1, 0, -1):
inverse_factorial[index - 1] = inverse_factorial[index] * index % MOD
coefficient = [0] * count
power = 1
for index in range(count):
result[index] = result[index] * inverse_factorial[index] % MOD
coefficient[index] = inverse_factorial[index] * power % MOD
power = power * constant % MOD
result = multiply(result, coefficient)[:count]
result.extend([0] * (count - len(result)))
return [result[index] * factorial[index] % MOD for index in range(count)]
def power_coefficient(polynomial, multiplier=None, count=None):
"""`[x^n]polynomial(x)^i multiplier(x)`を`0<=i<count`で返す。O(N log^2 N)。"""
degree = len(polynomial) - 1
if degree < 0:
return [] if count is None else [0] * count
if multiplier is None:
multiplier = [1]
if count is None:
count = degree + 1
weights = [0] * (degree + 1)
for exponent in range(degree + 1):
multiplier_index = degree - exponent
if multiplier_index < len(multiplier):
weights[exponent] = multiplier[multiplier_index]
return power_projection(polynomial, weights, count)
# 本体
"""998244353上でFPS合成と合成逆関数を計算する。
`fps_compose(outer, inner, degree)`は`outer(inner(x)) mod x^degree`、
`fps_compositional_inv(series, degree)`は`series(g(x))=x mod x^degree`
となる`g`の係数列を返す。
"""
def _add_constant(series, value):
if series:
series[0] = (series[0] + value) % MOD
else:
series.append(value % MOD)
def _compose_naive(outer, inner, degree):
result = []
inner = [value % MOD for value in inner[:degree]]
for coefficient in reversed(outer[:degree]):
result = multiply(result, inner)[:degree]
_add_constant(result, coefficient)
result.extend([0] * (degree - len(result)))
return result
def _build_frequency_q(series, n, height, blocks):
total = 4 * height * blocks
frequency = [0] * total
for block in range(blocks):
source = block * height
target = block * height * 2
frequency[target:target + n + 1] = series[source:source + n + 1]
frequency[blocks * height * 2] = (frequency[blocks * height * 2] + 1) % MOD
_convolution_ntt998_butterfly(frequency)
for index in range(0, total, 2):
(frequency[index], frequency[index + 1]) = (frequency[index + 1], frequency[index])
return frequency
def _descend_q(series, n, height, blocks):
frequency = _build_frequency_q(series, n, height, blocks)
half_total = 2 * height * blocks
reduced = [0] * half_total
for index in range(half_total):
reduced[index] = frequency[index << 1] * frequency[index << 1 | 1] % MOD
_convolution_ntt998__intt(reduced)
reduced[0] = (reduced[0] - 1) % MOD
child_height = height >> 1
child = [0] * (height * blocks)
child_degree = n >> 1
for block in range(blocks << 1):
source = block * height
target = block * child_height
child[target:target + child_degree + 1] = reduced[source:source + child_degree + 1]
return child
def _reverse_frequency_blocks(values):
start = 1
while start < len(values):
left = start
right = (start << 1) - 1
while left < right:
(values[left], values[right]) = (values[right], values[left])
left += 1
right -= 1
start <<= 1
def _ascend_p(child, series, n, height, blocks):
frequency_q = _build_frequency_q(series, n, height, blocks)
total = len(frequency_q)
frequency_p = [0] * total
child_height = height >> 1
child_degree = n >> 1
parity = n & 1
for block in range(blocks << 1):
source = block * child_height
target = block * height * 2 + parity
for index in range(child_degree + 1):
frequency_p[target + (index << 1)] = child[source + index]
_convolution_ntt998_butterfly(frequency_p)
_reverse_frequency_blocks(frequency_q)
for index in range(total):
frequency_p[index] = frequency_p[index] * frequency_q[index] % MOD
_convolution_ntt998__intt(frequency_p)
result = [0] * (height * blocks)
for block in range(blocks):
source = block * height * 2
target = block * height
result[target:target + n + 1] = frequency_p[source:source + n + 1]
return result
def _compose_ntt(outer, inner, degree):
original_degree = degree - 1
height = 1 << (degree - 1).bit_length()
_convolution_ntt998_check_length(height << 2)
fixed_size = height
outer_values = [value % MOD for value in outer[:degree]]
outer_values.extend([0] * (degree - len(outer_values)))
current = [0] * fixed_size
for (index, value) in enumerate(inner[:degree]):
current[index] = -value % MOD
frames = []
n = original_degree
block_height = height
blocks = 1
while n:
frames.append((current, n, block_height, blocks))
current = _descend_q(current, n, block_height, blocks)
n >>= 1
block_height >>= 1
blocks <<= 1
denominator = current[:blocks]
denominator.append(1)
denominator.reverse()
inverse = fps_inv(denominator, len(denominator))
inverse.reverse()
product = multiply(outer_values, inverse)
result = [0] * fixed_size
for index in range(degree):
result[blocks - 1 - index] = product[index + blocks]
while frames:
(series, n, block_height, blocks) = frames.pop()
result = _ascend_p(result, series, n, block_height, blocks)
result = result[:degree]
result.reverse()
return result
def fps_compose(outer, inner, degree=None):
"""`outer(inner(x)) mod x^degree`の係数を`degree`個返す。O(N log^2 N)。"""
if degree is None:
degree = max(len(outer), len(inner))
if degree < 0:
raise ValueError('degree must be nonnegative')
if degree == 0:
return []
if not outer:
return [0] * degree
if len(outer) == 1:
return [outer[0] % MOD] + [0] * (degree - 1)
if len(inner) <= 1:
point = inner[0] % MOD if inner else 0
value = 0
for coefficient in reversed(outer[:degree]):
value = (value * point + coefficient) % MOD
return [value] + [0] * (degree - 1)
if inner[0] % MOD == 0 and inner[1] % MOD == 1 and (len(inner) == 2):
result = [value % MOD for value in outer[:degree]]
result.extend([0] * (degree - len(result)))
return result
if degree <= 64:
return _compose_naive(outer, inner, degree)
return _compose_ntt(outer, inner, degree)
def fps_compositional_inv(series, degree=None):
"""`series(g(x))=x mod x^degree`となる`g`の係数を返す。O(N log^2 N)。"""
if degree is None:
degree = len(series)
if degree < 0:
raise ValueError('degree must be nonnegative')
if degree == 0:
return []
if not series or series[0] % MOD:
raise ValueError('compositional inverse requires series[0] = 0')
if len(series) < 2 or series[1] % MOD == 0:
raise ValueError('a nonzero linear coefficient is required')
inverse_linear = pow(series[1] % MOD, MOD - 2, MOD)
if degree == 1:
return [0]
last = min(len(series), degree) - 1
while last > 1 and series[last] % MOD == 0:
last -= 1
if last == 1:
return [0, inverse_linear] + [0] * (degree - 2)
order = degree - 1
source = [value % MOD for value in series[:degree]]
source.extend([0] * (degree - len(source)))
coefficients = power_coefficient(source, count=degree)
inverses = [0] * degree
inverses[1] = 1
for index in range(2, degree):
inverses[index] = -(MOD // index) * inverses[MOD % index] % MOD
for index in range(1, degree):
coefficients[index] = coefficients[index] * order * inverses[index] % MOD
coefficients.reverse()
scale = pow(coefficients[0], MOD - 2, MOD)
for index in range(degree):
coefficients[index] = coefficients[index] * scale % MOD
exponent = -pow(order, MOD - 2, MOD) % MOD
logarithm = fps_log(coefficients, degree - 1)
for index in range(degree - 1):
logarithm[index] = logarithm[index] * exponent % MOD
result = fps_exp(logarithm, degree - 1)
for index in range(degree - 1):
result[index] = result[index] * inverse_linear % MOD
return [0] + result
n, m = map(int, input().split())
f = list(map(int, input().split()))
g = list(map(int, input().split()))
ig = list(map(int, input().split()))
a1 = pow(f[1], m, MOD)
ans = fps_compose(ig, [x*a1%MOD for x in g])
print(*ans)