#1033571 (C++17) No.1907 DETERMINATION

提出ソース

結果

問題	No.1907 DETERMINATION
コンテスト
ユーザー	Benjamin Qi
提出日時	2024-12-29 11:56:10
言語	C++17 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	964 ms / 4,000 ms
コード長	16,766 bytes
コンパイル時間	2,568 ms
コンパイル使用メモリ	157,552 KB
最終ジャッジ日時	2025-02-26 17:08:04
ジャッジサーバーID （参考情報）	judge4 / judge2

このコードへのチャレンジ
（要ログイン）

ファイルパターン	結果
sample	AC * 4
other	AC * 63

権限があれば一括ダウンロードができます

コンパイルメッセージ

main.cpp: In function ‘void FileIO::setIn(str)’:
main.cpp:246:28: warning: ignoring return value of ‘FILE* freopen(const char*, const char*, FILE*)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  246 | void setIn(str s) { freopen(s.c_str(), "r", stdin); }
      |                     ~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~
main.cpp: In function ‘void FileIO::setOut(str)’:
main.cpp:247:29: warning: ignoring return value of ‘FILE* freopen(const char*, const char*, FILE*)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  247 | void setOut(str s) { freopen(s.c_str(), "w", stdout); }
      |                      ~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~

ソースコード

raw source code

#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <climits>
#include <cmath>
#include <complex>
#include <cstring>
#include <functional>
#include <iomanip>
#include <iostream>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <vector>
using namespace std;

using ll = long long;
using db = long double;  // or double, if TL is tight
using str = string;      // yay python!

// pairs
using pi = pair<int, int>;
using pl = pair<ll, ll>;
using pd = pair<db, db>;
#define mp make_pair
#define f first
#define s second

#define tcT template <class T
#define tcTU tcT, class U
// ^ lol this makes everything look weird but I'll try it
tcT > using V = vector<T>;
tcT, size_t SZ > using AR = array<T, SZ>;
using vi = V<int>;
using vb = V<bool>;
using vl = V<ll>;
using vd = V<db>;
using vs = V<str>;
using vpi = V<pi>;
using vpl = V<pl>;
using vpd = V<pd>;

// vectors
#define sz(x) int(size(x))
#define bg(x) begin(x)
#define all(x) bg(x), end(x)
#define rall(x) rbegin(x), rend(x)
#define sor(x) sort(all(x))
#define rsz resize
#define ins insert
#define pb push_back
#define eb emplace_back
#define ft front()
#define bk back()

#define lb lower_bound
#define ub upper_bound
tcT > int lwb(const V<T> &a, const T &b) { return int(lb(all(a), b) - bg(a)); }
tcT > int upb(const V<T> &a, const T &b) { return int(ub(all(a), b) - bg(a)); }

// loops
#define FOR(i, a, b) for (int i = (a); i < (b); ++i)
#define F0R(i, a) FOR(i, 0, a)
#define ROF(i, a, b) for (int i = (b)-1; i >= (a); --i)
#define R0F(i, a) ROF(i, 0, a)
#define rep(a) F0R(_, a)
#define each(a, x) for (auto &a : x)

const int MOD = 998244353;  // 1e9+7;
const int MX = (int)2e5 + 5;
const ll BIG = 1e18;  // not too close to LLONG_MAX
const db PI = acos((db)-1);
const int dx[4]{1, 0, -1, 0}, dy[4]{0, 1, 0, -1};  // for every grid problem!!
mt19937 rng((uint32_t)chrono::steady_clock::now().time_since_epoch().count());
template <class T> using pqg = priority_queue<T, vector<T>, greater<T>>;

// bitwise ops
// also see https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html
constexpr int pct(int x) { return __builtin_popcount(x); }  // # of bits set
constexpr int bits(int x) {  // assert(x >= 0); // make C++11 compatible until
	                         // USACO updates ...
	return x == 0 ? 0 : 31 - __builtin_clz(x);
}  // floor(log2(x))
constexpr int p2(int x) { return 1 << x; }
constexpr int msk2(int x) { return p2(x) - 1; }

