QOJ.ac
QOJ
ID | 题目 | 提交者 | 结果 | 用时 | 内存 | 语言 | 文件大小 | 提交时间 | 测评时间 |
---|---|---|---|---|---|---|---|---|---|
#889145 | #9042. Fast Bogosort | HoMaMaOvO (Riku Kawasaki, Masaki Nishimoto, Yui Hosaka)# | AC ✓ | 149ms | 27976kb | C++14 | 43.9kb | 2025-02-08 14:41:38 | 2025-02-08 14:41:42 |
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 N;
vector<int> A;
int main() {
for (; ~scanf("%d", &N); ) {
A.resize(N);
for (int i = 0; i < N; ++i) {
scanf("%d", &A[i]);
--A[i];
}
Poly P(N + 1);
for (int i = 1; i <= N; ++i) P[i] = fac[i];
const Poly C = P.div(Poly{1} + P, N + 1);
// cerr<<"C = "<<C<<endl;
// 1 / (1 - C(1+E)) = F+Gt
/*
Poly E(N + 1), F(N + 1), G(N + 1);
F[0] = 1;
for (int i = 1; i <= N; ++i) {
for (int j = 1; j <= i; ++j) F[i] += F[i - j] * C[j];
for (int j = 1; j <= i; ++j) G[i] += F[i - j] * C[j] * E[j];
for (int j = 1; j <= i; ++j) G[i] += G[i - j] * C[j];
// E[i] = 1 + (1/i!) (G[i] + C[i] E[i])
if (i >= 2) {
E[i] = (1 + invFac[i] * G[i]) / (1 - invFac[i] * C[i]);
}
G[i] += C[i] * E[i];
}
// cerr<<"E = "<<E<<endl;
// cerr<<"F = "<<F<<endl;
// cerr<<"G = "<<G<<endl;
*/
Poly E(N + 1), G(N + 1);
for (int i = 1; i <= N; ++i) {
const int w = i & -i;
Poly hs(w << 1);
{
Poly hs0(w << 1), hs1(w << 1);
for (int j = i - w; j < i; ++j) hs0[j - (i - w)] = C[j] * E[j];
for (int j = 0; j < w << 1; ++j) hs1[j] = fac[j];
fft(hs0);
fft(hs1);
for (int j = 0; j < w << 1; ++j) hs[j] += hs0[j] * hs1[j];
}
{
Poly hs0(w << 1), hs1(w << 1);
for (int j = i - w; j < i; ++j) hs0[j - (i - w)] = G[j];
for (int j = 0; j < w << 1 && j <= N; ++j) hs1[j] = C[j];
fft(hs0);
fft(hs1);
for (int j = 0; j < w << 1; ++j) hs[j] += hs0[j] * hs1[j];
}
invFft(hs);
for (int j = i; j < i + w && j <= N; ++j) G[j] += hs[j - (i - w)];
if (i >= 2) {
E[i] = (1 + invFac[i] * G[i]) / (1 - invFac[i] * C[i]);
}
G[i] += C[i] * E[i];
}
// cerr<<"E = "<<E<<endl;
// cerr<<"G = "<<G<<endl;
vector<int> AMax(N + 1, -1);
for (int i = 0; i < N; ++i) AMax[i + 1] = max(AMax[i], A[i]);
vector<int> cut;
for (int i = 0; i <= N; ++i) if (AMax[i] == i - 1) cut.push_back(i);
// cerr<<"cut = "<<cut<<endl;
Mint ans = 0;
for (int j = 0; j < (int)cut.size() - 1; ++j) ans += E[cut[j + 1] - cut[j]];
printf("%u\n", ans.x);
}
return 0;
}
/*
E = [0, 0, 2, 332748121, 725995898, 181498980, 231603911, 962358036, 612851697, 375472763, 352546511]
F = [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
G = [0, 0, 2, 14, 453747527, 816746067, 48008249, 813429526, 637912911, 185889237, 760013790]
*/
这程序好像有点Bug,我给组数据试试?
詳細信息
Test #1:
score: 100
Accepted
time: 12ms
memory: 20416kb
input:
5 2 1 5 3 4
output:
332748123
result:
ok 1 number(s): "332748123"
Test #2:
score: 0
Accepted
time: 9ms
memory: 21792kb
input:
10 4 2 3 1 6 5 9 7 10 8
output:
453747445
result:
ok 1 number(s): "453747445"
Test #3:
score: 0
Accepted
time: 12ms
memory: 22344kb
input:
1 1
output:
0
result:
ok 1 number(s): "0"
Test #4:
score: 0
Accepted
time: 11ms
memory: 20880kb
input:
100 76 73 11 3 60 77 27 93 12 86 40 66 36 2 49 65 17 39 4 20 5 72 67 61 6 35 43 64 58 96 24 68 88 99 94 81 52 13 71 45 32 15 78 85 46 44 55 95 21 53 9 48 26 47 75 33 34 14 70 90 100 62 23 29 22 41 63 84 59 80 38 56 69 25 92 31 50 74 51 87 30 19 91 28 42 89 82 83 97 16 37 54 1 98 10 18 8 57 7 79
output:
508737759
result:
ok 1 number(s): "508737759"
Test #5:
score: 0
Accepted
time: 12ms
memory: 21156kb
input:
1000 644 749 436 202 582 786 943 410 790 963 151 352 625 18 731 435 184 758 614 897 257 917 910 402 589 825 134 380 832 547 77 9 956 229 448 568 769 163 255 364 363 868 148 168 407 663 928 886 450 854 546 994 902 517 793 700 677 936 799 958 475 558 350 650 135 728 768 118 789 231 322 455 221 438 317...
output:
679753239
result:
ok 1 number(s): "679753239"
Test #6:
score: 0
Accepted
time: 142ms
memory: 25304kb
input:
100000 73995 97531 73688 16747 13958 91688 67443 27446 19321 76786 84775 67782 50687 46881 24663 30944 52408 88244 25123 73060 26112 21434 89167 40474 48955 7460 39368 85932 58361 5973 55773 52839 31407 22934 66635 51540 12512 2579 16556 74872 88172 39911 20556 99234 48914 97168 60821 47627 45539 20...
output:
890612070
result:
ok 1 number(s): "890612070"
Test #7:
score: 0
Accepted
time: 143ms
memory: 25104kb
input:
100000 75639 75018 56529 53744 97883 13810 42035 21043 11904 37776 66333 31553 82625 47621 18443 36067 36565 17637 37034 21085 10887 40242 76683 90743 58636 53510 13083 67050 36592 26238 24625 77890 36480 12131 9739 54629 22636 91688 93885 13773 59255 5228 32255 91075 95398 20112 6931 92388 78047 55...
output:
890612070
result:
ok 1 number(s): "890612070"
Test #8:
score: 0
Accepted
time: 144ms
memory: 25364kb
input:
100000 6365 78742 79932 37888 33087 27426 32649 32057 47518 99279 43260 73945 55233 82364 94282 40764 46484 47000 68011 87455 18440 38672 77423 91903 52861 5406 36397 57738 97560 89551 66432 6075 53328 50420 46057 49025 2212 5425 22322 83098 73958 18694 54172 30402 9877 47903 74697 32095 26271 46989...
output:
890612070
result:
ok 1 number(s): "890612070"
Test #9:
score: 0
Accepted
time: 13ms
memory: 21788kb
input:
2 2 1
output:
2
result:
ok 1 number(s): "2"
Test #10:
score: 0
Accepted
time: 14ms
memory: 21532kb
input:
3 3 1 2
output:
332748121
result:
ok 1 number(s): "332748121"
Test #11:
score: 0
Accepted
time: 12ms
memory: 20504kb
input:
4 2 4 1 3
output:
725995898
result:
ok 1 number(s): "725995898"
Test #12:
score: 0
Accepted
time: 13ms
memory: 21028kb
input:
5 4 1 3 5 2
output:
181498980
result:
ok 1 number(s): "181498980"
Test #13:
score: 0
Accepted
time: 9ms
memory: 20372kb
input:
6 5 3 2 1 6 4
output:
231603911
result:
ok 1 number(s): "231603911"
Test #14:
score: 0
Accepted
time: 13ms
memory: 20380kb
input:
7 5 1 4 2 7 3 6
output:
962358036
result:
ok 1 number(s): "962358036"
Test #15:
score: 0
Accepted
time: 14ms
memory: 21044kb
input:
8 4 2 8 7 3 1 6 5
output:
612851697
result:
ok 1 number(s): "612851697"
Test #16:
score: 0
Accepted
time: 12ms
memory: 20372kb
input:
9 2 9 3 6 1 5 4 8 7
output:
375472763
result:
ok 1 number(s): "375472763"
Test #17:
score: 0
Accepted
time: 11ms
memory: 21352kb
input:
10 10 4 5 1 8 9 2 6 3 7
output:
352546511
result:
ok 1 number(s): "352546511"
Test #18:
score: 0
Accepted
time: 142ms
memory: 27748kb
input:
100000 1 2 3 4 6 5 7 9 10 8 12 11 13 14 15 16 17 18 19 22 24 21 25 20 23 27 26 29 28 30 33 31 32 34 35 39 36 37 38 40 42 41 43 44 46 45 48 47 49 50 51 52 53 54 55 56 57 59 58 61 60 62 63 64 65 66 67 69 68 70 72 71 74 73 75 76 78 79 77 80 82 81 83 84 85 86 87 88 90 89 91 92 93 94 95 98 96 97 100 99 1...
output:
210320684
result:
ok 1 number(s): "210320684"
Test #19:
score: 0
Accepted
time: 145ms
memory: 24608kb
input:
100000 1 2 3 4 5 6 7 8 11 9 10 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 36 35 37 38 40 39 41 43 42 45 44 46 47 50 49 48 51 52 53 54 56 55 57 59 58 60 61 62 63 65 64 68 66 67 69 71 70 73 72 75 74 76 78 77 81 79 80 82 83 84 85 86 87 88 89 91 90 92 93 94 95 97 96 99 98 101 1...
output:
270035219
result:
ok 1 number(s): "270035219"
Test #20:
score: 0
Accepted
time: 144ms
memory: 27976kb
input:
100000 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 60 59 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 78 77 79 80 81 82 83 84 85 86 87 88 89 91 90 92 93 95 94 96 97 98 99 100 1...
output:
605008306
result:
ok 1 number(s): "605008306"
Test #21:
score: 0
Accepted
time: 144ms
memory: 25380kb
input:
100000 2 1 4 5 3 7 6 9 8 12 11 10 14 13 15 17 18 16 19 21 20 23 24 22 25 26 27 28 30 29 32 31 34 33 36 37 35 38 39 41 40 43 42 46 44 48 45 47 50 49 51 52 53 54 58 57 55 56 62 60 59 61 67 68 69 65 70 66 63 64 71 72 73 75 74 77 76 79 78 82 81 83 80 86 84 85 87 88 89 93 90 91 92 96 94 95 97 99 98 101 1...
output:
495537170
result:
ok 1 number(s): "495537170"
Test #22:
score: 0
Accepted
time: 147ms
memory: 27956kb
input:
100000 63073 28667 90975 30352 3761 94190 79972 43048 21778 33653 36197 15091 608 19719 989 52387 77322 60994 47624 43439 23816 67364 99594 78540 14592 23284 76709 93993 91131 92333 87166 76515 65187 93200 55352 73418 5417 62724 89224 44808 30562 77028 86252 62548 60089 72547 31563 68303 9267 37119 ...
output:
890612070
result:
ok 1 number(s): "890612070"
Test #23:
score: 0
Accepted
time: 143ms
memory: 24848kb
input:
100000 40205 6086 18833 72928 8675 2421 39538 10358 33411 19062 21955 49175 25071 356 53723 1249 39168 4166 41373 71509 7236 19815 51707 28958 19596 12000 34436 64228 24654 13516 37785 58526 64106 67407 22334 38939 8020 50688 75295 37521 32185 64781 4779 34223 22027 69356 44383 22618 23734 15468 270...
output:
88556225
result:
ok 1 number(s): "88556225"
Test #24:
score: 0
Accepted
time: 142ms
memory: 25332kb
input:
100000 13177 13025 19647 20671 21372 25297 23161 25709 9980 18962 26049 259 15667 5396 19714 25946 7421 19163 8310 5203 7356 22748 21994 16406 8927 4437 13123 11514 25572 20363 2776 8988 6884 10095 13026 16181 3703 321 5459 19285 17305 25129 13822 13685 23863 5389 9401 2056 22066 13275 13133 2023 90...
output:
249803513
result:
ok 1 number(s): "249803513"
Test #25:
score: 0
Accepted
time: 141ms
memory: 25436kb
input:
100000 118 359 351 462 171 279 312 125 315 362 249 425 194 356 226 436 204 122 87 331 398 127 73 278 384 262 377 299 129 242 264 212 215 75 82 392 36 417 407 382 102 210 4 21 130 136 361 461 410 432 319 420 427 292 443 180 261 197 63 370 459 92 151 214 454 158 396 74 376 435 287 386 199 437 348 441 ...
output:
224837018
result:
ok 1 number(s): "224837018"
Test #26:
score: 0
Accepted
time: 147ms
memory: 27812kb
input:
100000 16 60 89 57 102 10 94 39 1 62 49 76 61 82 105 17 42 63 9 47 69 7 31 87 50 21 71 43 3 114 54 25 32 93 115 81 36 118 79 34 58 88 29 91 5 56 66 117 80 104 83 55 108 75 95 18 112 101 113 100 41 96 6 14 35 107 84 85 38 45 33 103 59 22 28 68 8 11 26 98 110 116 2 78 97 12 86 46 111 53 92 15 65 77 44...
output:
831156714
result:
ok 1 number(s): "831156714"
Test #27:
score: 0
Accepted
time: 143ms
memory: 25348kb
input:
100000 10 22 8 13 20 3 9 14 11 12 15 7 19 16 21 17 23 18 1 4 2 6 5 25 24 29 27 28 26 31 30 44 34 32 38 37 41 43 36 39 40 42 46 33 35 45 48 47 53 49 50 51 52 54 57 58 61 60 56 64 65 55 66 63 59 62 68 67 72 70 69 71 77 76 80 74 75 73 81 78 79 82 85 84 83 88 89 86 87 102 106 103 98 93 100 94 96 105 101...
output:
374593697
result:
ok 1 number(s): "374593697"
Test #28:
score: 0
Accepted
time: 145ms
memory: 24888kb
input:
100000 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 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 1...
output:
0
result:
ok 1 number(s): "0"
Test #29:
score: 0
Accepted
time: 141ms
memory: 27712kb
input:
100000 40 48 75 84 93 34 1 12 16 8 22 96 68 55 74 70 59 52 24 19 78 2 63 43 39 18 54 32 89 42 57 47 6 72 28 35 9 95 20 29 67 64 37 71 38 41 46 77 61 88 3 5 87 80 51 65 17 86 82 90 97 98 73 36 44 62 60 49 25 27 83 13 10 11 58 14 66 15 4 7 56 30 33 81 99 31 23 76 85 21 91 92 45 26 50 94 69 79 53 135 1...
output:
577674591
result:
ok 1 number(s): "577674591"
Test #30:
score: 0
Accepted
time: 149ms
memory: 27496kb
input:
100000 3 1 2 5 6 4 10 9 8 7 11 19 17 14 16 13 12 15 18 24 20 22 21 25 23 27 26 28 30 29 31 36 34 35 33 37 32 39 40 38 43 45 41 42 46 44 48 47 50 49 51 53 52 56 54 58 55 59 57 60 62 64 61 63 67 66 65 69 72 68 75 73 74 71 70 98 78 102 77 93 86 92 85 96 79 99 84 83 81 91 87 97 89 94 100 82 103 80 95 10...
output:
516708075
result:
ok 1 number(s): "516708075"
Extra Test:
score: 0
Extra Test Passed