Finding the roots of a special kind of polynomial

find_root.cpp

// Mivik 2021.3.1
#include <mivik.h>

#include <algorithm>
#include <cassert>
#include <cctype>
#include <climits>
#include <random>
#include <vector>

MI cin;

typedef long long qe;

const int mod = 998244353;

inline int add(int x, int y) { return (x += y) >= mod? x - mod: x; }
inline void Add(int &x, int y) { if ((x += y) >= mod) x -= mod; }
inline int sub(int x, int y) { return (x -= y) < 0? x + mod: x; }
inline void Sub(int &x, int y) { if ((x -= y) < 0) x += mod; }
inline qe ksm(qe x, int p = mod - 2) {
	qe ret = 1;
	for (; p; p >>= 1, (x *= x) %= mod) if (p & 1) (ret *= x) %= mod;
	return ret;
}
// coco: preserve_begin
inline void ntt(int *v, int len, bool rev) {
	static const int N = 17;
	static struct _
	{ int rt[1 << N], inv[N + 1]; _() {
		for (int i = 0, j = 1, q = 0; i <= N; ++i, j = (q = j) << 1) {
			inv[i] = ksm(j);
			const int del = ksm(3, (mod - 1) / j);
			for (int k = 0, c = 1; k < q; ++k, c = (qe)c * del % mod) rt[q | k] = c;
		}
	} } _;
	static int last, r[1 << N];
	if (last != len) {
		for (int i = 1; i < len; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1)? (len >> 1): 0);
		last = len;
	}
	for (int i = 1; i < len; ++i)
		if (i < r[i]) std::swap(v[i], v[r[i]]);
	assert(!(len & (len - 1)));
	for (int i = 1, q = 2; i < len; i = q, q <<= 1) {
		for (int j = 0; j < len; j += q)
			for (int k = 0; k < i; ++k) {
				const int x = v[j | k], y = (qe)v[i | j | k] * _.rt[i | k] % mod;
				v[j | k] = add(x, y);
				v[i | j | k] = sub(x, y);
			}
	}
	if (!rev) return;
	std::reverse(v + 1, v + len);
	for (int i = 0, r = _.inv[__builtin_ctz(len)]; i < len; ++i) v[i] = (qe)v[i] * r % mod;
}
struct poly : public std::vector<int> {
#define v (*this)
#define F(n) \
	inline poly n(int len = -1) const { return poly().from_##n(v, len); } \
	inline poly& from_##n(const poly &a, int len = -1)
#define I(n) \
	inline poly n() const { return poly(v).inplace_##n(); } \
	inline poly& inplace_##n()

	static inline int round_up(int x) { return (x & (x - 1))? (1 << (32 - __builtin_clz(x))): x; }

	poly() {}
	poly(const std::initializer_list<int> &list): std::vector<int>(list) {}
	poly(int x) { push_back(x); }
	inline int deg() const { return empty()? 0: (size() - 1); }
	inline poly& from(const poly &a, int n) {
		resize(n); const int len = std::min(n, (int)a.size());
		std::copy(a.begin(), a.begin() + len, begin());
		std::fill(begin() + len, end(), 0);
		return v;
	}
	inline poly take(int n) const { return poly().from(v, n); }
	inline poly& ntt(int len, bool rev) { resize(round_up(len), 0); ::ntt(data(), size(), rev); return v; }
	inline poly& ntt(bool rev) { return ntt(size(), rev); }
	inline poly& trim() { while (!empty() && !back()) pop_back(); return v; }
	inline poly operator+(const poly &a) const {
		poly ret; ret.resize(std::max(size(), a.size()));
		for (int i = std::min(size(), a.size()) - 1; ~i; --i) ret[i] = add(v[i], a[i]);
		if (size() > a.size()) std::copy(begin() + a.size(), end(), ret.begin() + a.size());
		else std::copy(a.begin() + size(), a.end(), ret.begin() + size());
		return ret.trim();
	}
	inline poly operator-() const {
		poly r; r.resize(size());
		for (int i = size() - 1; ~i; --i) r[i] = sub(0, v[i]);
		return r;
	}
	inline poly& imul(const poly &a, int l) {
		static poly tmp;
		if (empty() || a.empty()) return clear(), v;
		const int len = round_up(l);
		(tmp = a).ntt(len, 0); ntt(len, 0);
		for (int i = 0; i < len; ++i) v[i] = (qe)v[i] * tmp[i] % mod;
		return ntt(len, 1), resize(l), v;
	}
	inline poly& operator*=(const poly &a) { return imul(a, size() + a.size() - 1); }
	inline poly operator*(const poly &t) const { return poly(v) *= t; }
	I(integral) {
		if (empty()) return clear(), v; push_back(0);
		for (int i = size() - 1; i; --i) v[i] = (qe)v[i - 1] * ksm(i) % mod;
		v[0] = 0; return v;
	}
	I(derivative) {
		if (size() <= 1) return clear(), v;
		for (int i = 1; i < size(); ++i) v[i - 1] = (qe)v[i] * i % mod;
		return pop_back(), v;
	}
	I(reverse) { return std::reverse(begin(), end()), v; }
	template<class Func>
	inline poly& newton(const poly &a, int len, int initial, const Func &trans) {
		static poly tmp;
		if (len == -1) len = a.size();
		if (a.empty()) return clear(), v;
		const int lim = round_up(len);
		clear(); push_back(initial);
		for (int l = 2; l <= lim; l <<= 1) { tmp.from(a, l); trans(l, tmp); }
		return resize(len), v;
	}
	F(inv) {
		assert(!a.empty());
		return newton(a, len, ksm(a[0]), [this](int l, poly &tmp) {
			const int q = l << 1;
			tmp.ntt(q, 0); ntt(q, 0);
			for (int i = 0; i < q; ++i)
				v[i] = (qe)v[i] * sub(2, (qe)v[i] * tmp[i] % mod) % mod;
			ntt(1); resize(l);
		});
	}
	F(ln) {
		if (len == -1) len = a.size();
		assert(!a.empty() && a[0] == 1);
		from_inv(a, len); v *= a.derivative();
		return resize(len - 1), inplace_integral();
	}
	F(exp) {
		static poly t2;
		if (a.empty()) return clear(), push_back(1), v;
		assert(!a[0]);
		return newton(a, len, 1, [this](int l, poly &t1) {
			t2.from_ln(v, l);
			for (int i = 0; i < l; ++i) t2[i] = sub(t1[i], t2[i]);
			Add(t2[0], 1);
			v *= t2; resize(l);
		});
	}
	inline poly& from_division(const poly &a, const poly &b) {
		const int len = a.size() - b.size() + 1; assert(len >= 0);
		resize(len); std::reverse_copy(a.end() - len, a.end(), begin());
		v *= b.reverse().inv(len);
		return resize(len), inplace_reverse(), trim();
	}
	inline poly& from_remainder(const poly &a, const poly &b, const poly &q) {
		const int m = b.size() - 1;
		from(b, m) *= q.take(m); resize(m);
		for (int i = 0; i < m; ++i) v[i] = sub(a[i], v[i]);
		return v.trim();
	}
	inline poly operator/(const poly &a) const { return poly().from_division(v, a); }
	inline poly operator%(const poly &a) const { return poly().from_remainder(v, a, v / a); }
	inline std::pair<poly, poly> divide(const poly &a) const { std::pair<poly, poly> ret; ret.second.from_remainder(v, a, ret.first.from_division(v, a)); return ret; }
	inline friend std::ostream& operator<(std::ostream &out, const poly &t) { out < '['; for (int d : t) out < ' ' < d; return out < " ]"; }
#undef v
#undef F
#undef I
};
// coco: preserve_end

const int T = (mod - 1) / 2;

inline poly monic(const poly &A) {
	poly r; r.resize(A.size()); const int inv = ksm(A.back());
	for (int i = r.size() - 1; ~i; --i) r[i] = (qe)A[i] * inv % mod;
	return r;
}
struct mat { poly v[2][2]; };
inline mat operator*(const mat &A, const mat &B) {
#define C(a, b) (A.v[a][0] * B.v[0][b] + A.v[a][1] * B.v[1][b])
	return { C(0, 0), C(0, 1), C(1, 0), C(1, 1) };
#undef C
}
inline std::pair<poly, poly> operator*(const mat &A, const std::pair<poly, poly> &B) {
	return {
		A.v[0][0] * B.first + A.v[0][1] * B.second,
		A.v[1][0] * B.first + A.v[1][1] * B.second
	};
}
inline poly operator>>(const poly &A, size_t p) {
	if (p >= A.size()) return {};
	poly r; r.resize(A.size() - p);
	std::copy(A.begin() + p, A.end(), r.begin());
	return r;
}
inline mat half_gcd(const poly &A, const poly &B) {
	assert(A.size() > B.size());
	const int m = (A.deg() + 1) >> 1;
	if (B.empty() || B.deg() < m) return { 1, 0, 0, 1 };
	auto M1 = half_gcd(A >> m, B >> m);
	const auto [s, t] = M1 * std::make_pair(A, B);
	if (t.deg() < m) return M1;
	const auto [q, r] = s.divide(t);
	const int k = (m << 1) - t.deg();
	mat M2 = { 0, 1, 1, -q };
	return half_gcd(t >> k, r >> k) * M2 * M1;
}
inline mat gcd_mat(const poly &A, const poly &B) {
	if (A.empty()) return { 0, 1, 1, 0 };
	if (B.empty()) return { 1, 0, 0, 1 };
	if (A.size() <= B.size()) {
		const auto [q, r] = B.divide(A);
		return gcd_mat(A, r) * (mat){ 1, 0, -q, 1 };
	} else {
		const auto M = half_gcd(A, B);
		const auto [s, t] = M * std::make_pair(A, B);
		return gcd_mat(t, s) * (mat){ 0, 1, 1, 0 } * M;
	}
}
inline poly gcd(const poly &A, const poly &B) {
	if (A.empty()) return B;
	if (B.empty()) return A;
	auto M = gcd_mat(A, B);
	return monic(M.v[0][0] * A + M.v[0][1] * B);
}

poly F;
std::vector<int> fac, ifac;
inline poly diff(const poly &A, int k) {
	if (k > A.size()) return {};
	poly r; r.resize(A.size() - k);
	for (int i = r.size() - 1; ~i; --i) r[i] = (qe)A[i + k] * fac[i + k] % mod * ifac[i] % mod;
	return r;
}
inline void init(int n) {
	fac.resize(n + 1); ifac.resize(n + 1);
	for (int i = fac[0] = ifac[0] = 1; i <= n; ++i) fac[i] = (qe)fac[i - 1] * i % mod;
	ifac[n] = ksm(fac[n]); for (int i = n; i > 1; --i) ifac[i - 1] = (qe)ifac[i] * i % mod;
}
void offset(poly &A, int c) {
	static poly F, G;
	const int n = A.size();
	F.resize(n); G.resize(n);
	for (int i = 0, cur = 1; i < n; ++i) {
		G[i] = (qe)A[i] * fac[i] % mod;
		F[n - i - 1] = (qe)cur * ifac[i] % mod;
		cur = (qe)cur * c % mod;
	}
	F *= G; assert(F.size() == n * 2 - 1);
	for (int i = 0; i < n; ++i) A[i] = (qe)F[i + n - 1] * ifac[i] % mod;
}
inline poly ediv(const poly &A, const poly &B) { // Cond: A*Q=B, calculates Q
	auto r = A * B.inv(A.size());
	r.resize(A.size());
	return r.trim();
}
inline poly x_k_mod(int k, const poly &A) {
	if (k < A.size()) { poly r; r.resize(k + 1); r[k] = 1; return r; }
	auto try_mod = [&A](poly &B) { if (B.size() >= A.size()) (B = B % A).trim(); };
	poly cur = { 0, 1 }, ret = { 1 };
	for (qe p = k; p; p >>= 1, try_mod(cur *= cur)) if (p & 1) try_mod(ret *= cur);
	return ret;
}
inline poly gcd(const poly &A) {
	auto B = x_k_mod(T, A); Sub(B[0], 1);
	B.trim(); return gcd(monic(A), B);
}

inline poly pow(const poly &A, int t) {
	if (t == 1) return A;
	poly r; r.resize(A.deg() * t + 1); const int w = ksm(A[0]);
	for (int i = A.size() - 1; ~i; --i) r[i] = (qe)A[i] * w % mod;
	r = r.ln(); for (int &v : r) v = (qe)v * t % mod;
	return r.exp();
}
std::vector<int> solve(poly A, int imp = 1) {
	static std::mt19937 rng(std::random_device{}());
	static std::uniform_int_distribution<int> dist(1, mod - 1);

	if (A.size() == 1) return {};
	if (A.size() == 2) {
		static int cnt; cout < (cnt += imp) < " found\n";
		return { (int)((qe)sub(0, A[0]) * ksm(A[1]) % mod) };
	}
	int off = 0;
	while (1) {
		const int c = dist(rng); offset(A, c); Add(off, c);
		const auto g = gcd(A);
		if (g.size() == 1 || g == A) continue;
		int coe = 1; A = ediv(A, g);
		if (A.size() > g.size() && (A % g).empty()) {
			int l = 1, r = A.deg() / g.deg(), ret = 0;
			while (l <= r) {
				const int mid = (l + r) >> 1;
				if ((A % pow(g, mid)).empty()) l = (ret = mid) + 1;
				else r = mid - 1;
			}
			assert(ret); coe += ret; A = ediv(A, pow(g, ret));
		}
		auto p1 = solve(A, imp), p2 = solve(g, imp * coe);
		p1.reserve(p1.size() + p2.size() * coe);
		while (coe--) p1.insert(p1.end(), p2.begin(), p2.end());
		for (int &v : p1) Add(v, off);
		return p1;
	}
}
int main() {
	const int n = R; init(n); F.resize(n + 1);
	for (int &v : F) cin > v;
	auto rts = solve(F); assert(rts.size() == n);
	std::sort(rts.begin(), rts.end());
	for (int i = 0; i < rts.size(); ++i) cout < rts[i] < " \n"[i == rts.size() - 1];
}