ll cdiv(ll a, ll b) {
	return a / b + ((a ^ b) > 0 && a % b);
}  // divide a by b rounded up
ll fdiv(ll a, ll b) {
	return a / b - ((a ^ b) < 0 && a % b);
}  // divide a by b rounded down

tcT > bool ckmin(T &a, const T &b) {
	return b < a ? a = b, 1 : 0;
}  // set a = min(a,b)
tcT > bool ckmax(T &a, const T &b) {
	return a < b ? a = b, 1 : 0;
}  // set a = max(a,b)

tcTU > T fstTrue(T lo, T hi, U f) {
	++hi;
	assert(lo <= hi);  // assuming f is increasing
	while (lo < hi) {  // find first index such that f is true
		T mid = lo + (hi - lo) / 2;
		f(mid) ? hi = mid : lo = mid + 1;
	}
	return lo;
}
tcTU > T lstTrue(T lo, T hi, U f) {
	--lo;
	assert(lo <= hi);  // assuming f is decreasing
	while (lo < hi) {  // find first index such that f is true
		T mid = lo + (hi - lo + 1) / 2;
		f(mid) ? lo = mid : hi = mid - 1;
	}
	return lo;
}
tcT > void remDup(vector<T> &v) {  // sort and remove duplicates
	sort(all(v));
	v.erase(unique(all(v)), end(v));
}
tcTU > void safeErase(T &t, const U &u) {
	auto it = t.find(u);
	assert(it != end(t));
	t.erase(it);
}

inline namespace IO {
#define SFINAE(x, ...)                                                         \
	template <class, class = void> struct x : std::false_type {};              \
	template <class T> struct x<T, std::void_t<__VA_ARGS__>> : std::true_type {}

SFINAE(DefaultI, decltype(std::cin >> std::declval<T &>()));
SFINAE(DefaultO, decltype(std::cout << std::declval<T &>()));
SFINAE(IsTuple, typename std::tuple_size<T>::type);
SFINAE(Iterable, decltype(std::begin(std::declval<T>())));

template <auto &is> struct Reader {
	template <class T> void Impl(T &t) {
		if constexpr (DefaultI<T>::value) is >> t;
		else if constexpr (Iterable<T>::value) {
			for (auto &x : t) Impl(x);
		} else if constexpr (IsTuple<T>::value) {
			std::apply([this](auto &...args) { (Impl(args), ...); }, t);
		} else static_assert(IsTuple<T>::value, "No matching type for read");
	}
	template <class... Ts> void read(Ts &...ts) { ((Impl(ts)), ...); }
};

template <class... Ts> void re(Ts &...ts) { Reader<cin>{}.read(ts...); }
#define def(t, args...)                                                        \
	t args;                                                                    \
	re(args);

template <auto &os, bool debug, bool print_nd> struct Writer {
	string comma() const { return debug ? "," : ""; }
	template <class T> constexpr char Space(const T &) const {
		return print_nd && (Iterable<T>::value or IsTuple<T>::value) ? '\n'
		                                                             : ' ';
	}
	template <class T> void Impl(T const &t) const {
		if constexpr (DefaultO<T>::value) os << t;
		else if constexpr (Iterable<T>::value) {
			if (debug) os << '{';
			int i = 0;
			for (auto &&x : t)
				((i++) ? (os << comma() << Space(x), Impl(x)) : Impl(x));
			if (debug) os << '}';
		} else if constexpr (IsTuple<T>::value) {
			if (debug) os << '(';
			std::apply(
			    [this](auto const &...args) {
				    int i = 0;
				    (((i++) ? (os << comma() << " ", Impl(args)) : Impl(args)),
				     ...);
			    },
			    t);
			if (debug) os << ')';
		} else static_assert(IsTuple<T>::value, "No matching type for print");
	}
	template <class T> void ImplWrapper(T const &t) const {
		if (debug) os << "\033[0;31m";
		Impl(t);
		if (debug) os << "\033[0m";
	}
	template <class... Ts> void print(Ts const &...ts) const {
		((Impl(ts)), ...);
	}
	template <class F, class... Ts>
	void print_with_sep(const std::string &sep, F const &f,
	                    Ts const &...ts) const {
		ImplWrapper(f), ((os << sep, ImplWrapper(ts)), ...), os << '\n';
	}
	void print_with_sep(const std::string &) const { os << '\n'; }
};

template <class... Ts> void pr(Ts const &...ts) {
	Writer<cout, false, true>{}.print(ts...);
}
template <class... Ts> void ps(Ts const &...ts) {
	Writer<cout, false, true>{}.print_with_sep(" ", ts...);
}
}  // namespace IO

