#951833 (C++17) No.2633 Subsequence Combination Score

提出ソース
結果

問題	No.2633 Subsequence Combination Score
コンテスト
ユーザー	ゆにぽけ
提出日時	2024-02-17 00:07:47
言語	C++17 (gcc 13.3.0 + boost 1.87.0)
結果	AC
実行時間	509 ms / 2,000 ms
コード長	15,328 bytes
コンパイル時間	1,920 ms
コンパイル使用メモリ	147,932 KB
最終ジャッジ日時	2025-02-19 15:11:10
ジャッジサーバーID （参考情報）	judge3 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン	結果
sample	AC * 3
other	AC * 38
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <iostream>
#include <vector>
#include <algorithm>
#include <array>
#include <iterator>
#include <string>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <cmath>
#include <ctime>
#include <iomanip>
#include <numeric>
#include <stack>
#include <queue>
#include <map>
#include <unordered_map>
#include <set>
#include <unordered_set>
#include <bitset>
#include <random>
#include <utility>
#include <functional>
using namespace std;
template<class T>
inline bool chmin(T &a,const T &b)
{
	if(a > b)
	{
		a = b;
		return true;
	}
	return false;
}
template<class T>
inline bool chmax(T &a,const T &b)
{
	if(a < b)
	{
		a = b;
		return true;
	}
	return false;
}
template<class T>
void print(const vector<T> &V)
{
	for(int i = 0;i < (int)V.size();i++)
	{
		cerr << V[i] << (i + 1 == (int)V.size() ? "\n":" ");
	}
}
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdlib>
#include <cassert>
#include <cmath>
using namespace std;

namespace geometry
{
	using real = long double;
	const real EPS = 1e-9;
	bool EQ(real a,real b)
	{
		return abs(a - b) < EPS;
	}

	struct Point
	{
		real x,y;

		Point(real x_ = 0,real y_ = 0) : x(x_),y(y_) {}

		Point operator-() const
		{
			return Point(-x,-y);
		}

		Point operator+(const Point &rhs) const
		{
			return Point(x + rhs.x,y + rhs.y);
		}

		Point operator-(const Point &rhs) const
		{
			return Point(x - rhs.x,y - rhs.y);
		}

		Point operator*(const real k) const
		{
			return Point(x * k,y * k);
		}

		Point operator/(const real k) const
		{
			assert(!EQ(0,k));
			return Point(x / k,y / k);
		}

		bool operator<(const Point &rhs) const
		{
			return EQ(x,rhs.x) ? y < rhs.y : x < rhs.x;
		}

		bool operator==(const Point &rhs) const
		{
			return EQ(x,rhs.x) && EQ(y,rhs.y);
		}
	};

	istream &operator>>(istream &is,Point &p)
	{
		return is >> p.x >> p.y;
	}

	ostream &operator<<(ostream &os,const Point &p)
	{
		return os << p.x << " " << p.y;
	}

	struct Line
	{
		Point p1,p2;
		Line(Point p1_ = Point(),Point p2_ = Point()) : p1(p1_),p2(p2_) {}
	};

	struct Segment : Line
	{
		Segment(Point p1_ = Point(),Point p2_ = Point()) : Line(p1_,p2_) {}
	};

	struct Circle
	{
		Point O;
		real r;
		Circle(Point O_ = Point(),real r_ = 0) : O(O_),r(r_) {}
	};

	using Polygon = vector<Point>;

	Point vec(const Line &l)
	{
		return l.p2 - l.p1;
	}
	real norm2(const Point &p)
	{
		return p.x * p.x + p.y * p.y;
	}
	real abs(const Point &p)
	{
		return hypot(p.x,p.y);
	}

	real dot(const Point &a,const Point &b)
	{
		return a.x * b.x + a.y * b.y;
	}
	
	real cross(const Point &a,const Point &b)
	{
		return a.x * b.y - a.y * b.x;
	}

	Point rotate(const Point &p,const real &theta)
	{
		return Point(p.x * cos(theta) - p.y * sin(theta), p.x * sin(theta) + p.y * cos(theta));
	}

	Point rotate(const Point &a,const Point &p,const real &theta)
	{
		Point q = rotate(p - a,theta);
		return a + q;
	}

