QOJ.ac
QOJ
ID | 题目 | 提交者 | 结果 | 用时 | 内存 | 语言 | 文件大小 | 提交时间 | 测评时间 |
---|---|---|---|---|---|---|---|---|---|
#101842 | #5448. 另一个欧拉数问题 | hos_lyric# | 80 | 950ms | 34348kb | C++14 | 39.6kb | 2023-05-01 14:00:14 | 2024-05-26 02:49:30 |
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 <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; }
////////////////////////////////////////////////////////////////////////////////
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;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// 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)
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];
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 pow(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.pow(Mint(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);
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));
}
};
////////////////////////////////////////////////////////////////////////////////
Mint binom(Int n, Int k) {
if (n < 0) {
if (k >= 0) {
return ((k & 1) ? -1 : +1) * binom(-n + k - 1, k);
} else if (n - k >= 0) {
return (((n - k) & 1) ? -1 : +1) * binom(-k - 1, n - k);
} else {
return 0;
}
} else {
if (0 <= k && k <= n) {
assert(n < LIM_INV);
return fac[n] * invFac[k] * invFac[n - k];
} else {
return 0;
}
}
}
int A, N, L, M;
vector<int> P;
int K;
namespace brute {
Mint f[2010][2010];
Mint run() {
memset(f, 0, sizeof(f));
f[M][K] = 1;
for (int i = M; i < N; ++i) {
for (int j = 0; j <= i; ++j) {
f[i + 1][j] += f[i][j] * j;
f[i + 1][j + 1] += f[i][j] * (A * i + 1 - j);
}
}
// for(int i=M;i<=N;++i){cerr<<i<<": ";pv(f[i],f[i]+i+1);}
return f[N][L];
}
} // brute
/*
f[i][j]: (j-1) descents
f[i][j] = (A(i-1)+1 - (j-1)) f[i-1][j-1] + j f[i-1][j]
f[i] = ((A(i-1)+1)x - x(1-x) (d/dx)) f[i-1]
g[i] := (1-x)^(Ai+1) f[i]
g[i] = x (1-x)^-A (d/dx) g[i-1]
*/
namespace a1 {
Mint run() {
// (1-x)^-(M+1) x^K
vector<Mint> gs(L + 1, 0);
for (int j = K; j <= L; ++j) {
gs[j] = (((j-K)&1)?-1:+1) * binom(-(M+1), j - K);
}
for (int j = K; j <= L; ++j) {
gs[j] *= Mint(j).pow(N-M);
}
// *= (1-x)^(N+1)
Mint ans = 0;
for (int j = 0; j <= L; ++j) {
ans += ((j&1)?-1:+1) * binom(N+1, j) * gs[L - j];
}
return ans;
}
} // a1
/*
multiset of *=k, ..., *=l
\sum[0<=k<=l<=L] b[k] ([y^(N-M)] (1-ky)^-1 ... (1-ly)^-1) c[l]
*/
namespace a2m1 {
vector<Mint> bs, cs;
struct Info {
Poly h, hB, hC, hBC;
};
Info rec(int l, int r) {
Info ret;
if (l + 1 == r) {
ret.h = {1, -l};
ret.hB = {bs[l]};
ret.hC = {cs[l]};
ret.hBC = {bs[l] * cs[l]};
} else {
const int mid = (l + r) / 2;
const auto resL = rec(l, mid);
const auto resR = rec(mid, r);
ret.h = resL.h * resR.h;
ret.hB = resL.hB + resL.h * resR.hB;
ret.hC = resL.hC * resR.h + resR.hC;
ret.hBC = resL.hBC * resR.h + resL.hB * resR.hC + resL.h * resR.hBC;
}
// cerr<<"rec "<<l<<" "<<r<<" = "<<ret.h<<" "<<ret.hB<<" "<<ret.hC<<" "<<ret.hBC<<endl;
return ret;
}
Mint run() {
// (1-x)^-(2M+1) x^K
// undo (1-x)^-1
bs.resize(L + 1);
for (int j = 0; j <= L - K; ++j) {
bs[K + j] = ((j&1)?-1:+1) * binom(-(2*M), j);
}
// *= (1-x)^(2N+1)
cs.resize(L + 1);
for (int j = 0; j <= L; ++j) {
cs[L - j] = ((j&1)?-1:+1) * binom(2*N+1, j);
}
// cerr<<"bs = "<<bs<<endl;
// cerr<<"cs = "<<cs<<endl;
const auto res = rec(0, L + 1);
const Mint ans = res.hBC.divAt(res.h, N-M);
return ans;
}
} // a2m1
int main() {
for (; ~scanf("%d%d%d%d", &A, &N, &L, &M); ) {
P.resize(A * M);
for (int i = 0; i < A * M; ++i) {
scanf("%d", &P[i]);
}
++L;
K = 1;
for (int i = 0; i < A * M - 1; ++i) if (P[i] > P[i + 1]) {
++K;
}
cerr<<"A = "<<A<<"; ("<<M<<", "<<K<<") -> ("<<N<<", "<<L<<")"<<endl;
Mint ans = 0;
if (A == 1) {
ans = a1::run();
} else if (A == 2) {
ans = a2m1::run();
} else {
ans = brute::run();
}
printf("%u\n", ans.x);
#ifdef LOCAL
if(N<=2000){
const Mint brt=brute::run();
cerr<<"brt = "<<brt<<endl;
assert(brt==ans);
}
#endif
}
return 0;
}
詳細信息
Subtask #1:
score: 10
Accepted
Test #1:
score: 10
Accepted
time: 27ms
memory: 33640kb
input:
584 1985 1017 186 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
output:
683567838
result:
ok 1 number(s): "683567838"
Test #2:
score: 0
Accepted
time: 21ms
memory: 33228kb
input:
549 1352 215 144 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1...
output:
514220824
result:
ok 1 number(s): "514220824"
Test #3:
score: 0
Accepted
time: 23ms
memory: 32112kb
input:
280 1833 1153 253 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 ...
output:
801426026
result:
ok 1 number(s): "801426026"
Test #4:
score: 0
Accepted
time: 21ms
memory: 33380kb
input:
886 1885 1736 131 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ...
output:
732050911
result:
ok 1 number(s): "732050911"
Test #5:
score: 0
Accepted
time: 28ms
memory: 32928kb
input:
464 1859 1766 385 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 55 150 150 150 150 150 150 150 150 150 150 150 150 150 150 150...
output:
268474724
result:
ok 1 number(s): "268474724"
Test #6:
score: 0
Accepted
time: 23ms
memory: 33436kb
input:
468 1281 718 285 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1...
output:
7958328
result:
ok 1 number(s): "7958328"
Test #7:
score: 0
Accepted
time: 25ms
memory: 33132kb
input:
973 1966 417 99 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ...
output:
943626008
result:
ok 1 number(s): "943626008"
Test #8:
score: 0
Accepted
time: 15ms
memory: 33700kb
input:
224 1781 1724 358 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ...
output:
550100613
result:
ok 1 number(s): "550100613"
Test #9:
score: 0
Accepted
time: 21ms
memory: 33552kb
input:
952 1793 939 19 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 1...
output:
357850463
result:
ok 1 number(s): "357850463"
Test #10:
score: 0
Accepted
time: 25ms
memory: 33352kb
input:
978 1415 776 108 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1...
output:
284767772
result:
ok 1 number(s): "284767772"
Subtask #2:
score: 10
Accepted
Test #11:
score: 10
Accepted
time: 21ms
memory: 18432kb
input:
1 199913 100055 1 1
output:
856065368
result:
ok 1 number(s): "856065368"
Test #12:
score: 0
Accepted
time: 20ms
memory: 18156kb
input:
1 197151 139288 1 1
output:
933745116
result:
ok 1 number(s): "933745116"
Test #13:
score: 0
Accepted
time: 19ms
memory: 18224kb
input:
1 198165 124475 1 1
output:
650379731
result:
ok 1 number(s): "650379731"
Test #14:
score: 0
Accepted
time: 23ms
memory: 18360kb
input:
1 199407 183908 1 1
output:
213038974
result:
ok 1 number(s): "213038974"
Test #15:
score: 0
Accepted
time: 19ms
memory: 18152kb
input:
1 195940 66608 1 1
output:
203000250
result:
ok 1 number(s): "203000250"
Subtask #3:
score: 30
Accepted
Dependency #2:
100%
Accepted
Test #16:
score: 30
Accepted
time: 13ms
memory: 18136kb
input:
1 199968 100029 81 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
output:
882662058
result:
ok 1 number(s): "882662058"
Test #17:
score: 0
Accepted
time: 25ms
memory: 18360kb
input:
1 199237 160399 740 534 70 635 736 56 408 68 506 354 47 101 96 66 699 497 89 642 133 179 211 318 261 205 321 300 363 165 643 236 21 176 212 232 702 195 687 14 382 207 216 149 374 222 438 50 256 427 79 38 147 215 605 678 46 645 610 538 118 22 11 548 194 628 590 375 26 180 591 264 366 226 185 188 406 ...
output:
119843424
result:
ok 1 number(s): "119843424"
Test #18:
score: 0
Accepted
time: 28ms
memory: 18332kb
input:
1 197117 182307 1145 649 641 894 352 882 609 560 160 360 183 542 220 1059 750 489 977 810 964 48 812 540 346 210 345 284 969 371 274 978 1083 282 986 632 484 889 768 533 226 922 58 870 888 63 850 321 548 471 212 643 233 1109 994 1090 1145 979 405 418 555 394 781 1000 709 819 652 1117 703 397 715 102...
output:
321584073
result:
ok 1 number(s): "321584073"
Test #19:
score: 0
Accepted
time: 20ms
memory: 18596kb
input:
1 198080 196910 32807 13007 28354 16100 16920 19352 19610 14746 9160 6498 3691 32453 7386 16081 16509 7087 30770 27569 32126 31416 15799 12675 25754 20433 93 21896 24361 20050 15789 1724 5623 25658 31528 30773 30698 21920 5237 2778 29667 25560 237 8327 892 12072 9109 30004 22873 24074 29546 25501 26...
output:
0
result:
ok 1 number(s): "0"
Test #20:
score: 0
Accepted
time: 19ms
memory: 18140kb
input:
1 198858 125784 938 257 251 742 779 56 449 69 7 370 454 811 343 831 233 532 626 253 216 398 675 851 677 690 705 34 414 28 309 685 102 786 835 67 351 612 357 399 530 859 130 732 540 275 467 469 311 642 664 574 802 581 20 42 916 52 375 124 154 572 342 730 231 134 692 561 431 452 30 747 302 723 367 386...
output:
568776181
result:
ok 1 number(s): "568776181"
Subtask #4:
score: 15
Accepted
Test #21:
score: 15
Accepted
time: 466ms
memory: 30472kb
input:
2 199984 99989 1 1 1
output:
169573504
result:
ok 1 number(s): "169573504"
Test #22:
score: 0
Accepted
time: 950ms
memory: 32696kb
input:
2 198253 176846 1 1 1
output:
206376024
result:
ok 1 number(s): "206376024"
Test #23:
score: 0
Accepted
time: 496ms
memory: 29720kb
input:
2 197635 128767 1 1 1
output:
517345888
result:
ok 1 number(s): "517345888"
Test #24:
score: 0
Accepted
time: 892ms
memory: 31884kb
input:
2 198367 139324 1 1 1
output:
421558104
result:
ok 1 number(s): "421558104"
Test #25:
score: 0
Accepted
time: 949ms
memory: 34348kb
input:
2 197383 177188 1 1 1
output:
47912879
result:
ok 1 number(s): "47912879"
Subtask #5:
score: 15
Accepted
Dependency #4:
100%
Accepted
Test #26:
score: 15
Accepted
time: 470ms
memory: 29000kb
input:
2 199972 99897 63 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 ...
output:
646686125
result:
ok 1 number(s): "646686125"
Test #27:
score: 0
Accepted
time: 918ms
memory: 33852kb
input:
2 198230 150849 8875 2376 5489 5489 2376 3711 3711 5679 5679 4601 5973 5973 4601 554 4256 4596 4596 4256 554 4800 6366 6366 4800 2879 5348 5348 2879 1693 4230 6262 6993 8595 8595 6993 6262 4230 1693 2663 3954 3954 2663 14 14 7910 7910 4020 4020 3246 6907 6907 3246 8873 8873 4762 4762 520 520 1366 54...
output:
321924755
result:
ok 1 number(s): "321924755"
Test #28:
score: 0
Accepted
time: 913ms
memory: 32980kb
input:
2 199806 150052 89583 49845 49845 25156 25156 5539 20884 20884 5539 58128 84993 84993 58128 2142 47997 47997 65672 65672 22214 22214 10670 39754 74544 74544 48064 48064 39754 10670 7926 74018 74018 7926 33809 33809 14019 14019 38860 38860 2142 3555 71750 71750 26496 32708 32708 56732 56732 82049 820...
output:
160305146
result:
ok 1 number(s): "160305146"
Test #29:
score: 0
Accepted
time: 496ms
memory: 29016kb
input:
2 199723 124579 36954 26124 26124 18080 18080 7015 13358 13358 7015 300 15945 18773 18773 30036 30036 33473 33473 15945 2386 11858 11858 12166 12166 13792 13792 14189 14189 29559 29559 4323 4323 7264 23756 28872 28872 31656 31656 30918 30918 28504 28504 23756 7821 7821 10550 10550 14561 14561 15913 ...
output:
43656653
result:
ok 1 number(s): "43656653"
Test #30:
score: 0
Accepted
time: 481ms
memory: 33728kb
input:
2 199535 116407 31243 37 18798 22113 22113 18798 37 5174 19335 28778 28778 20844 20844 19335 5174 5511 5511 3022 14132 14132 6436 6436 3022 3531 4883 4883 5347 30026 30026 19094 19094 7086 7086 7300 15455 18190 18190 15455 7300 8333 14085 14085 8333 5347 5295 5295 15994 15994 27205 27205 3531 24459 ...
output:
966854518
result:
ok 1 number(s): "966854518"
Subtask #6:
score: 0
Runtime Error
Dependency #1:
100%
Accepted
Dependency #2:
100%
Accepted
Dependency #3:
100%
Accepted
Dependency #4:
100%
Accepted
Dependency #5:
100%
Accepted
Test #31:
score: 0
Runtime Error
input:
993 199975 99931 71 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...