inline namespace Debug {
template <typename... Args> void err(Args... args) {
	Writer<cerr, true, false>{}.print_with_sep(" | ", args...);
}
template <typename... Args> void errn(Args... args) {
	Writer<cerr, true, true>{}.print_with_sep(" | ", args...);
}

void err_prefix(str func, int line, string args) {
	cerr << "\033[0;31m\u001b[1mDEBUG\033[0m"
	     << " | "
	     << "\u001b[34m" << func << "\033[0m"
	     << ":"
	     << "\u001b[34m" << line << "\033[0m"
	     << " - "
	     << "[" << args << "] = ";
}

#ifdef LOCAL
#define dbg(args...) err_prefix(__FUNCTION__, __LINE__, #args), err(args)
#define dbgn(args...) err_prefix(__FUNCTION__, __LINE__, #args), errn(args)
#else
#define dbg(...)
#define dbgn(args...)
#endif

const auto beg_time = std::chrono::high_resolution_clock::now();
// https://stackoverflow.com/questions/47980498/accurate-c-c-clock-on-a-multi-core-processor-with-auto-overclock?noredirect=1&lq=1
double time_elapsed() {
	return chrono::duration<double>(std::chrono::high_resolution_clock::now() -
	                                beg_time)
	    .count();
}
}  // namespace Debug

inline namespace FileIO {
void setIn(str s) { freopen(s.c_str(), "r", stdin); }
void setOut(str s) { freopen(s.c_str(), "w", stdout); }
void setIO(str s = "") {
	cin.tie(0)->sync_with_stdio(0);  // unsync C / C++ I/O streams
	cout << fixed << setprecision(12);
	// cin.exceptions(cin.failbit);
	// throws exception when do smth illegal
	// ex. try to read letter into int
	if (sz(s)) setIn(s + ".in"), setOut(s + ".out");  // for old USACO
}
}  // namespace FileIO

/**
 * Description: modular arithmetic operations
 * Source:
 * KACTL
 * https://codeforces.com/blog/entry/63903
 * https://codeforces.com/contest/1261/submission/65632855 (tourist)
 * https://codeforces.com/contest/1264/submission/66344993 (ksun)
 * also see https://github.com/ecnerwala/cp-book/blob/master/src/modnum.hpp
 * (ecnerwal) Verification: https://open.kattis.com/problems/modulararithmetic
 */