	enum
	{
		ONLINE_FRONT = -2,
		CLOCKWISE= -1,
		ON_SEGMENT = 0,
		COUNTER_CLOCKWISE = 1,
		ONLINE_BACK = 2
	};

	int ccw(const Point &a,const Point &b)
	{
		real C = cross(a,b);
		return C > EPS ? COUNTER_CLOCKWISE : C < -EPS ? CLOCKWISE : dot(a,b) < -EPS ? ONLINE_BACK : norm2(b) - norm2(a) > EPS ? ONLINE_FRONT : ON_SEGMENT;
	}

	int ccw(const Point &a,const Point &b,const Point &c)
	{
		return ccw(b - a,c - a);
	}

	bool orthogonal(const Point &a,const Point &b)
	{
		return EQ(dot(a,b),0);
	}
	
	bool orthogonal(const Line &a,const Line &b)
	{
		return orthogonal(vec(a),vec(b));
	}

	bool parallel(const Point &a,const Point &b)
	{
		return EQ(cross(a,b),0);
	}

	bool parallel(const Line &a,const Line &b)
	{
		return parallel(vec(a),vec(b));
	}

	bool intersect(const Line &l,const Point &p)
	{
		return parallel(vec(l),p - l.p1);
	}

	bool intersect(const Segment &s,const Point &p)
	{
		return ccw(s.p1,s.p2,p) == ON_SEGMENT;
	}

	bool intersect(const Segment &a,const Segment &b)
	{
		return ccw(a.p1,a.p2,b.p1) * ccw(a.p1,a.p2,b.p2) <= 0 && ccw(b.p1,b.p2,a.p1) * ccw(b.p1,b.p2,a.p2) <= 0;
	}

	Point cross_point(const Line &a,const Line &b)
	{
		real s1 = cross(vec(a),b.p1 - a.p1);
		real s2 = -cross(vec(a),b.p2 - a.p1);
		return b.p1 + vec(b) * (s1 / (s1 + s2));
	}

	Point crossPoint(const Line &s, const Line &t) {
		real d1 = cross(s.p2 - s.p1, t.p2 - t.p1);
		real d2 = cross(s.p2 - s.p1, s.p2 - t.p1);
		if(EQ(abs(d1), 0) && EQ(abs(d2), 0)) {
			return t.p1;
		}
		return t.p1 + (t.p2 - t.p1) * (d2 / d1);
	}


	enum
	{
		OUT,
		ON,
		IN
	};

	Polygon convex_hull(Polygon P,bool ONLINE = false,bool SORT = false)
	{
		if((int)P.size() <= 2)
		{
			return P;
		}
		sort(P.begin(),P.end());
		Polygon res(2 * P.size());
		int sz = 0;
		real threshold = EPS;
		if(ONLINE)
		{
			threshold = -EPS;
		}

		for(int i = 0;i < (int)P.size();i++)
		{
			while(sz >= 2 && cross(res[sz - 1] - res[sz - 2],P[i] - res[sz - 1]) < threshold)
			{
				sz--;
			}
			res[sz++] = P[i];
		}
		for(int i = (int)P.size() - 2,t = sz + 1;i >= 0;i--)
		{
			while(sz >= t && cross(res[sz - 1] - res[sz - 2],P[i] - res[sz - 1]) < threshold)
			{
				sz--;
			}
			res[sz++] = P[i];
		}
		res.resize(sz - 1);
		if(SORT)
		{
			int mi = 0;
			for(int i = 1;i < (int)res.size();i++)
			{
				if((EQ(res[mi].y,res[i].y) && res[mi].x > res[i].x) || res[mi].y > res[i].y)
				{
					mi = i;
				}
			}
			rotate(res.begin(),res.begin() + mi,res.end());
		}
		return res;
	}

	int convex_contain(const Polygon &P,const Point &p)
	{
		if(P[0] == p)
		{
			return ON;
		}

		int L = 0,R = (int)P.size();
		while(R - L > 1)
		{
			int M = (L + R) / 2;
			if(ccw(P[0],P[M],p) == CLOCKWISE)
			{
				R = M;
			}
			else
			{
				L = M;
			}
		}
		
		if(R == 1)
		{
			return OUT;
		}
		
		if(L + 1 == (int)P.size())
		{
			if(intersect(Segment(P[0],P[L]),p))
			{
				return ON;
			}
			return OUT;
		}

		if(L == 1)
		{
			if(intersect(Segment(P[0],P[L]),p))
			{
				return ON;
			}
		}

		real tri = cross(P[L] - p,P[R] - p);
		return EQ(tri,0) ? ON : tri < -EPS ? OUT : IN;
	}
}; //namespace geometry

