QOJ.ac
QOJ
ID | Problem | Submitter | Result | Time | Memory | Language | File size | Submit time | Judge time |
---|---|---|---|---|---|---|---|---|---|
#667554 | #9493. 路径计数 | hos_lyric# | 26 | 773ms | 38840kb | C++14 | 43.3kb | 2024-10-23 00:08:35 | 2024-10-23 00:08:44 |
Judging History
answer
#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <bitset>
#include <complex>
#include <deque>
#include <functional>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <sstream>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
using namespace std;
using Int = long long;
template <class T1, class T2> ostream &operator<<(ostream &os, const pair<T1, T2> &a) { return os << "(" << a.first << ", " << a.second << ")"; };
template <class T> ostream &operator<<(ostream &os, const vector<T> &as) { const int sz = as.size(); os << "["; for (int i = 0; i < sz; ++i) { if (i >= 256) { os << ", ..."; break; } if (i > 0) { os << ", "; } os << as[i]; } return os << "]"; }
template <class T> void pv(T a, T b) { for (T i = a; i != b; ++i) cerr << *i << " "; cerr << endl; }
template <class T> bool chmin(T &t, const T &f) { if (t > f) { t = f; return true; } return false; }
template <class T> bool chmax(T &t, const T &f) { if (t < f) { t = f; return true; } return false; }
#define COLOR(s) ("\x1b[" s "m")
////////////////////////////////////////////////////////////////////////////////
template <unsigned M_> struct ModInt {
static constexpr unsigned M = M_;
unsigned x;
constexpr ModInt() : x(0U) {}
constexpr ModInt(unsigned x_) : x(x_ % M) {}
constexpr ModInt(unsigned long long x_) : x(x_ % M) {}
constexpr ModInt(int x_) : x(((x_ %= static_cast<int>(M)) < 0) ? (x_ + static_cast<int>(M)) : x_) {}
constexpr ModInt(long long x_) : x(((x_ %= static_cast<long long>(M)) < 0) ? (x_ + static_cast<long long>(M)) : x_) {}
ModInt &operator+=(const ModInt &a) { x = ((x += a.x) >= M) ? (x - M) : x; return *this; }
ModInt &operator-=(const ModInt &a) { x = ((x -= a.x) >= M) ? (x + M) : x; return *this; }
ModInt &operator*=(const ModInt &a) { x = (static_cast<unsigned long long>(x) * a.x) % M; return *this; }
ModInt &operator/=(const ModInt &a) { return (*this *= a.inv()); }
ModInt pow(long long e) const {
if (e < 0) return inv().pow(-e);
ModInt a = *this, b = 1U; for (; e; e >>= 1) { if (e & 1) b *= a; a *= a; } return b;
}
ModInt inv() const {
unsigned a = M, b = x; int y = 0, z = 1;
for (; b; ) { const unsigned q = a / b; const unsigned c = a - q * b; a = b; b = c; const int w = y - static_cast<int>(q) * z; y = z; z = w; }
assert(a == 1U); return ModInt(y);
}
ModInt operator+() const { return *this; }
ModInt operator-() const { ModInt a; a.x = x ? (M - x) : 0U; return a; }
ModInt operator+(const ModInt &a) const { return (ModInt(*this) += a); }
ModInt operator-(const ModInt &a) const { return (ModInt(*this) -= a); }
ModInt operator*(const ModInt &a) const { return (ModInt(*this) *= a); }
ModInt operator/(const ModInt &a) const { return (ModInt(*this) /= a); }
template <class T> friend ModInt operator+(T a, const ModInt &b) { return (ModInt(a) += b); }
template <class T> friend ModInt operator-(T a, const ModInt &b) { return (ModInt(a) -= b); }
template <class T> friend ModInt operator*(T a, const ModInt &b) { return (ModInt(a) *= b); }
template <class T> friend ModInt operator/(T a, const ModInt &b) { return (ModInt(a) /= b); }
explicit operator bool() const { return x; }
bool operator==(const ModInt &a) const { return (x == a.x); }
bool operator!=(const ModInt &a) const { return (x != a.x); }
friend std::ostream &operator<<(std::ostream &os, const ModInt &a) { return os << a.x; }
};
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
constexpr unsigned MO = 998244353U;
constexpr unsigned MO2 = 2U * MO;
constexpr int FFT_MAX = 23;
using Mint = ModInt<MO>;
constexpr Mint FFT_ROOTS[FFT_MAX + 1] = {1U, 998244352U, 911660635U, 372528824U, 929031873U, 452798380U, 922799308U, 781712469U, 476477967U, 166035806U, 258648936U, 584193783U, 63912897U, 350007156U, 666702199U, 968855178U, 629671588U, 24514907U, 996173970U, 363395222U, 565042129U, 733596141U, 267099868U, 15311432U};
constexpr Mint INV_FFT_ROOTS[FFT_MAX + 1] = {1U, 998244352U, 86583718U, 509520358U, 337190230U, 87557064U, 609441965U, 135236158U, 304459705U, 685443576U, 381598368U, 335559352U, 129292727U, 358024708U, 814576206U, 708402881U, 283043518U, 3707709U, 121392023U, 704923114U, 950391366U, 428961804U, 382752275U, 469870224U};
constexpr Mint FFT_RATIOS[FFT_MAX] = {911660635U, 509520358U, 369330050U, 332049552U, 983190778U, 123842337U, 238493703U, 975955924U, 603855026U, 856644456U, 131300601U, 842657263U, 730768835U, 942482514U, 806263778U, 151565301U, 510815449U, 503497456U, 743006876U, 741047443U, 56250497U, 867605899U};
constexpr Mint INV_FFT_RATIOS[FFT_MAX] = {86583718U, 372528824U, 373294451U, 645684063U, 112220581U, 692852209U, 155456985U, 797128860U, 90816748U, 860285882U, 927414960U, 354738543U, 109331171U, 293255632U, 535113200U, 308540755U, 121186627U, 608385704U, 438932459U, 359477183U, 824071951U, 103369235U};
// as[rev(i)] <- \sum_j \zeta^(ij) as[j]
void fft(Mint *as, int n) {
assert(!(n & (n - 1))); assert(1 <= n); assert(n <= 1 << FFT_MAX);
int m = n;
if (m >>= 1) {
for (int i = 0; i < m; ++i) {
const unsigned x = as[i + m].x; // < MO
as[i + m].x = as[i].x + MO - x; // < 2 MO
as[i].x += x; // < 2 MO
}
}
if (m >>= 1) {
Mint prod = 1U;
for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
for (int i = i0; i < i0 + m; ++i) {
const unsigned x = (prod * as[i + m]).x; // < MO
as[i + m].x = as[i].x + MO - x; // < 3 MO
as[i].x += x; // < 3 MO
}
prod *= FFT_RATIOS[__builtin_ctz(++h)];
}
}
for (; m; ) {
if (m >>= 1) {
Mint prod = 1U;
for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
for (int i = i0; i < i0 + m; ++i) {
const unsigned x = (prod * as[i + m]).x; // < MO
as[i + m].x = as[i].x + MO - x; // < 4 MO
as[i].x += x; // < 4 MO
}
prod *= FFT_RATIOS[__builtin_ctz(++h)];
}
}
if (m >>= 1) {
Mint prod = 1U;
for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
for (int i = i0; i < i0 + m; ++i) {
const unsigned x = (prod * as[i + m]).x; // < MO
as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x; // < 2 MO
as[i + m].x = as[i].x + MO - x; // < 3 MO
as[i].x += x; // < 3 MO
}
prod *= FFT_RATIOS[__builtin_ctz(++h)];
}
}
}
for (int i = 0; i < n; ++i) {
as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x; // < 2 MO
as[i].x = (as[i].x >= MO) ? (as[i].x - MO) : as[i].x; // < MO
}
}
// as[i] <- (1/n) \sum_j \zeta^(-ij) as[rev(j)]
void invFft(Mint *as, int n) {
assert(!(n & (n - 1))); assert(1 <= n); assert(n <= 1 << FFT_MAX);
int m = 1;
if (m < n >> 1) {
Mint prod = 1U;
for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
for (int i = i0; i < i0 + m; ++i) {
const unsigned long long y = as[i].x + MO - as[i + m].x; // < 2 MO
as[i].x += as[i + m].x; // < 2 MO
as[i + m].x = (prod.x * y) % MO; // < MO
}
prod *= INV_FFT_RATIOS[__builtin_ctz(++h)];
}
m <<= 1;
}
for (; m < n >> 1; m <<= 1) {
Mint prod = 1U;
for (int h = 0, i0 = 0; i0 < n; i0 += (m << 1)) {
for (int i = i0; i < i0 + (m >> 1); ++i) {
const unsigned long long y = as[i].x + MO2 - as[i + m].x; // < 4 MO
as[i].x += as[i + m].x; // < 4 MO
as[i].x = (as[i].x >= MO2) ? (as[i].x - MO2) : as[i].x; // < 2 MO
as[i + m].x = (prod.x * y) % MO; // < MO
}
for (int i = i0 + (m >> 1); i < i0 + m; ++i) {
const unsigned long long y = as[i].x + MO - as[i + m].x; // < 2 MO
as[i].x += as[i + m].x; // < 2 MO
as[i + m].x = (prod.x * y) % MO; // < MO
}
prod *= INV_FFT_RATIOS[__builtin_ctz(++h)];
}
}
if (m < n) {
for (int i = 0; i < m; ++i) {
const unsigned y = as[i].x + MO2 - as[i + m].x; // < 4 MO
as[i].x += as[i + m].x; // < 4 MO
as[i + m].x = y; // < 4 MO
}
}
const Mint invN = Mint(n).inv();
for (int i = 0; i < n; ++i) {
as[i] *= invN;
}
}
void fft(vector<Mint> &as) {
fft(as.data(), as.size());
}
void invFft(vector<Mint> &as) {
invFft(as.data(), as.size());
}
vector<Mint> convolve(vector<Mint> as, vector<Mint> bs) {
if (as.empty() || bs.empty()) return {};
const int len = as.size() + bs.size() - 1;
int n = 1;
for (; n < len; n <<= 1) {}
as.resize(n); fft(as);
bs.resize(n); fft(bs);
for (int i = 0; i < n; ++i) as[i] *= bs[i];
invFft(as);
as.resize(len);
return as;
}
vector<Mint> square(vector<Mint> as) {
if (as.empty()) return {};
const int len = as.size() + as.size() - 1;
int n = 1;
for (; n < len; n <<= 1) {}
as.resize(n); fft(as);
for (int i = 0; i < n; ++i) as[i] *= as[i];
invFft(as);
as.resize(len);
return as;
}
// m := |as|, n := |bs|
// cs[k] = \sum[i-j=k] as[i] bs[j] (0 <= k <= m-n)
// transpose of ((multiply by bs): K^[0,m-n] -> K^[0,m-1])
vector<Mint> middle(vector<Mint> as, vector<Mint> bs) {
const int m = as.size(), n = bs.size();
assert(m >= n); assert(n >= 1);
int len = 1;
for (; len < m; len <<= 1) {}
as.resize(len, 0);
fft(as);
std::reverse(bs.begin(), bs.end());
bs.resize(len, 0);
fft(bs);
for (int i = 0; i < len; ++i) as[i] *= bs[i];
invFft(as);
as.resize(m);
as.erase(as.begin(), as.begin() + (n - 1));
return as;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// inv: log, exp, pow
// fac: shift
// invFac: shift
constexpr int LIM_INV = 1 << 20; // @
Mint inv[LIM_INV], fac[LIM_INV], invFac[LIM_INV];
struct ModIntPreparator {
ModIntPreparator() {
inv[1] = 1;
for (int i = 2; i < LIM_INV; ++i) inv[i] = -((Mint::M / i) * inv[Mint::M % i]);
fac[0] = 1;
for (int i = 1; i < LIM_INV; ++i) fac[i] = fac[i - 1] * i;
invFac[0] = 1;
for (int i = 1; i < LIM_INV; ++i) invFac[i] = invFac[i - 1] * inv[i];
}
} preparator;
// polyWork0: *, inv, div, divAt, log, exp, pow, sqrt, shift
// polyWork1: inv, div, divAt, log, exp, pow, sqrt, shift
// polyWork2: divAt, exp, pow, sqrt
// polyWork3: exp, pow, sqrt
static constexpr int LIM_POLY = 1 << 20; // @
static_assert(LIM_POLY <= 1 << FFT_MAX, "Poly: LIM_POLY <= 1 << FFT_MAX must hold.");
static Mint polyWork0[LIM_POLY], polyWork1[LIM_POLY], polyWork2[LIM_POLY], polyWork3[LIM_POLY];
struct Poly : public vector<Mint> {
Poly() {}
explicit Poly(int n) : vector<Mint>(n) {}
Poly(const vector<Mint> &vec) : vector<Mint>(vec) {}
Poly(std::initializer_list<Mint> il) : vector<Mint>(il) {}
int size() const { return vector<Mint>::size(); }
Mint at(long long k) const { return (0 <= k && k < size()) ? (*this)[k] : 0U; }
int ord() const { for (int i = 0; i < size(); ++i) if ((*this)[i]) return i; return -1; }
int deg() const { for (int i = size(); --i >= 0; ) if ((*this)[i]) return i; return -1; }
Poly mod(int n) const { return Poly(vector<Mint>(data(), data() + min(n, size()))); }
friend std::ostream &operator<<(std::ostream &os, const Poly &fs) {
os << "[";
for (int i = 0; i < fs.size(); ++i) { if (i > 0) os << ", "; os << fs[i]; }
return os << "]";
}
Poly &operator+=(const Poly &fs) {
if (size() < fs.size()) resize(fs.size());
for (int i = 0; i < fs.size(); ++i) (*this)[i] += fs[i];
return *this;
}
Poly &operator-=(const Poly &fs) {
if (size() < fs.size()) resize(fs.size());
for (int i = 0; i < fs.size(); ++i) (*this)[i] -= fs[i];
return *this;
}
// 3 E(|t| + |f|)
Poly &operator*=(const Poly &fs) {
if (empty() || fs.empty()) return *this = {};
const int nt = size(), nf = fs.size();
int n = 1;
for (; n < nt + nf - 1; n <<= 1) {}
assert(n <= LIM_POLY);
resize(n);
fft(data(), n); // 1 E(n)
memcpy(polyWork0, fs.data(), nf * sizeof(Mint));
memset(polyWork0 + nf, 0, (n - nf) * sizeof(Mint));
fft(polyWork0, n); // 1 E(n)
for (int i = 0; i < n; ++i) (*this)[i] *= polyWork0[i];
invFft(data(), n); // 1 E(n)
resize(nt + nf - 1);
return *this;
}
// 13 E(deg(t) - deg(f) + 1)
// rev(t) = rev(f) rev(q) + x^(deg(t)-deg(f)+1) rev(r)
Poly &operator/=(const Poly &fs) {
const int m = deg(), n = fs.deg();
assert(n != -1);
if (m < n) return *this = {};
Poly tsRev(m - n + 1), fsRev(min(m - n, n) + 1);
for (int i = 0; i <= m - n; ++i) tsRev[i] = (*this)[m - i];
for (int i = 0, i0 = min(m - n, n); i <= i0; ++i) fsRev[i] = fs[n - i];
const Poly qsRev = tsRev.div(fsRev, m - n + 1); // 13 E(m - n + 1)
resize(m - n + 1);
for (int i = 0; i <= m - n; ++i) (*this)[i] = qsRev[m - n - i];
return *this;
}
// 13 E(deg(t) - deg(f) + 1) + 3 E(|t|)
Poly &operator%=(const Poly &fs) {
const Poly qs = *this / fs; // 13 E(deg(t) - deg(f) + 1)
*this -= fs * qs; // 3 E(|t|)
resize(deg() + 1);
return *this;
}
Poly &operator*=(const Mint &a) {
for (int i = 0; i < size(); ++i) (*this)[i] *= a;
return *this;
}
Poly &operator/=(const Mint &a) {
const Mint b = a.inv();
for (int i = 0; i < size(); ++i) (*this)[i] *= b;
return *this;
}
Poly operator+() const { return *this; }
Poly operator-() const {
Poly fs(size());
for (int i = 0; i < size(); ++i) fs[i] = -(*this)[i];
return fs;
}
Poly operator+(const Poly &fs) const { return (Poly(*this) += fs); }
Poly operator-(const Poly &fs) const { return (Poly(*this) -= fs); }
Poly operator*(const Poly &fs) const { return (Poly(*this) *= fs); }
Poly operator/(const Poly &fs) const { return (Poly(*this) /= fs); }
Poly operator%(const Poly &fs) const { return (Poly(*this) %= fs); }
Poly operator*(const Mint &a) const { return (Poly(*this) *= a); }
Poly operator/(const Mint &a) const { return (Poly(*this) /= a); }
friend Poly operator*(const Mint &a, const Poly &fs) { return fs * a; }
// 10 E(n)
// f <- f - (t f - 1) f
Poly inv(int n) const {
assert(!empty()); assert((*this)[0]); assert(1 <= n);
assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= LIM_POLY);
Poly fs(n);
fs[0] = (*this)[0].inv();
for (int m = 1; m < n; m <<= 1) {
memcpy(polyWork0, data(), min(m << 1, size()) * sizeof(Mint));
memset(polyWork0 + min(m << 1, size()), 0, ((m << 1) - min(m << 1, size())) * sizeof(Mint));
fft(polyWork0, m << 1); // 2 E(n)
memcpy(polyWork1, fs.data(), min(m << 1, n) * sizeof(Mint));
memset(polyWork1 + min(m << 1, n), 0, ((m << 1) - min(m << 1, n)) * sizeof(Mint));
fft(polyWork1, m << 1); // 2 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 2 E(n)
memset(polyWork0, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 2 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 2 E(n)
for (int i = m, i0 = min(m << 1, n); i < i0; ++i) fs[i] = -polyWork0[i];
}
return fs;
}
// 9 E(n)
// Need (4 m)-th roots of unity to lift from (mod x^m) to (mod x^(2m)).
// f <- f - (t f - 1) f
// (t f^2) mod ((x^(2m) - 1) (x^m - 1^(1/4)))
/*
Poly inv(int n) const {
assert(!empty()); assert((*this)[0]); assert(1 <= n);
assert(n == 1 || 3 << (31 - __builtin_clz(n - 1)) <= LIM_POLY);
assert(n <= 1 << (FFT_MAX - 1));
Poly fs(n);
fs[0] = (*this)[0].inv();
for (int h = 2, m = 1; m < n; ++h, m <<= 1) {
const Mint a = FFT_ROOTS[h], b = INV_FFT_ROOTS[h];
memcpy(polyWork0, data(), min(m << 1, size()) * sizeof(Mint));
memset(polyWork0 + min(m << 1, size()), 0, ((m << 1) - min(m << 1, size())) * sizeof(Mint));
{
Mint aa = 1;
for (int i = 0; i < m; ++i) { polyWork0[(m << 1) + i] = aa * polyWork0[i]; aa *= a; }
for (int i = 0; i < m; ++i) { polyWork0[(m << 1) + i] += aa * polyWork0[m + i]; aa *= a; }
}
fft(polyWork0, m << 1); // 2 E(n)
fft(polyWork0 + (m << 1), m); // 1 E(n)
memcpy(polyWork1, fs.data(), min(m << 1, n) * sizeof(Mint));
memset(polyWork1 + min(m << 1, n), 0, ((m << 1) - min(m << 1, n)) * sizeof(Mint));
{
Mint aa = 1;
for (int i = 0; i < m; ++i) { polyWork1[(m << 1) + i] = aa * polyWork1[i]; aa *= a; }
for (int i = 0; i < m; ++i) { polyWork1[(m << 1) + i] += aa * polyWork1[m + i]; aa *= a; }
}
fft(polyWork1, m << 1); // 2 E(n)
fft(polyWork1 + (m << 1), m); // 1 E(n)
for (int i = 0; i < (m << 1) + m; ++i) polyWork0[i] *= polyWork1[i] * polyWork1[i];
invFft(polyWork0, m << 1); // 2 E(n)
invFft(polyWork0 + (m << 1), m); // 1 E(n)
// 2 f0 + (-f2), (-f1) + (-f3), 1^(1/4) (-f1) - (-f2) - 1^(1/4) (-f3)
{
Mint bb = 1;
for (int i = 0, i0 = min(m, n - m); i < i0; ++i) {
unsigned x = polyWork0[i].x + (bb * polyWork0[(m << 1) + i]).x + MO2 - (fs[i].x << 1); // < 4 MO
fs[m + i] = Mint(static_cast<unsigned long long>(FFT_ROOTS[2].x) * x) - polyWork0[m + i];
fs[m + i].x = ((fs[m + i].x & 1) ? (fs[m + i].x + MO) : fs[m + i].x) >> 1;
bb *= b;
}
}
}
return fs;
}
*/
// 13 E(n)
// g = (1 / f) mod x^m
// h <- h - (f h - t) g
Poly div(const Poly &fs, int n) const {
assert(!fs.empty()); assert(fs[0]); assert(1 <= n);
if (n == 1) return {at(0) / fs[0]};
// m < n <= 2 m
const int m = 1 << (31 - __builtin_clz(n - 1));
assert(m << 1 <= LIM_POLY);
Poly gs = fs.inv(m); // 5 E(n)
gs.resize(m << 1);
fft(gs.data(), m << 1); // 1 E(n)
if (size()) memcpy(polyWork0, data(), min(m, size()) * sizeof(Mint));
memset(polyWork0 + min(m, size()), 0, ((m << 1) - min(m, size())) * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= gs[i];
invFft(polyWork0, m << 1); // 1 E(n)
Poly hs(n);
memcpy(hs.data(), polyWork0, m * sizeof(Mint));
memset(polyWork0 + m, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
memcpy(polyWork1, fs.data(), min(m << 1, fs.size()) * sizeof(Mint));
memset(polyWork1 + min(m << 1, fs.size()), 0, ((m << 1) - min(m << 1, fs.size())) * sizeof(Mint));
fft(polyWork1, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 1 E(n)
memset(polyWork0, 0, m * sizeof(Mint));
for (int i = m, i0 = min(m << 1, size()); i < i0; ++i) polyWork0[i] -= (*this)[i];
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= gs[i];
invFft(polyWork0, m << 1); // 1 E(n)
for (int i = m; i < n; ++i) hs[i] = -polyWork0[i];
return hs;
}
// (4 (floor(log_2 k) - ceil(log_2 |f|)) + 16) E(|f|) for |t| < |f|
// [x^k] (t(x) / f(x)) = [x^k] ((t(x) f(-x)) / (f(x) f(-x))
// polyWork0: half of (2 m)-th roots of unity, inversed, bit-reversed
Mint divAt(const Poly &fs, long long k) const {
assert(k >= 0);
if (size() >= fs.size()) {
const Poly qs = *this / fs; // 13 E(deg(t) - deg(f) + 1)
Poly rs = *this - fs * qs; // 3 E(|t|)
rs.resize(rs.deg() + 1);
return qs.at(k) + rs.divAt(fs, k);
}
int h = 0, m = 1;
for (; m < fs.size(); ++h, m <<= 1) {}
if (k < m) {
const Poly gs = fs.inv(k + 1); // 10 E(|f|)
Mint sum;
for (int i = 0, i0 = min<int>(k + 1, size()); i < i0; ++i) sum += (*this)[i] * gs[k - i];
return sum;
}
assert(m << 1 <= LIM_POLY);
polyWork0[0] = Mint(2U).inv();
for (int hh = 0; hh < h; ++hh) for (int i = 0; i < 1 << hh; ++i) polyWork0[1 << hh | i] = polyWork0[i] * INV_FFT_ROOTS[hh + 2];
const Mint a = FFT_ROOTS[h + 1];
if (size()) memcpy(polyWork2, data(), size() * sizeof(Mint));
memset(polyWork2 + size(), 0, ((m << 1) - size()) * sizeof(Mint));
fft(polyWork2, m << 1); // 2 E(|f|)
memcpy(polyWork1, fs.data(), fs.size() * sizeof(Mint));
memset(polyWork1 + fs.size(), 0, ((m << 1) - fs.size()) * sizeof(Mint));
fft(polyWork1, m << 1); // 2 E(|f|)
for (; ; ) {
if (k & 1) {
for (int i = 0; i < m; ++i) polyWork2[i] = polyWork0[i] * (polyWork2[i << 1 | 0] * polyWork1[i << 1 | 1] - polyWork2[i << 1 | 1] * polyWork1[i << 1 | 0]);
} else {
for (int i = 0; i < m; ++i) {
polyWork2[i] = polyWork2[i << 1 | 0] * polyWork1[i << 1 | 1] + polyWork2[i << 1 | 1] * polyWork1[i << 1 | 0];
polyWork2[i].x = ((polyWork2[i].x & 1) ? (polyWork2[i].x + MO) : polyWork2[i].x) >> 1;
}
}
for (int i = 0; i < m; ++i) polyWork1[i] = polyWork1[i << 1 | 0] * polyWork1[i << 1 | 1];
if ((k >>= 1) < m) {
invFft(polyWork2, m); // 1 E(|f|)
invFft(polyWork1, m); // 1 E(|f|)
// Poly::inv does not use polyWork2
const Poly gs = Poly(vector<Mint>(polyWork1, polyWork1 + k + 1)).inv(k + 1); // 10 E(|f|)
Mint sum;
for (int i = 0; i <= k; ++i) sum += polyWork2[i] * gs[k - i];
return sum;
}
memcpy(polyWork2 + m, polyWork2, m * sizeof(Mint));
invFft(polyWork2 + m, m); // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
memcpy(polyWork1 + m, polyWork1, m * sizeof(Mint));
invFft(polyWork1 + m, m); // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
Mint aa = 1;
for (int i = m; i < m << 1; ++i) { polyWork2[i] *= aa; polyWork1[i] *= aa; aa *= a; }
fft(polyWork2 + m, m); // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
fft(polyWork1 + m, m); // (floor(log_2 k) - ceil(log_2 |f|)) E(|f|)
}
}
// 13 E(n)
// D log(t) = (D t) / t
Poly log(int n) const {
assert(!empty()); assert((*this)[0].x == 1U); assert(n <= LIM_INV);
Poly fs = mod(n);
for (int i = 0; i < fs.size(); ++i) fs[i] *= i;
fs = fs.div(*this, n);
for (int i = 1; i < n; ++i) fs[i] *= ::inv[i];
return fs;
}
// (16 + 1/2) E(n)
// f = exp(t) mod x^m ==> (D f) / f == D t (mod x^m)
// g = (1 / exp(t)) mod x^m
// f <- f - (log f - t) / (1 / f)
// = f - (I ((D f) / f) - t) f
// == f - (I ((D f) / f + (f g - 1) ((D f) / f - D (t mod x^m))) - t) f (mod x^(2m))
// = f - (I (g (D f - f D (t mod x^m)) + D (t mod x^m)) - t) f
// g <- g - (f g - 1) g
// polyWork1: DFT(f, 2 m), polyWork2: g, polyWork3: DFT(g, 2 m)
Poly exp(int n) const {
assert(!empty()); assert(!(*this)[0]); assert(1 <= n);
assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= min(LIM_INV, LIM_POLY));
if (n == 1) return {1U};
if (n == 2) return {1U, at(1)};
Poly fs(n);
fs[0].x = polyWork1[0].x = polyWork1[1].x = polyWork2[0].x = 1U;
int m;
for (m = 1; m << 1 < n; m <<= 1) {
for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork0[i] = i * (*this)[i];
memset(polyWork0 + min(m, size()), 0, (m - min(m, size())) * sizeof(Mint));
fft(polyWork0, m); // (1/2) E(n)
for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m); // (1/2) E(n)
for (int i = 0; i < m; ++i) polyWork0[i] -= i * fs[i];
memset(polyWork0 + m, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
memcpy(polyWork3, polyWork2, m * sizeof(Mint));
memset(polyWork3 + m, 0, m * sizeof(Mint));
fft(polyWork3, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
invFft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m; ++i) polyWork0[i] *= ::inv[m + i];
for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork0[i] += (*this)[m + i];
memset(polyWork0 + m, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 1 E(n)
memcpy(fs.data() + m, polyWork0, m * sizeof(Mint));
memcpy(polyWork1, fs.data(), (m << 1) * sizeof(Mint));
memset(polyWork1 + (m << 1), 0, (m << 1) * sizeof(Mint));
fft(polyWork1, m << 2); // 2 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] = polyWork1[i] * polyWork3[i];
invFft(polyWork0, m << 1); // 1 E(n)
memset(polyWork0, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
invFft(polyWork0, m << 1); // 1 E(n)
for (int i = m; i < m << 1; ++i) polyWork2[i] = -polyWork0[i];
}
for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork0[i] = i * (*this)[i];
memset(polyWork0 + min(m, size()), 0, (m - min(m, size())) * sizeof(Mint));
fft(polyWork0, m); // (1/2) E(n)
for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m); // (1/2) E(n)
for (int i = 0; i < m; ++i) polyWork0[i] -= i * fs[i];
memcpy(polyWork0 + m, polyWork0 + (m >> 1), (m >> 1) * sizeof(Mint));
memset(polyWork0 + (m >> 1), 0, (m >> 1) * sizeof(Mint));
memset(polyWork0 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
fft(polyWork0, m); // (1/2) E(n)
fft(polyWork0 + m, m); // (1/2) E(n)
memcpy(polyWork3 + m, polyWork2 + (m >> 1), (m >> 1) * sizeof(Mint));
memset(polyWork3 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
fft(polyWork3 + m, m); // (1/2) E(n)
for (int i = 0; i < m; ++i) polyWork0[m + i] = polyWork0[i] * polyWork3[m + i] + polyWork0[m + i] * polyWork3[i];
for (int i = 0; i < m; ++i) polyWork0[i] *= polyWork3[i];
invFft(polyWork0, m); // (1/2) E(n)
invFft(polyWork0 + m, m); // (1/2) E(n)
for (int i = 0; i < m >> 1; ++i) polyWork0[(m >> 1) + i] += polyWork0[m + i];
for (int i = 0; i < m; ++i) polyWork0[i] *= ::inv[m + i];
for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork0[i] += (*this)[m + i];
memset(polyWork0 + m, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 1 E(n)
memcpy(fs.data() + m, polyWork0, (n - m) * sizeof(Mint));
return fs;
}
// (29 + 1/2) E(n)
// g <- g - (log g - a log t) g
Poly pow1(Mint a, int n) const {
assert(!empty()); assert((*this)[0].x == 1U); assert(1 <= n);
return (a * log(n)).exp(n); // 13 E(n) + (16 + 1/2) E(n)
}
// (29 + 1/2) E(n - a ord(t))
Poly pow(long long a, int n) const {
assert(a >= 0); assert(1 <= n);
if (a == 0) { Poly gs(n); gs[0].x = 1U; return gs; }
const int o = ord();
if (o == -1 || o > (n - 1) / a) return Poly(n);
const Mint b = (*this)[o].inv(), c = (*this)[o].pow(a);
const int ntt = min<int>(n - a * o, size() - o);
Poly tts(ntt);
for (int i = 0; i < ntt; ++i) tts[i] = b * (*this)[o + i];
tts = tts.pow1(a, n - a * o); // (29 + 1/2) E(n - a ord(t))
Poly gs(n);
for (int i = 0; i < n - a * o; ++i) gs[a * o + i] = c * tts[i];
return gs;
}
// (10 + 1/2) E(n)
// f = t^(1/2) mod x^m, g = 1 / t^(1/2) mod x^m
// f <- f - (f^2 - h) g / 2
// g <- g - (f g - 1) g
// polyWork1: DFT(f, m), polyWork2: g, polyWork3: DFT(g, 2 m)
Poly sqrt(int n) const {
assert(!empty()); assert((*this)[0].x == 1U); assert(1 <= n);
assert(n == 1 || 1 << (32 - __builtin_clz(n - 1)) <= LIM_POLY);
if (n == 1) return {1U};
if (n == 2) return {1U, at(1) / 2};
Poly fs(n);
fs[0].x = polyWork1[0].x = polyWork2[0].x = 1U;
int m;
for (m = 1; m << 1 < n; m <<= 1) {
for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork1[i];
invFft(polyWork1, m); // (1/2) E(n)
for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork1[i] -= (*this)[i];
for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork1[i] -= (*this)[m + i];
memset(polyWork1 + m, 0, m * sizeof(Mint));
fft(polyWork1, m << 1); // 1 E(n)
memcpy(polyWork3, polyWork2, m * sizeof(Mint));
memset(polyWork3 + m, 0, m * sizeof(Mint));
fft(polyWork3, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork1[i] *= polyWork3[i];
invFft(polyWork1, m << 1); // 1 E(n)
for (int i = 0; i < m; ++i) { polyWork1[i] = -polyWork1[i]; fs[m + i].x = ((polyWork1[i].x & 1) ? (polyWork1[i].x + MO) : polyWork1[i].x) >> 1; }
memcpy(polyWork1, fs.data(), (m << 1) * sizeof(Mint));
fft(polyWork1, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] = polyWork1[i] * polyWork3[i];
invFft(polyWork0, m << 1); // 1 E(n)
memset(polyWork0, 0, m * sizeof(Mint));
fft(polyWork0, m << 1); // 1 E(n)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork3[i];
invFft(polyWork0, m << 1); // 1 E(n)
for (int i = m; i < m << 1; ++i) polyWork2[i] = -polyWork0[i];
}
for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork1[i];
invFft(polyWork1, m); // (1/2) E(n)
for (int i = 0, i0 = min(m, size()); i < i0; ++i) polyWork1[i] -= (*this)[i];
for (int i = 0, i0 = min(m, size() - m); i < i0; ++i) polyWork1[i] -= (*this)[m + i];
memcpy(polyWork1 + m, polyWork1 + (m >> 1), (m >> 1) * sizeof(Mint));
memset(polyWork1 + (m >> 1), 0, (m >> 1) * sizeof(Mint));
memset(polyWork1 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
fft(polyWork1, m); // (1/2) E(n)
fft(polyWork1 + m, m); // (1/2) E(n)
memcpy(polyWork3 + m, polyWork2 + (m >> 1), (m >> 1) * sizeof(Mint));
memset(polyWork3 + m + (m >> 1), 0, (m >> 1) * sizeof(Mint));
fft(polyWork3 + m, m); // (1/2) E(n)
// for (int i = 0; i < m << 1; ++i) polyWork1[i] *= polyWork3[i];
for (int i = 0; i < m; ++i) polyWork1[m + i] = polyWork1[i] * polyWork3[m + i] + polyWork1[m + i] * polyWork3[i];
for (int i = 0; i < m; ++i) polyWork1[i] *= polyWork3[i];
invFft(polyWork1, m); // (1/2) E(n)
invFft(polyWork1 + m, m); // (1/2) E(n)
for (int i = 0; i < m >> 1; ++i) polyWork1[(m >> 1) + i] += polyWork1[m + i];
for (int i = 0; i < n - m; ++i) { polyWork1[i] = -polyWork1[i]; fs[m + i].x = ((polyWork1[i].x & 1) ? (polyWork1[i].x + MO) : polyWork1[i].x) >> 1; }
return fs;
}
// (10 + 1/2) E(n)
// modSqrt must return a quadratic residue if exists, or anything otherwise.
// Return {} if *this does not have a square root.
template <class F> Poly sqrt(int n, F modSqrt) const {
assert(1 <= n);
const int o = ord();
if (o == -1) return Poly(n);
if (o & 1) return {};
const Mint c = modSqrt((*this)[o]);
if (c * c != (*this)[o]) return {};
if (o >> 1 >= n) return Poly(n);
const Mint b = (*this)[o].inv();
const int ntt = min(n - (o >> 1), size() - o);
Poly tts(ntt);
for (int i = 0; i < ntt; ++i) tts[i] = b * (*this)[o + i];
tts = tts.sqrt(n - (o >> 1)); // (10 + 1/2) E(n)
Poly gs(n);
for (int i = 0; i < n - (o >> 1); ++i) gs[(o >> 1) + i] = c * tts[i];
return gs;
}
// 6 E(|t|)
// x -> x + a
Poly shift(const Mint &a) const {
if (empty()) return {};
const int n = size();
int m = 1;
for (; m < n; m <<= 1) {}
for (int i = 0; i < n; ++i) polyWork0[i] = fac[i] * (*this)[i];
memset(polyWork0 + n, 0, ((m << 1) - n) * sizeof(Mint));
fft(polyWork0, m << 1); // 2 E(|t|)
{
Mint aa = 1;
for (int i = 0; i < n; ++i) { polyWork1[n - 1 - i] = invFac[i] * aa; aa *= a; }
}
memset(polyWork1 + n, 0, ((m << 1) - n) * sizeof(Mint));
fft(polyWork1, m << 1); // 2 E(|t|)
for (int i = 0; i < m << 1; ++i) polyWork0[i] *= polyWork1[i];
invFft(polyWork0, m << 1); // 2 E(|t|)
Poly fs(n);
for (int i = 0; i < n; ++i) fs[i] = invFac[i] * polyWork0[n - 1 + i];
return fs;
}
};
Mint linearRecurrenceAt(const vector<Mint> &as, const vector<Mint> &cs, long long k) {
assert(!cs.empty()); assert(cs[0]);
const int d = cs.size() - 1;
assert(as.size() >= static_cast<size_t>(d));
return (Poly(vector<Mint>(as.begin(), as.begin() + d)) * cs).mod(d).divAt(cs, k);
}
struct SubproductTree {
int logN, n, nn;
vector<Mint> xs;
// [DFT_4((X-xs[0])(X-xs[1])(X-xs[2])(X-xs[3]))] [(X-xs[0])(X-xs[1])(X-xs[2])(X-xs[3])mod X^4]
// [ DFT_4((X-xs[0])(X-xs[1])) ] [ DFT_4((X-xs[2])(X-xs[3])) ]
// [ DFT_2(X-xs[0]) ] [ DFT_2(X-xs[1]) ] [ DFT_2(X-xs[2]) ] [ DFT_2(X-xs[3]) ]
vector<Mint> buf;
vector<Mint *> gss;
// (1 - xs[0] X) ... (1 - xs[nn-1] X)
Poly all;
// (ceil(log_2 n) + O(1)) E(n)
SubproductTree(const vector<Mint> &xs_) {
n = xs_.size();
for (logN = 0, nn = 1; nn < n; ++logN, nn <<= 1) {}
xs.assign(nn, 0U);
memcpy(xs.data(), xs_.data(), n * sizeof(Mint));
buf.assign((logN + 1) * (nn << 1), 0U);
gss.assign(nn << 1, nullptr);
for (int h = 0; h <= logN; ++h) for (int u = 1 << h; u < 1 << (h + 1); ++u) {
gss[u] = buf.data() + (h * (nn << 1) + ((u - (1 << h)) << (logN - h + 1)));
}
for (int i = 0; i < nn; ++i) {
gss[nn + i][0] = -xs[i] + 1;
gss[nn + i][1] = -xs[i] - 1;
}
if (nn == 1) gss[1][1] += 2;
for (int h = logN; --h >= 0; ) {
const int m = 1 << (logN - h);
for (int u = 1 << (h + 1); --u >= 1 << h; ) {
for (int i = 0; i < m; ++i) gss[u][i] = gss[u << 1][i] * gss[u << 1 | 1][i];
memcpy(gss[u] + m, gss[u], m * sizeof(Mint));
invFft(gss[u] + m, m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
if (h > 0) {
gss[u][m] -= 2;
const Mint a = FFT_ROOTS[logN - h + 1];
Mint aa = 1;
for (int i = m; i < m << 1; ++i) { gss[u][i] *= aa; aa *= a; };
fft(gss[u] + m, m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
}
}
}
all.resize(nn + 1);
all[0] = 1;
for (int i = 1; i < nn; ++i) all[i] = gss[1][nn + nn - i];
all[nn] = gss[1][nn] - 1;
}
// ((3/2) ceil(log_2 n) + O(1)) E(n) + 10 E(|f|) + 3 E(|f| + 2^(ceil(log_2 n)))
vector<Mint> multiEval(const Poly &fs) const {
vector<Mint> work0(nn), work1(nn), work2(nn);
{
const int m = max(fs.size(), 1);
auto invAll = all.inv(m); // 10 E(|f|)
std::reverse(invAll.begin(), invAll.end());
int mm;
for (mm = 1; mm < m - 1 + nn; mm <<= 1) {}
invAll.resize(mm, 0U);
fft(invAll); // E(|f| + 2^(ceil(log_2 n)))
vector<Mint> ffs(mm, 0U);
if (fs.size()) memcpy(ffs.data(), fs.data(), fs.size() * sizeof(Mint));
fft(ffs); // E(|f| + 2^(ceil(log_2 n)))
for (int i = 0; i < mm; ++i) ffs[i] *= invAll[i];
invFft(ffs); // E(|f| + 2^(ceil(log_2 n)))
memcpy(((logN & 1) ? work1 : work0).data(), ffs.data() + m - 1, nn * sizeof(Mint));
}
for (int h = 0; h < logN; ++h) {
const int m = 1 << (logN - h);
for (int u = 1 << h; u < 1 << (h + 1); ++u) {
Mint *hs = (((logN - h) & 1) ? work1 : work0).data() + ((u - (1 << h)) << (logN - h));
Mint *hs0 = (((logN - h) & 1) ? work0 : work1).data() + ((u - (1 << h)) << (logN - h));
Mint *hs1 = hs0 + (m >> 1);
fft(hs, m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
for (int i = 0; i < m; ++i) work2[i] = gss[u << 1 | 1][i] * hs[i];
invFft(work2.data(), m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
memcpy(hs0, work2.data() + (m >> 1), (m >> 1) * sizeof(Mint));
for (int i = 0; i < m; ++i) work2[i] = gss[u << 1][i] * hs[i];
invFft(work2.data(), m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
memcpy(hs1, work2.data() + (m >> 1), (m >> 1) * sizeof(Mint));
}
}
work0.resize(n);
return work0;
}
// ((5/2) ceil(log_2 n) + O(1)) E(n)
Poly interpolate(const vector<Mint> &ys) const {
assert(static_cast<int>(ys.size()) == n);
Poly gs(n);
for (int i = 0; i < n; ++i) gs[i] = (i + 1) * all[n - (i + 1)];
const vector<Mint> denoms = multiEval(gs); // ((3/2) ceil(log_2 n) + O(1)) E(n)
vector<Mint> work(nn << 1, 0U);
for (int i = 0; i < n; ++i) {
// xs[0], ..., xs[n - 1] are not distinct
assert(denoms[i]);
work[i << 1] = work[i << 1 | 1] = ys[i] / denoms[i];
}
for (int h = logN; --h >= 0; ) {
const int m = 1 << (logN - h);
for (int u = 1 << (h + 1); --u >= 1 << h; ) {
Mint *hs = work.data() + ((u - (1 << h)) << (logN - h + 1));
for (int i = 0; i < m; ++i) hs[i] = gss[u << 1 | 1][i] * hs[i] + gss[u << 1][i] * hs[m + i];
if (h > 0) {
memcpy(hs + m, hs, m * sizeof(Mint));
invFft(hs + m, m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
const Mint a = FFT_ROOTS[logN - h + 1];
Mint aa = 1;
for (int i = m; i < m << 1; ++i) { hs[i] *= aa; aa *= a; };
fft(hs + m, m); // ((1/2) ceil(log_2 n) + O(1)) E(n)
}
}
}
invFft(work.data(), nn); // E(n)
return Poly(vector<Mint>(work.data() + nn - n, work.data() + nn));
}
};
// q: rev([0, m]) * [0, n], [t^0] q(t, x) = 1 omitted
// ret: [0, m-1] * [0, n]
vector<Mint> comRec(int m, int n, const vector<Mint> &as, const vector<Mint> &qss) {
if (!n) { auto ret = as; ret.resize(m, 0); return ret; }
// reuse DFT(q(t, -x)); (2n+2) instead of (2n+1)
int len;
for (len = 2; len < (2*m) * (2*n+2); len <<= 1) {}
vector<Mint> qs(len, 0);
for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) qs[i * (2*n+2) + j] = qss[i * (n+1) + j];
fft(qs);
vector<Mint> work(len >> 1, 0);
for (int k = 0; k < len >> 1; ++k) { work[k] = qs[k << 1] * qs[k << 1 | 1]; swap(qs[k << 1], qs[k << 1 | 1]); }
invFft(work);
vector<Mint> qqss((2*m) * (n/2+1), 0);
for (int i = 0; i < 2*m-1; ++i) for (int j = 0; j <= n/2; ++j) qqss[i * (n/2+1) + j] = work[i * (n+1) + j];
for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) qqss[(m+i) * (n/2+1) + j] += qss[i * (n+1) + 2*j] + qss[i * (n+1) + 2*j];
const auto res = comRec(2*m, n/2, as, qqss);
vector<Mint> ps(len, 0);
for (int i = 0; i < 2*m; ++i) for (int j = 0; j <= n/2; ++j) ps[i * (2*n+2) + (2*n+1) - (2*j+(n&1))] = res[i * (n/2+1) + j];
fft(ps);
for (int k = 0; k < len; ++k) ps[k] *= qs[k];
invFft(ps);
vector<Mint> ret(m * (n+1));
for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) ret[i * (n+1) + j] = ps[(m+i) * (2*n+2) + (2*n+1) - j];
for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) ret[i * (n+1) + (2*j+(n&1))] += res[i * (n/2+1) + j];
return ret;
}
// a(b(x)) mod x^n
// transpose and rev: p(x) -> [x^(n-1)] p(x) b(x)^i for each 0 <= i < n
// [x^(n-1)] p(x) / (1 - t b(x))
vector<Mint> com(const vector<Mint> &as, const vector<Mint> &bs, int n) {
assert(bs.size() == 0 || !bs[0]);
if (n == 0) return {};
vector<Mint> qss(n, 0);
for (int j = 0; j < min<int>(bs.size(), n); ++j) qss[j] = -bs[j];
auto cs = comRec(1, n - 1, as, qss);
reverse(cs.begin(), cs.end());
return cs;
}
// [x^(n-1)] a(x) b(x)^i for each 0 <= i < n
// [x^(n-1)] a(x) / (1 - t b(x))
vector<Mint> powProj(const vector<Mint> &as, const vector<Mint> &bs, int n) {
assert(bs.size() == 0 || !bs[0]);
assert(n >= 1);
// p(t, x): [0, m-1] * [0, n]
// q(t, x): rev([0, m]) * [0, n], [t^0] q(t, x) = 1 omitted
vector<Mint> pss(n, 0), qss(n, 0);
for (int j = 0; j < min<int>(as.size(), n); ++j) pss[j] = as[j];
for (int j = 0; j < min<int>(bs.size(), n); ++j) qss[j] = -bs[j];
const int n0 = n--;
for (int m = 1; n; m *= 2, n /= 2) {
// reuse DFT(q(t, -x)); (2n+2) instead of (2n+1)
int len;
for (len = 2; len < (m*2) * (n*2+2); len <<= 1) {}
vector<Mint> qs(len, 0);
for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) qs[i * (n*2+2) + j] = qss[i * (n+1) + j];
fft(qs);
vector<Mint> work(len >> 1, 0);
for (int k = 0; k < len >> 1; ++k) { work[k] = qs[k << 1] * qs[k << 1 | 1]; swap(qs[k << 1], qs[k << 1 | 1]); }
invFft(work);
vector<Mint> qqss((m*2) * (n/2+1), 0);
for (int i = 0; i <= m*2-2; ++i) for (int j = 0; j <= n/2; ++j) qqss[i * (n/2+1) + j] = work[i * (n+1) + j];
for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) qqss[(m+i) * (n/2+1) + j] += qss[i * (n+1) + j*2] + qss[i * (n+1) + j*2];
vector<Mint> ps(len, 0);
for (int i = 0; i < m; ++i) for (int j = 0; j <= n; ++j) ps[i * (n*2+2) + j] = pss[i * (n+1) + j];
fft(ps);
for (int k = 0; k < len; ++k) ps[k] *= qs[k];
invFft(ps);
vector<Mint> ppss((m*2) * (n/2+1), 0);
for (int i = 0; i <= m*2-2; ++i) for (int j = 0; j <= n/2; ++j) ppss[i * (n/2+1) + j] = ps[i * (n*2+2) + (j*2+(n&1))];
for (int i = 0; i < m; ++i) for (int j = 0; j <= n/2; ++j) ppss[(m+i) * (n/2+1) + j] += pss[i * (n+1) + (j*2+(n&1))];
pss.swap(ppss);
qss.swap(qqss);
}
std::reverse(pss.begin(), pss.end());
pss.resize(n0, 0);
return pss;
}
// a^<-1>(x) mod x^n
// (n-1) [x^(n-1)] a(x)^i = i [x^(n-1-i)] (x/a^<-1>(x))^(n-1)
Poly comInv(const Poly &as, int n) {
assert(as.size() >= 2); assert(!as[0]); assert(as[1]);
assert(n >= 0);
if (n <= 1) return Poly(n);
// reduce to [x^1] a(x) = 1 in order to take (n-1)-th root
const Mint t = as[1].inv();
const auto res = powProj({1}, t * as, n);
Poly ret(n - 1);
for (int i = 1; i < n; ++i) ret[n - 1 - i] = inv[i] * (n - 1) * res[i];
ret = ret.pow1(-inv[n - 1], n - 1);
ret.insert(ret.begin(), 0);
{ Mint tt = 1; for (int i = 0; i < n; ++i) { ret[i] *= tt; tt *= t; } }
return ret;
}
////////////////////////////////////////////////////////////////////////////////
int SUB, N, M;
vector<Mint> A, B, C, D, E, F;
// \sum[0<=i<=N] E[i] (1+B[0]x)...(1+B[i-1]x) x^(N-i)
pair<Poly, Poly> dfsI(int l, int r) {
if (l + 1 == r) {
return make_pair(Poly{1, B[l]}, Poly{E[l]});
} else {
const int m = (l + r) / 2;
const auto resL = dfsI(l, m);
const auto resR = dfsI(m, r);
Poly fs = resL.second;
fs.insert(fs.begin(), r - m, 0);
return make_pair(resL.first * resR.first, fs + resL.first * resR.second);
}
}
// \sum[0<=j<=M] F[j] (A[0]x)...(A[j-1]x) (1/(1-C[0]x))...(1/(1-C[j]x))
vector<Mint> AProd;
pair<Poly, Poly> dfsJ(int l, int r) {
if (l + 1 == r) {
return make_pair(Poly{1, -C[l]}, Poly{F[l] * AProd[l]});
} else {
const int m = (l + r) / 2;
const auto resL = dfsJ(l, m);
const auto resR = dfsJ(m, r);
Poly fs = resR.second;
fs.insert(fs.begin(), m - l, 0);
return make_pair(resL.first * resR.first, resL.second * resR.first + fs);
}
}
int main() {
for (; ~scanf("%d%d%d%*d", &SUB, &N, &M); ) {
A.assign(M + 1, 0); for (int j = 0; j < M; ++j) scanf("%u", &A[j].x);
B.assign(N + 1, 0); for (int i = 0; i < N; ++i) scanf("%u", &B[i].x);
C.assign(M + 1, 0); for (int j = 0; j <= M; ++j) scanf("%u", &C[j].x);
D.assign(M + 1, 0); for (int j = 1; j <= M; ++j) scanf("%u", &D[j].x);
E.assign(N + 1, 0); for (int i = 0; i <= N; ++i) scanf("%u", &E[i].x);
F.assign(M + 1, 0); for (int j = 0; j <= M; ++j) scanf("%u", &F[j].x);
const auto resI = dfsI(0, N + 1);
const Poly I = resI.second;
AProd.resize(M + 1);
AProd[0] = 1;
for (int j = 0; j < M; ++j) AProd[j + 1] = AProd[j] * A[j];
const auto resJ = dfsJ(0, M + 1);
const Poly J = resJ.second.div(resJ.first, N + 1);
Mint ans = 0;
for (int k = 0; k <= N; ++k) ans += I[k] * J[N - k];
printf("%u\n", ans.x);
}
return 0;
}
Details
Tip: Click on the bar to expand more detailed information
Subtask #1:
score: 0
Wrong Answer
Test #1:
score: 0
Wrong Answer
time: 27ms
memory: 21568kb
input:
1 5000 5000 998244353 121811167 379924090 361631583 174189813 559424693 889647308 193102812 469875055 32237660 96186933 624360154 404866088 859165067 748410791 926700198 368632735 476560636 798138824 17883437 712872225 448819400 33122704 572152288 627010263 336722521 775918346 465722630 681224329 60...
output:
893602296
result:
wrong answer 1st lines differ - expected: '698779876', found: '893602296'
Subtask #2:
score: 5
Accepted
Test #4:
score: 5
Accepted
time: 755ms
memory: 38420kb
input:
2 200000 200000 998244353 903563506 433239074 632684755 970178226 831892753 932120646 157832416 517140217 296101978 998244343 850946564 2367119 708278025 376919151 752106478 994362560 806760771 672325565 9 958330492 343658496 153627310 330649130 983441587 829707074 135609242 706388812 325297147 4972...
output:
108905794
result:
ok single line: '108905794'
Test #5:
score: 5
Accepted
time: 751ms
memory: 38840kb
input:
2 199910 194100 998244353 587911377 573048398 832688590 809066619 524442920 218487661 649170169 8 150333233 204150153 800582862 464558080 291668841 361834956 998244344 998244349 806341682 775965963 459031329 867640103 425129750 7 998244343 274941091 809744915 443910210 859200100 998244350 725497297 ...
output:
597176160
result:
ok single line: '597176160'
Test #6:
score: 5
Accepted
time: 735ms
memory: 35012kb
input:
2 150810 200000 998244353 288330007 105173193 991831123 698131025 301828280 273289882 387551340 542768677 115972971 425381688 811911805 962095963 566257196 435928108 337873530 109252306 933737641 967531573 29209951 787608009 497111219 315932660 878605444 903737754 260092904 447237039 37123388 594371...
output:
78643368
result:
ok single line: '78643368'
Subtask #3:
score: 8
Accepted
Test #7:
score: 8
Accepted
time: 773ms
memory: 36176kb
input:
3 200000 200000 998244353 493605813 622283646 579332647 528537957 211102509 336893985 106292723 166710741 443808575 180234697 744331593 252016663 693453924 975202110 519712748 174749950 250990381 584528219 619047448 641168296 914918741 785005029 433175528 400547992 845115512 278159630 227330862 6407...
output:
144815343
result:
ok single line: '144815343'
Test #8:
score: 8
Accepted
time: 759ms
memory: 35372kb
input:
3 200000 199996 998244353 742 699221406 301 364485093 804 294 873 282584633 204 882 889 438104412 349 737559671 152908512 490 206 11 613 175 155 659898955 983 206756334 273919067 843 898 359154783 360 612 201572138 207697004 396 482 136 361 339 831490112 279585283 516749439 318 198386232 69807322 72...
output:
758876599
result:
ok single line: '758876599'
Test #9:
score: 8
Accepted
time: 364ms
memory: 32456kb
input:
3 1 200000 998244353 734 796 543458399 920002767 990425370 338 107321680 90089647 557777271 434728777 503 709 616 850 448 884 143526194 299065437 45361809 666577969 923 492965257 850 229 769 183867189 675 463141085 791 281 718 763 240353338 711 671313626 203252664 888743426 450 223 446 542052402 565...
output:
53421494
result:
ok single line: '53421494'
Test #10:
score: 8
Accepted
time: 419ms
memory: 29988kb
input:
3 200000 1 998244353 277786488 264242127 2 319061241 5 8 2 1 51038317 673171629 267552091 920379441 5 880796442 618282833 514210179 9 8 67744028 9 1 691195881 92394487 719900388 638641507 5 663287704 4 579997250 985749664 984980853 4 519857232 4 735668111 3 7 360023171 7 7 234677046 359140782 8 5227...
output:
140453572
result:
ok single line: '140453572'
Subtask #4:
score: 8
Accepted
Dependency #2:
100%
Accepted
Test #11:
score: 8
Accepted
time: 754ms
memory: 35292kb
input:
4 200000 200000 998244353 199752876 397435299 166254444 779278347 182416898 901682085 69090988 738092954 712758505 681482091 950425501 897402866 870967623 36849705 515322829 172446310 900429861 34423824 785446470 553333476 938968413 86270822 973639556 431435054 741273084 677434142 222677171 61928101...
output:
65201275
result:
ok single line: '65201275'
Test #12:
score: 8
Accepted
time: 410ms
memory: 29256kb
input:
4 199996 2000 998244353 865 499 31904732 259114092 967550739 255 785 143545371 699 573 216 789 553653892 329548961 451423632 886 741679575 78720133 464611720 972 242 869 294933191 351 788310280 697 861990346 731 734 701 407 401383418 835 189335197 668 758851074 250436747 962105421 562203389 11011899...
output:
85153184
result:
ok single line: '85153184'
Test #13:
score: 8
Accepted
time: 448ms
memory: 30472kb
input:
4 40000 190000 998244353 953 586890334 906 161 524210942 327745937 100 231 970 63891350 28 370040454 628796638 275 239621470 708518643 995 351 376514813 814869235 66 742 561 483786281 5379525 603 608 274128226 256 844 662 445297819 123782163 920 957935561 621339377 26129412 763585646 451 911 5350051...
output:
99807982
result:
ok single line: '99807982'
Subtask #5:
score: 5
Accepted
Dependency #3:
100%
Accepted
Dependency #4:
100%
Accepted
Test #14:
score: 5
Accepted
time: 767ms
memory: 35624kb
input:
5 200000 200000 998244353 357 270 635 81 18806969 807 433 242 242 154 347094190 96 953 893447557 497 265754515 220901127 205 958 488 212420388 457 799989458 726134433 547661079 946786926 891 955400819 664511484 701 155 791 4 845 954462573 536216065 136 54 251 481 584310172 389122260 34 527707943 300...
output:
366841924
result:
ok single line: '366841924'
Test #15:
score: 5
Accepted
time: 762ms
memory: 38712kb
input:
5 199999 190000 998244353 930048626 861740403 534792638 856 178 908419417 742 685 40557520 132 779489457 792 468990716 153965694 395450527 688 398541680 651 789664707 703746126 579 588989999 936017101 652728925 314526790 297 962 771009949 626 192587923 377598397 654 459 187173523 857466261 567 73215...
output:
755407887
result:
ok single line: '755407887'
Subtask #6:
score: 0
Wrong Answer
Test #16:
score: 15
Accepted
time: 757ms
memory: 38744kb
input:
6 200000 200000 998244353 401806059 107033001 530043262 862506025 263940497 48524969 232075248 849682830 420464058 64900333 394986647 954304290 957385577 86269798 579307969 896967626 230611640 527078096 39773429 402432856 495204529 272090833 100466767 562115973 196636941 736050044 580541546 81233872...
output:
562220526
result:
ok single line: '562220526'
Test #17:
score: 0
Wrong Answer
time: 751ms
memory: 38136kb
input:
6 190000 180000 998244353 331419431 585273774 326021268 911984877 504951700 663 667 967180428 535316997 888573230 112066317 286249963 593 994 230 367194808 906865758 946973557 921 1 302880572 377068830 444057418 953 1000 946 906 274661283 90 345429012 719 547222223 806 856 911 970 110923264 26390908...
output:
215419983
result:
wrong answer 1st lines differ - expected: '467130256', found: '215419983'
Subtask #7:
score: 0
Skipped
Dependency #1:
0%
Subtask #8:
score: 0
Skipped
Dependency #6:
0%
Subtask #9:
score: 0
Skipped
Dependency #5:
100%
Accepted
Dependency #7:
0%
Subtask #10:
score: 0
Skipped
Dependency #1:
0%