template <int MOD, int RT> struct mint {
	static const int mod = MOD;
	static constexpr mint rt() { return RT; }  // primitive root for FFT
	int v;
	explicit operator int() const {
		return v;
	}  // explicit -> don't silently convert to int
	mint() : v(0) {}
	mint(ll _v) {
		v = int((-MOD < _v && _v < MOD) ? _v : _v % MOD);
		if (v < 0) v += MOD;
	}
	bool operator==(const mint &o) const { return v == o.v; }
	friend bool operator!=(const mint &a, const mint &b) { return !(a == b); }
	friend bool operator<(const mint &a, const mint &b) { return a.v < b.v; }
	friend istream &operator>>(istream &is, mint &a) {
		ll x;
		is >> x;
		a = mint(x);
		return is;
	}
	friend ostream &operator<<(ostream &os, mint a) {
		os << int(a);
		return os;
	}

	mint &operator+=(const mint &o) {
		if ((v += o.v) >= MOD) v -= MOD;
		return *this;
	}
	mint &operator-=(const mint &o) {
		if ((v -= o.v) < 0) v += MOD;
		return *this;
	}
	mint &operator*=(const mint &o) {
		v = int((ll)v * o.v % MOD);
		return *this;
	}
	mint &operator/=(const mint &o) { return (*this) *= inv(o); }
	friend mint pow(mint a, ll p) {
		mint ans = 1;
		assert(p >= 0);
		for (; p; p /= 2, a *= a)
			if (p & 1) ans *= a;
		return ans;
	}
	friend mint inv(const mint &a) {
		assert(a.v != 0);
		return pow(a, MOD - 2);
	}

	mint operator-() const { return mint(-v); }
	mint &operator++() { return *this += 1; }
	mint &operator--() { return *this -= 1; }
	friend mint operator+(mint a, const mint &b) { return a += b; }
	friend mint operator-(mint a, const mint &b) { return a -= b; }
	friend mint operator*(mint a, const mint &b) { return a *= b; }
	friend mint operator/(mint a, const mint &b) { return a /= b; }
};

using mi = mint<MOD, 5>;  // 5 is primitive root for both common mods
using vmi = V<mi>;
using pmi = pair<mi, mi>;
using vpmi = V<pmi>;

V<vmi> scmb;  // small combinations
void genComb(int SZ) {
	scmb.assign(SZ, vmi(SZ));
	scmb[0][0] = 1;
	FOR(i, 1, SZ)
	F0R(j, i + 1) scmb[i][j] = scmb[i - 1][j] + (j ? scmb[i - 1][j - 1] : 0);
}

void hessenberg(V<vmi> &v) {
	// https://www.phys.uri.edu/nigh/NumRec/bookfpdf/f11-5.pdf#page=3
	int N = sz(v);
	F0R(c, N) {
		bool found = 0;
		FOR(r, c + 1, N) if (v.at(r).at(c) != 0) {
			swap(v.at(c + 1), v.at(r));
			F0R(k, N) swap(v.at(k).at(c + 1), v.at(k).at(r));
			found = true;
			break;
		}
		if (!found) continue;
		mi iv = 1 / v.at(c + 1).at(c);
		FOR(r, c + 2, N) {
			mi quo = v.at(r).at(c) * iv;
			FOR(co, c, N) v.at(r).at(co) -= quo * v.at(c + 1).at(co);
			F0R(ro, N) v.at(ro).at(c + 1) += quo * v.at(ro).at(r);
			assert(v.at(r).at(c) == 0);
		}
	}
}

/**
 * Description: Basic poly ops including division. Can replace \texttt{T} with
 * double, complex. Source: Own. Also see
 * https://github.com/kth-competitive-programming/kactl/blob/master/content/numerical/PolyInterpolate.h
 * https://github.com/ecnerwala/icpc-book/blob/master/content/numerical/fft.cpp
 * Verification: see FFT
 */

// #include "../../number-theory (11.1)/Modular Arithmetic/ModInt.h"