#include<algorithm>
#include<vector>
#include<cassert>
using namespace std;
struct UnionFind
{
	private:
	int n;
	vector<int> par,siz;

	public:
	UnionFind(int n) :n(n),par(n,-1),siz(n,1) {}

	int root(int u) 
	{
		assert(0 <= u && u < n);
		return (par[u] < 0 ? u:par[u] = root(par[u]));
	}

	bool same(int u,int v)
	{
		assert(0 <= u && u < n && 0 <= v && v < n);
		return root(u) == root(v);
	}

	bool unite(int u,int v)
	{
		assert(0 <= u && u < n && 0 <= v && v < n);
		u = root(u),v = root(v);
		if(u == v) return false;

		if(siz[u] < siz[v]) swap(u,v);

		siz[u] += siz[v];
		par[v] = u;

		return true;
	}

	int size(int u)
	{
		assert(0 <= u && u < n);
		return siz[root(u)];
	}

	vector<vector<int>> components()
	{
		vector<vector<int>> ret(n);
		for(int u = 0;u < n;u++) ret[root(u)].push_back(u);

		ret.erase(remove_if(ret.begin(),ret.end(),[](vector<int> v) { return v.empty();}),ret.end());

		return ret;
	}
};
template<int m> struct modint
{
	private:
	unsigned int value;
	static constexpr int mod() {return m;}

	public:
	constexpr modint(const long long x = 0) noexcept
	{
		long long y = x;
		if(y < 0 || y >= mod())
		{
			y %= mod();
			if(y < 0) y += mod();
		}
		value = (unsigned int)y;
	}
	static constexpr int get_mod() noexcept {return m;}
	static constexpr int primitive_root() noexcept
	{
		assert(m == 998244353);
		return 3;
	}
	constexpr unsigned int val() noexcept {return value;}
	constexpr modint &operator+=(const modint &other) noexcept
	{
		value += other.value;
		if(value >= mod()) value -= mod();
		return *this;
	}
	constexpr modint &operator-=(const modint &other) noexcept
	{
		unsigned int x = value;
		if(x < other.value) x += mod();
		x -= other.value;
		value = x;
		return *this;
	}
	constexpr modint &operator*=(const modint &other) noexcept
	{
		unsigned long long x = value;
		x *= other.value;
		value = (unsigned int) (x % mod());
		return *this;
	}
	constexpr modint &operator/=(const modint &other) noexcept
	{
		return *this *= other.inverse();
	}
	constexpr modint inverse() const noexcept
	{
		assert(value);
		long long a = value,b = mod(),x = 1,y = 0;
		while(b)
		{
			long long q = a/b;
			a -= q*b; swap(a,b);
			x -= q*y; swap(x,y);
		}
		return modint(x);
	}
	constexpr modint power(long long N) const noexcept
	{
		assert(N >= 0);
		modint p = *this,ret = 1;
		while(N)
		{
			if(N & 1) ret *= p;
			p *= p;
			N >>= 1;
		}
		return ret;
	}
	constexpr modint operator+() {return *this;}
	constexpr modint operator-() {return modint() - *this;}
	constexpr modint &operator++(int) noexcept {return *this += 1;}
	constexpr modint &operator--(int) noexcept {return *this -= 1;}
	friend modint operator+(const modint& lhs, const modint& rhs) {return modint(lhs) += rhs;}
	friend modint operator-(const modint& lhs, const modint& rhs) {return modint(lhs) -= rhs;}
	friend modint operator*(const modint& lhs, const modint& rhs) {return modint(lhs) *= rhs;}
	friend modint operator/(const modint& lhs, const modint& rhs) {return modint(lhs) /= rhs;}
	friend ostream &operator<<(ostream &os,const modint &x) {return os << x.value;}
};
using mint = modint<998244353>;
/* using mint = modint<1000000007>; */

template<class S>
struct combination
{
	private:
	vector<S> f,invf;

	public:
	combination(int N = 0) : f(1,1),invf(1,1)
	{
		update(N);
	}

