結果
| 問題 | No.3620 Compositional Power with Schröder Coordinate 2 |
| コンテスト | |
| ユーザー |
 lif4635
|
| 提出日時 | 2026-08-10 06:49:46 |
| 言語 | PyPy3 (7.3.17) |
| 結果 |
AC
|
| 実行時間 | 5,057 ms / 10,000 ms |
| + 755µs | |
| コード長 | 41,109 bytes |
| 記録 | |
| コンパイル時間 | 245 ms |
| コンパイル使用メモリ | 95,852 KB |
| 実行使用メモリ | 317,060 KB |
| 最終ジャッジ日時 | 2026-08-10 23:59:30 |
| 合計ジャッジ時間 | 27,664 ms |
|
ジャッジサーバーID （参考情報） |
judge1_1 / judge3_1 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 2 |
| other | AC * 7 |
ソースコード
"""
power projection にも ラグランジュ反転 にも帰着できなかったので変形こねこね
h = g^-1
b_k = \sum_i [x^n] h_i (a^k g)^i
[x^n] g^i が i に対して列挙できれば良い
これは powerprojection
"""
# 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]
_fps998_fps_SPARSE_INV_THRESHOLD = 160
_fps998_fps_SPARSE_DIV_THRESHOLD = 200
_fps998_fps_SPARSE_LOG_THRESHOLD = 200
_fps998_fps_SPARSE_EXP_THRESHOLD = 320
_fps998_fps_SPARSE_POWER_THRESHOLD = 32
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 _fps998_fps_sparse_terms(series, degree, threshold):
terms = []
for index in range(1, min(len(series), degree)):
value = series[index] % MOD
if value:
terms.append((index, value))
if len(terms) > threshold:
return None
return terms
def _fps998_fps_fps_inv_sparse(series, degree, first_inverse, terms):
result = [0] * degree
result[0] = first_inverse
for index in range(1, degree):
total = 0
for offset, coefficient in terms:
if offset > index:
break
total += coefficient * result[index - offset]
result[index] = -total * first_inverse % MOD
return result
def _fps998_fps_fps_div_sparse(numerator, degree, first_inverse, terms):
result = [0] * degree
for index in range(degree):
total = numerator[index] % MOD if index < len(numerator) else 0
for offset, coefficient in terms:
if offset > index:
break
total -= coefficient * result[index - offset]
result[index] = total * first_inverse % MOD
return result
def _fps998_fps_fps_log_sparse(series, degree, terms):
inverse = _fps998_fps_inverses(degree)
result = [0] * degree
for index in range(1, degree):
total = index * (series[index] % MOD) if index < len(series) else 0
for offset, coefficient in terms:
if offset >= index:
break
total -= (index - offset) * coefficient * result[index - offset]
result[index] = total * inverse[index] % MOD
return result
def _fps998_fps_fps_exp_sparse(degree, terms):
inverse = _fps998_fps_inverses(degree)
result = [0] * degree
result[0] = 1
for index in range(1, degree):
total = 0
for offset, coefficient in terms:
if offset > index:
break
total += offset * coefficient * result[index - offset]
result[index] = total * inverse[index] % MOD
return result
def _fps998_fps_fps_power_unit_sparse(degree, exponent, terms):
inverse = _fps998_fps_inverses(degree)
result = [0] * degree
result[0] = 1
exponent %= MOD
for index in range(1, degree):
total = 0
for offset, coefficient in terms:
if offset > index:
break
factor = (exponent * offset - index + offset) % MOD
total += factor * coefficient * result[index - offset]
result[index] = total * inverse[index] % MOD
return result
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')
first_inverse = pow(series[0] % MOD, MOD - 2, MOD)
terms = _fps998_fps_sparse_terms(series, degree, _fps998_fps_SPARSE_INV_THRESHOLD)
if terms is not None:
return _fps998_fps_fps_inv_sparse(series, degree, first_inverse, terms)
result = [first_inverse]
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 []
if not denominator or denominator[0] % MOD == 0:
raise ZeroDivisionError('fps division requires nonzero denominator constant')
first_inverse = pow(denominator[0] % MOD, MOD - 2, MOD)
terms = _fps998_fps_sparse_terms(denominator, degree, _fps998_fps_SPARSE_DIV_THRESHOLD)
if terms is not None:
return _fps998_fps_fps_div_sparse(numerator, degree, first_inverse, terms)
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')
terms = _fps998_fps_sparse_terms(series, degree, _fps998_fps_SPARSE_LOG_THRESHOLD)
if terms is not None:
return _fps998_fps_fps_log_sparse(series, degree, terms)
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')
terms = _fps998_fps_sparse_terms(series, degree, _fps998_fps_SPARSE_EXP_THRESHOLD)
if terms is not None:
return _fps998_fps_fps_exp_sparse(degree, terms)
_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:]]
terms = _fps998_fps_sparse_terms(normalized, needed, _fps998_fps_SPARSE_POWER_THRESHOLD)
if terms is not None:
result = _fps998_fps_fps_power_unit_sparse(needed, exponent, terms)
else:
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_constant = pow(source[0], MOD - 2, MOD)
normalized = [value * inverse_constant % MOD for value in source]
terms = _fps998_fps_sparse_terms(normalized, needed, _fps998_fps_SPARSE_POWER_THRESHOLD)
if terms is not None:
result = _fps998_fps_fps_power_unit_sparse(needed, MOD + 1 >> 1, terms)
return [0] * shift + [value * root % MOD for value in result]
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
degree = original_size - 1
height = size
blocks = 1
numerator = [0] * height
denominator = [0] * height
padded_weights = [value % MOD for value in reversed(weights)]
padded_weights.extend([0] * (size - original_size))
numerator[:size] = padded_weights
limit = min(len(polynomial), original_size)
for index in range(1, limit):
denominator[index] = -polynomial[index] % MOD
while degree:
total = 4 * height * blocks
frequency_p = [0] * total
frequency_q = [0] * total
for block in range(blocks):
source = block * height
target = block * height * 2
frequency_p[target:target + degree + 1] = numerator[source:source + degree + 1]
frequency_q[target:target + degree + 1] = denominator[source:source + degree + 1]
frequency_q[blocks * height * 2] = (frequency_q[blocks * height * 2] + 1) % MOD
_convolution_ntt998_butterfly(frequency_p)
_convolution_ntt998_butterfly(frequency_q)
reduced_q = [0] * (total >> 1)
for index in range(0, total, 2):
frequency_q[index], frequency_q[index + 1] = (frequency_q[index + 1], frequency_q[index])
frequency_p[index] = frequency_p[index] * frequency_q[index] % MOD
frequency_p[index + 1] = frequency_p[index + 1] * frequency_q[index + 1] % MOD
reduced_q[index >> 1] = frequency_q[index] * frequency_q[index + 1] % MOD
_convolution_ntt998__intt(frequency_p)
_convolution_ntt998__intt(reduced_q)
reduced_q[0] = (reduced_q[0] - 1) % MOD
child_height = height >> 1
child_degree = degree >> 1
parity = degree & 1
child_size = height * blocks
child_p = [0] * child_size
child_q = [0] * child_size
for block in range(blocks << 1):
source = block * height * 2
target = block * child_height
for index in range(child_degree + 1):
child_p[target + index] = frequency_p[source + (index << 1) + parity]
child_q[target + index] = reduced_q[block * height + index]
numerator = child_p
denominator = child_q
degree >>= 1
height >>= 1
blocks <<= 1
result = numerator[:blocks]
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`の係数列を返す。
"""
from array import array
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, array('I', frequency))
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, frequency_q, 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:
current, frequency_q = _descend_q(current, n, block_height, blocks)
frames.append((frequency_q, 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:
frequency_q, n, block_height, blocks = frames.pop()
result = _ascend_p(result, frequency_q, 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
mod = MOD
def chirp_z(polynomial, ratio, count=None, start=1):
"""Return f(start*ratio^i) for i=0..count-1 by Bluestein."""
polynomial = [value % mod for value in polynomial]
if count is None:
count = len(polynomial)
if count < 0:
raise ValueError("count must be nonnegative")
if not polynomial or count == 0:
return [0] * count
if start % mod != 1:
power = 1
for i in range(len(polynomial)):
polynomial[i] = polynomial[i] * power % mod
power = power * start % mod
ratio %= mod
if ratio == 0:
result = [polynomial[0]] * count
result[0] = sum(polynomial) % mod
return result
length = len(polynomial)
triangular = [1] * (count + length)
inverse_triangular = [1] * max(count, length)
step = 1
for i in range(1, len(triangular)):
triangular[i] = triangular[i - 1] * step % mod
step = step * ratio % mod
inverse_ratio = pow(ratio, -1, mod)
step = 1
for i in range(1, len(inverse_triangular)):
inverse_triangular[i] = inverse_triangular[i - 1] * step % mod
step = step * inverse_ratio % mod
for i in range(length):
polynomial[i] = polynomial[i] * inverse_triangular[i] % mod
polynomial.reverse()
product = multiply(polynomial, triangular)
return [product[length - 1 + i] * inverse_triangular[i] % mod
for i in range(count)]
n, m = map(int, input().split())
f = list(map(int, input().split()))
g = list(map(int, input().split()))
ig = list(map(int, input().split()))
t = power_projection(g, [0]*(n-1) + [1], n)
t = [t[i]*ig[i]%MOD for i in range(n)]
ans = chirp_z(t, f[1], m)
print(*ans)