using T = mi;
using poly = V<T>;
void remz(poly &p) {
	while (sz(p) && p.bk == T(0)) p.pop_back();
}
poly REMZ(poly p) {
	remz(p);
	return p;
}
poly rev(poly p) {
	reverse(all(p));
	return p;
}
poly shift(poly p, int x) {
	if (x >= 0) p.insert(begin(p), x, 0);
	else assert(sz(p) + x >= 0), p.erase(begin(p), begin(p) - x);
	return p;
}
poly RSZ(const poly &p, int x) {
	if (x <= sz(p)) return poly(begin(p), begin(p) + x);
	poly q = p;
	q.rsz(x);
	return q;
}
T eval(const poly &p, T x) {  // evaluate at point x
	T res = 0;
	R0F(i, sz(p)) res = x * res + p[i];
	return res;
}
poly dif(const poly &p) {  // differentiate
	poly res;
	FOR(i, 1, sz(p)) res.pb(T(i) * p[i]);
	return res;
}
poly integ(const poly &p) {  // integrate
	static poly invs{0, 1};
	for (int i = sz(invs); i <= sz(p); ++i) invs.pb(-MOD / i * invs[MOD % i]);
	poly res(sz(p) + 1);
	F0R(i, sz(p)) res[i + 1] = p[i] * invs[i + 1];
	return res;
}

poly &operator+=(poly &l, const poly &r) {
	l.rsz(max(sz(l), sz(r)));
	F0R(i, sz(r)) l[i] += r[i];
	return l;
}
poly &operator-=(poly &l, const poly &r) {
	l.rsz(max(sz(l), sz(r)));
	F0R(i, sz(r)) l[i] -= r[i];
	return l;
}
poly &operator*=(poly &l, const T &r) {
	each(t, l) t *= r;
	return l;
}
poly &operator/=(poly &l, const T &r) {
	each(t, l) t /= r;
	return l;
}
poly operator+(poly l, const poly &r) { return l += r; }
poly operator-(poly l, const poly &r) { return l -= r; }
poly operator-(poly l) {
	each(t, l) t *= -1;
	return l;
}
poly operator*(poly l, const T &r) { return l *= r; }
poly operator*(const T &r, const poly &l) { return l * r; }
poly operator/(poly l, const T &r) { return l /= r; }
poly operator*(const poly &l, const poly &r) {
	if (!min(sz(l), sz(r))) return {};
	poly x(sz(l) + sz(r) - 1);
	F0R(i, sz(l)) F0R(j, sz(r)) x[i + j] += l[i] * r[j];
	return x;
}
poly &operator*=(poly &l, const poly &r) { return l = l * r; }

pair<poly, poly> quoRemSlow(poly a, poly b) {
	remz(a);
	remz(b);
	assert(sz(b));
	T lst = b.bk, B = T(1) / lst;
	each(t, a) t *= B;
	each(t, b) t *= B;
	poly q(max(sz(a) - sz(b) + 1, 0));
	for (int dif; (dif = sz(a) - sz(b)) >= 0; remz(a)) {
		q[dif] = a.bk;
		F0R(i, sz(b)) a[i + dif] -= q[dif] * b[i];
	}
	each(t, a) t *= lst;
	return {q, a};  // quotient, remainder
}
poly operator%(const poly &a, const poly &b) { return quoRemSlow(a, b).s; }
/**poly operator/(const poly& a, const poly& b) {
    return quoRemSlow(a,b).f; }
poly a = {1,3,5,8,6,0,0,0,0}, b = {1,5,1};
ps(quoRemSlow(a,b)); a = 2*a, b = 2*b; ps(quoRemSlow(a,b));
poly gcd(poly a, poly b) { return b == poly{} ? a : gcd(b,a%b); }*/
T resultant(poly a, poly b) {  // R(A,B)
	// =b_m^n*prod_{j=1}^mA(mu_j)
	// =b_m^na_n^m*prod_{i=1}^nprod_{j=1}^m(mu_j-lambda_i)
	// =(-1)^{mn}a_n^m*prod_{i=1}^nB(lambda_i)
	// =(-1)^{nm}R(B,A)
	// Also, R(A,B)=b_m^{deg(A)-deg(A-CB)}R(A-CB,B)
	int ad = sz(a) - 1, bd = sz(b) - 1;
	if (bd <= 0) return bd < 0 ? 0 : pow(b.bk, ad);
	int pw = ad;
	a = a % b;
	pw -= (ad = sz(a) - 1);
	return resultant(b, a) * pow(b.bk, pw) * T((bd & ad & 1) ? -1 : 1);
}