	void update(int N)
	{
		if((int)f.size() > N) return;
		int pi = (int)f.size();
		N = max(N,pi*2);

		f.resize(N+1),invf.resize(N+1);

		for(int i = pi;i <= N;i++) f[i] = f[i-1]*i;
		invf[N] = S(1)/f[N];
		for(int i = N-1;i >= pi;i--) invf[i] = invf[i+1]*(i+1);
	}

	S factorial(int N)
	{
		update(N);
		return f[N];
	}

	S invfactorial(int N)
	{
		update(N);
		return invf[N];
	}

	S P(int N,int K)
	{
		assert(0 <= K && K <= N);
		update(N);
		return f[N]*invf[N-K];
	}

	S C(int N,int K)
	{
		assert(0 <= K && K <= N);
		update(N);
		return f[N]*invf[K]*invf[N-K];
	}
};
combination<mint> C;

int ceil_log2(int n)
{
	int res = 0;
	while((1U << res) < (unsigned int)n)
	{
		res++;
	}
	return res;
}

template<class mint> void Butterfly(vector<mint> &a,bool inverse = false)
{
	int N = (int)a.size();
	int H = __builtin_ctz(N);
	assert(N == (1 << H));

	static constexpr int pr = mint::primitive_root();
	static bool first_call = true;
	static vector<mint> w(30),iw(30);
	if(first_call)
	{
		first_call = false;
		int cnt = __builtin_ctz(mint::get_mod() - 1);
		mint e = mint(pr).power((mint::get_mod() - 1) >> cnt);
		mint ie = e.inverse();
		for(int i = cnt;i >= 1;i--)
		{
			w[i] = e;
			iw[i] = ie;
			e *= e;
			ie *= ie;
		}
	}

	if(!inverse)
	{
		int width = N;
		int log = H;
		const mint im = w[2];
		while(width > 1)
		{
			mint cur = w[log];

			if(width == 2)
			{
				int offset = width >> 1;
				for(int i = 0;i < N;i += width)
				{
					mint root = 1;
					for(int j = i;j < i + offset;j++)
					{
						mint s = a[j],t = a[j + offset];
						a[j] = s + t;
						a[j + offset] = (s - t) * root;
						root *= cur;
					}
				}
				width >>= 1;
				log--;
			}
			else
			{
				int offset = width >> 2;
				for(int i = 0;i < N;i += width)
				{
					mint root = 1;
					for(int j = i;j < i + offset;j++)
					{
						mint root2 = root * root;
						mint root3 = root2 * root;
						mint s = a[j],t = a[j + offset],u = a[j + offset * 2],v = a[j + offset * 3];
						mint spu = s + u;
						mint smu = s - u;
						mint tpv = t + v;
						mint tmvim = (t - v) * im;
						a[j] = spu + tpv;
						a[j + offset] = (spu - tpv) * root2;
						a[j + offset * 2] = (smu + tmvim) * root;
						a[j + offset * 3] = (smu - tmvim) * root3;
						root *= cur;
					}
				}
				width >>= 2;
				log -= 2;
			}

		}
	}
	else
	{
		int width = H & 1 ? 2 : 4;
		int log = H & 1 ? 1 : 2;
		const mint im = iw[2];
		while(width <= N)
		{
			mint cur = iw[log];

			if(width == 2)
			{
				int offset = width >> 1;
				for(int i = 0;i < N;i += width)
				{
					mint root = 1;
					for(int j = i;j < i + offset;j++)
					{
						mint s = a[j],t = a[j + offset] * root;
						a[j] = s + t;
						a[j + offset] = s - t;
						root *= cur;
					}
				}
			}
			else
			{
				int offset = width >> 2;
				for(int i = 0;i < N;i += width)
				{
					mint root = 1;
					for(int j = i;j < i + offset;j++)
					{
						mint root2 = root * root;
						mint root3 = root2 * root;
						mint s = a[j],t = a[j + offset] * root2,u = a[j + offset * 2] * root,v = a[j + offset * 3] * root3;
						mint spt = s + t;
						mint smt = s - t;
						mint upv = u + v;
						mint umvim = (u - v) * im;
						a[j] = spt + upv;
						a[j + offset] = smt + umvim;
						a[j + offset * 2] = spt - upv;
						a[j + offset * 3] = smt - umvim;
						root *= cur;
					}
				}
			}

			width <<= 2;
			log += 2;
		}
	}
}

template<class mint> vector<mint> Convolution(vector<mint> a,vector<mint> b)
{
	int N = (int)a.size(),M = (int)b.size();

	if(min(N,M) <= 60)
	{
		vector<mint> res(N + M - 1);
		if(N < M)
		{
			swap(N,M);
			swap(a,b);
		}
		for(int i = 0;i < N;i++)
		{
			for(int j = 0;j < M;j++)
			{
				res[i + j] += a[i] * b[j];
			}
		}
		return res;
	}

	int L = 1 << ceil_log2(N + M - 1);

	a.resize(L);
	b.resize(L);

	Butterfly(a);
	Butterfly(b);

	for(int i = 0;i < L;i++)
	{
		a[i] *= b[i];
	}
	Butterfly(a,true);
	a.resize(N + M - 1);
	const mint invL = mint(L).inverse();
	for(int i = 0;i < N + M - 1;i++)
	{
		a[i] *= invL;
	}
	return a;
}

void Main()
{
	int N;
	cin >> N;
	vector<int> A(N);
	for(int i = 0;i < N;i++)
	{
		cin >> A[i];
	}

	vector<mint> dp(N);
	mint ans = 0;
	auto dfs = [&](auto dfs,int l,int r) -> void
	{
		if(l + 1 == r)
		{
			dp[l] += C.factorial(A[l]);
			ans += dp[l] * C.invfactorial(A[l]);
			return;
		}

		int mid = (l + r) / 2;
		dfs(dfs,l,mid);
		vector<mint> P(A[l] - A[mid - 1] + 1),Q(A[l] - A[r - 1] + 1);
		for(int i = l;i < mid;i++)
		{
			P[A[l] - A[i]] += dp[i];
		}
		for(int i = 0;i < (int)Q.size();i++)
		{
			Q[i] = C.invfactorial(i);
		}
		P = Convolution(P,Q);
		for(int i = mid;i < r;i++)
		{
			dp[i] += P[A[l] - A[i]];
		}
		dfs(dfs,mid,r);
	};
	dfs(dfs,0,N);
	cout << ans << "\n";
	/* vector<mint> dp(N); */
	/* mint ans = 0; */
	/* for(int i = N - 1;i >= 0;i--) */
	/* { */
	/* 	for(int j = i + 1;j < N;j++) */
	/* 	{ */
	/* 		dp[i] += dp[j] * C.C(A[i],A[j]); */
	/* 	} */
	/* 	dp[i]++; */
	/* 	ans += dp[i]; */
	/* 	cout << dp[i] << (i ? " ":"\n"); */
	/* } */
	/* cout << ans << "\n"; */
	/* const int M = (int)1e5; */
	/* for(int i = 0;i < N;i++) */
	/* { */
	/* 	vector<mint> dp(M + 1); */
	/* 	dp[A[i]] += C.factorial(A[i]); */
	/* 	for(int j = i + 1;j < N;j++) */
	/* 	{ */
	/* 		for(int k = A[j];k <= M;k++) */
	/* 		{ */
	/* 			dp[A[j]] += dp[k] * C.invfactorial(k - A[j]); */
	/* 		} */
	/* 	} */
	/* 	mint ans = 0; */
	/* 	for(int j = 0;j <= M;j++) */
	/* 	{ */
	/* 		ans += dp[j] * C.invfactorial(j); */
	/* 	} */
	/* 	cout << ans << (i + 1 == N ? "\n":" "); */
	/* } */
	/* vector<mint> dp(M + 1); */
	/* for(int i = 0;i < N;i++) */
	/* { */
	/* 	for(int j = A[i];j <= M;j++) */
	/* 	{ */
	/* 		dp[A[i]] += dp[j] * C.invfactorial(j - A[i]); */
	/* 	} */
	/* 	dp[A[i]] += C.factorial(A[i]); */
	/* } */
	/* mint ans = 0; */
	/* for(int i = 0;i <= M;i++) */
	/* { */
	/* 	ans += dp[i] * C.invfactorial(i); */
	/* } */
	/* cout << ans << "\n"; */
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int tt = 1;
	/* cin >> tt; */
	while(tt--) Main();
}
yukicoder

結果

ソースコード