vmi charpoly(V<vmi> A) {  // det(xI - A)
	int N = sz(A);
	hessenberg(A);
	V<vmi> charpoly;
	charpoly.pb({1});
	F0R(i, N) {
		charpoly.pb(charpoly.bk * vmi{-A.at(i).at(i), 1});
		F0R(j, i) {
			charpoly.bk +=
			    ((i - j) & 1 ? -1 : 1) * charpoly.at(j) * -A.at(j).at(i);
		}
		if (i + 1 < N) F0R(j, i + 1) charpoly.at(j) *= -A.at(i + 1).at(i);
	}
	return charpoly.bk;
}

vmi det(V<vmi> A, V<vmi> B) {  // det(A + Bx)
	int N = sz(A);
	int off = 0;
	mi prod = 1;
	auto swap_cols = [&](int c1, int c2) {
		assert(c1 != c2);
		F0R(r, N) {
			swap(A.at(r).at(c1), A.at(r).at(c2));
			swap(B.at(r).at(c1), B.at(r).at(c2));
		}
		prod *= -1;
	};
	auto sub_rows = [&](int r1, int r2, mi coef) {
		assert(r1 != r2);
		F0R(c, N) {
			A.at(r1).at(c) -= coef * A.at(r2).at(c);
			B.at(r1).at(c) -= coef * B.at(r2).at(c);
		}
	};
	auto mul_row = [&](int r, mi coef) {
		F0R(c, N) {
			A.at(r).at(c) *= coef;
			B.at(r).at(c) *= coef;
		}
	};
	F0R(r, N) {
		while (B[r][r] == 0) {
			FOR(c, r + 1, N) if (B.at(r).at(c) != 0) {
				swap_cols(r, c);
				break;
			}
			if (B[r][r] != 0) break;
			swap(A.at(r), B.at(r));
			++off;
			// assert(false);
			if (off > N) return vmi(N + 1);
			F0R(ro, r) sub_rows(r, ro, B.at(r).at(ro));
		}
		prod *= B[r][r];
		mul_row(r, 1 / B[r][r]);
		F0R(ro, N) if (r != ro) sub_rows(ro, r, B.at(ro).at(r));
	}
	F0R(i, N) assert(B[i][i] == 1);
	each(a, A) each(b, a) b *= -1;
	auto cA = charpoly(A);
	vmi ret(N + 1);
	F0R(i, N + 1 - off) ret.at(i) = cA.at(i + off) * prod;
	return ret;
}

int main() {
	// read read read
	setIO();
	def(int, N);
	V<vmi> A(N, vmi(N)), B(N, vmi(N));
	re(A, B);
	auto ret = det(A, B);
	each(t, ret) ps(t);
	// vi P(N);
	// re(P);
	// V<vmi> A(N - 1, vmi(N - 1)), B(N - 1, vmi(N - 1));
	// auto edge = [&](V<vmi> &v, int x, int y) {
	// 	if (x < N - 1) ++v.at(x).at(x);
	// 	if (x < N - 1 && y < N - 1) --v.at(x).at(y);
	// };
	// F0R(i, N) FOR(j, i + 1, N) {
	// 	if (P[i] > P[j]) {
	// 		edge(B, i, j);
	// 		edge(B, j, i);
	// 	} else {
	// 		edge(A, i, j);
	// 		edge(A, j, i);
	// 	}
	// }
	// dbg(A);
	// dbg(B);
	// auto ret = det(A, B);
	// ps(ret);

	// you should actually read the stuff at the bottom
}

/* stuff you should look for
 * int overflow, array bounds
 * special cases (n=1?)
 * do smth instead of nothing and stay organized
 * WRITE STUFF DOWN
 * DON'T GET STUCK ON ONE APPROACH
 */

yukicoder

結果

コンパイルメッセージ

ソースコード