loj556 咱们去烧菜吧 polyExp+调和级数枚举

最新推荐文章于 2021-07-30 15:34:44 发布

BPMThor(BPM136)

最新推荐文章于 2021-07-30 15:34:44 发布

阅读量193

点赞数

分类专栏：生成函数多项式全家桶

本文链接：https://blog.csdn.net/qq_20139423/article/details/99405858

版权

生成函数同时被 2 个专栏收录

5 篇文章 0 订阅

订阅专栏

多项式全家桶

2 篇文章 0 订阅

订阅专栏

咋一看就是要算 $\frac{\prod_{i=1}^n (1-x^{a_i(b_i+1)})}{\prod_{i=1}^n (1-x^{a_i})}$
大力分治NTT？完蛋发现好像是 $O (n m l o g n l o g m)$ 的
看到一堆乘法，大力上个EXP和LN
发现变成了
$e^{\sum_{i=1}^n ln(1-x^{a_i(b_i+1)}) - \sum_{i=1}^n ln(1-x^{a_i})}$
仔细想想EGF里的东西，大力展开LN
就变成了
$Exp\{-\sum_{i=1}^n \sum_{j=1}^{ \infin } \frac{ x^{ a_i(b_i+1)j }}{j} + \sum_{i=1}^n \sum_{j=1}^{ \infin } \frac{ x^{a_i j} }{j} \}$
然后这东西我不会搞了
beginend大爷看了一眼，这不就先去个重，然后枚举j调和级数搞搞就好了吗
哦有道理x

#include <bits/stdc++.h>
#include <sys/timeb.h>

using namespace std;

#define SZ(x) ((int)(x).size())

int random(int l, int r) {
	static std::random_device rd;
	struct timeb timeSeed;
	ftime(&timeSeed);
	size_t seed = timeSeed.time * 1000 + timeSeed.millitm;  // milli time
	static std::mt19937 gen(seed);
	std::uniform_int_distribution<> u(l, r);
	return u(gen);
}

typedef long long ll;
typedef double db;
typedef long double ld;
typedef unsigned int ui;
typedef unsigned long long ull;
typedef pair<int, int> pii;


/eval多点求值,interpolation插值
//不可传入一个系数为负数的多项式，可能在减法的时候会爆mod
//常数稍微有点大
//G是TT的原根
//NTT前记得getRevRoot

int const TT = 998244353, G = 3;

inline int KSM (int a, int k) {
	int ret = 1 % TT;
	for (; k; k >>= 1, a = 1ll * a * a % TT)
		if (k & 1)
			ret = 1ll * ret * a % TT;
	return ret;
}

inline int add (int x, int y) {
    if (x + y >= TT)
        return x + y - TT;
    else
        return x + y;
}

inline int sub (int x, int y) {
    if (x - y >= 0)
        return x - y;
    else
        return x - y + TT;
}

namespace Polynom
{
	vector<int> rev, rt, e, one = {1};

	vector<int> operator + (vector<int> f, vector<int> g) {
		int n = (int)max(f.size(), g.size());
		f.resize(n);
		g.resize(n);
		for (int i = 0; i < n; ++i)
			f[i] = add(f[i], g[i]);
		return f;
	}

	vector<int> operator - (vector<int> f, vector<int> g) {
		int n = (int)max(f.size(), g.size());
		f.resize(n);
		g.resize(n);
		for (int i = 0; i < n; ++i)
			f[i] = sub(f[i], g[i]);
		return f;
	}

	void getRevRoot (int n) {
		int m = log(n) / log(2) + 1e-7;
		rev.resize(n);
		rt.resize(n);
		for (int i = 1; i < n; ++i)
			rev[i] = rev[i >> 1] >> 1 | (i & 1) << (m - 1);
		for (int len = 1, uni; len < n; len <<= 1) {
			uni = KSM(G, (TT ^ 1) / (len << 1));
			rt[len] = 1;
			for (int i = 1; i < len; ++i)
				rt[i + len] = 1ll * rt[i + len - 1] * uni % TT;
		}
	}

	void NTT (vector<int> &f, int n) {
		f.resize(n);
		for (int i = 0; i < n; ++i) 
			if (i < rev[i])
				swap(f[i], f[rev[i]]);
		for (int len = 1; len < n; len <<= 1)
			for (int i = 0; i < n; i += len << 1)
				for (int j = 0, x, y; j < len; j++) {
					x = f[i + j];
					y = 1ll * f[i + j + len] * rt[j + len] % TT;
					f[i + j] = add(x, y);
					f[i + j + len] = sub(x, y);
				}
	}

	vector<int> operator * (vector<int> f, vector<int> g) {
		int n = 1, m = (int)f.size() + (int)g.size() - 1, ivn;
		while (n < m)
			n <<= 1;
		ivn = KSM(n, TT - 2);
		getRevRoot(n);
		NTT(f, n);
		NTT(g, n);
		for (int i = 0; i < n; ++i)
			f[i] = 1ll * f[i] * g[i] % TT;
		reverse(f.begin() + 1, f.end());
		NTT(f, n);
		f.resize(m);
		for (int i = 0; i < m; ++i)
			f[i] = 1ll * f[i] * ivn % TT;
		return f;
	}

	vector<int> polyInv (vector<int> f, int m) {
		if (m == 1)
			return { KSM(f[0], TT - 2) };
		f.resize(m);
		vector<int> g = polyInv(f, m >> 1), h;
		g.resize(m);
		h.resize(m);
		for (int i = (m >> 1) - 1; ~i; --i)
			h[i] = (g[i] << 1) % TT;
		int n = m << 1, ivn = KSM(n, TT - 2);
		getRevRoot(n);
		NTT(f, n);
		NTT(g, n);
		for (int i = 0; i < n; ++i)
			f[i] = 1ll * f[i] * g[i] % TT * g[i] % TT;
		reverse(f.begin() + 1, f.end());
		NTT(f, n);
		for (int i = 0; i < m; ++i)
			h[i] = sub(h[i], 1ll * f[i] * ivn % TT);
		return h;
	}

	vector<int> operator ~ (vector<int> f) {
		int n = 1, m = (int)f.size();
		while (n < m)
			n <<= 1;
		f = polyInv(f, n);
		f.resize(m);
		return f;
	}

	vector<int> polyDeri (vector<int> f) {
		if (f.empty())
			return f;
		int n = (int)f.size();
		for (int i = 1; i < n; ++i)
			f[i - 1] = 1ll * f[i] * i % TT;
		f.resize(n - 1);
		return f;
	}

	vector<int> polyInte (vector<int> f) {
		int n = (int)f.size();
		f.resize(n + 1);
		for (int i = n; ~i; --i)
			f[i] = 1ll * f[i - 1] * KSM(i, TT - 2) % TT;
		return f;
	}

	vector<int> polyLn (vector<int> f) {
		int n = (int)f.size();
		f = polyInte((~f) * polyDeri(f));
		f.resize(n);
		return f;
	}

	vector<int> polyExp (vector<int> f, int n) {
		if (n == 1)
			return one;
		f.resize(n);
		vector<int> g = polyExp(f, n >> 1);
		g.resize(n);
		return g * (one - polyLn(g) + f);
	}

	vector<int> polyExp (vector<int> f) {
		int n = 1, m = (int)f.size();
		while (n < m)
			n <<= 1;
		f = polyExp(f, n);
		f.resize(m);
		return f;
	}

	vector<int> operator ^ (vector<int> f, int x) {
		int n = 0, m = (int)f.size();
		while (n < m - 1 && !f[n])
			n++;
		f.erase(f.begin(), f.begin() + n);
		e = {KSM(f.empty() ? 0 : f[0], x)};
		f = polyLn(f);
		f.resize(n = max(0ll, m - 1ll * n * x));
		for (int i = 1; i < n; ++i)
			f[i] = 1ll * f[i] * x % TT;
		f = polyExp(f);
		reverse(f.begin(), f.end());
		for (int i = m - n; i; --i)
			f.push_back(0);
		reverse(f.begin(), f.end());
		f.resize(m);
		return f;
	}

	vector<int> operator / (vector<int> f, vector<int> g) {
		int n = f.size(), m = g.size();
		reverse(f.begin(), f.end());
		reverse(g.begin(), g.end());
		f.resize(n - m + 1);
		g.resize(n - m + 1);
		vector<int> ret = f * (~g);
		ret.resize(n - m + 1);
		reverse(ret.begin(), ret.end());
		return ret;
	}

	vector<int> operator % (vector<int> f, vector<int> g) {
		vector<int> ret = f / g;
		ret = f - (g * ret);
		ret.resize(g.size() - 1);
		return ret;
	}
} // namespace Polynom

using namespace Polynom;

vector<int> a_mul[70000 << 2];   

void getMula (vector<int> const& a, int k, int l, int r) {
	if (l == r) {
		a_mul[k] = { -a[l] + TT, 1 };
		return;
	}
	int mid = (l + r) >> 1;
	getMula(a, k << 1, l, mid);
	getMula(a, k << 1 | 1, mid + 1, r);
	a_mul[k] = a_mul[k << 1] * a_mul[k << 1 | 1];
}

void eval_solve (vector<int> f, int k, int l, int r, vector<int>& ans) {
	if (l == r) {
		ans[l] = f[0];
		return;
	}
	int mid = (l + r) >> 1;
	eval_solve(f % a_mul[k << 1], k << 1, l, mid, ans);
	eval_solve(f % a_mul[k << 1 | 1], k << 1 | 1, mid + 1, r, ans);
}

vector<int> eval (vector<int> f, vector<int> const& a) {
	int m = a.size();
	getMula(a, 1, 0, m - 1);
	vector<int> ret;
	ret.resize(m);
	eval_solve(f, 1, 0, m - 1, ret);
	return ret;
}

vector<int> interpolation_solve (vector<int> const& pr_x, vector<int> const& y, int k, int l, int r) {
	if (l == r) 
		return { int(1LL * y[l] * KSM(pr_x[l], TT - 2) % TT) };
	int mid = (l + r) >> 1;
	vector<int> f1 = interpolation_solve(pr_x, y, k << 1, l, mid);
	if (f1.size() > y.size())
		f1.resize(y.size());
	vector<int> f2 = interpolation_solve(pr_x, y, k << 1 | 1, mid + 1, r);
	if (f2.size() > y.size())
		f2.resize(y.size());
	return f1 * a_mul[k << 1 | 1] + f2 * a_mul[k << 1];
}

vector<int> interpolation (vector<int> const& x, vector<int> const& y) {
	assert(x.size() == y.size() && x.size());
	int n = x.size();
	getMula(x, 1, 0, n - 1);
	vector<int> prod = a_mul[1];
	prod = polyDeri(prod);
	vector<int> prod_x;
	prod_x.resize(n);
	eval_solve(prod, 1, 0, n - 1, prod_x);
	return interpolation_solve(prod_x, y, 1, 0, n - 1);
}

int const N = 100005;

int inv[N];

void solve (vector<int> const& a, vector<int>& f, int flag, int n) {
	auto mp = map<int, int>();
	for (auto& v : a)
		mp[v]++;
	for (auto const& p : mp) 
		for (ll j = 1; j * p.first <= n; ++j) {
			ll now = j * p.first;
			f[now] += 1LL * flag * p.second * inv[j] % TT + (flag == -1 ? TT : 0);
			if (f[now] >= TT)
				f[now] -= TT;
		}
}

int main() {
	freopen("a.in", "r", stdin);
	freopen("a.out", "w", stdout);
	int n, m;
	cin >> n >> m;

	inv[0] = inv[1] = 1;
	for (int i = 2; i <= n; ++i)
		inv[i] = 1LL * (TT - TT / i) * inv[TT % i] % TT;

	auto a = vector<int>(), b = vector<int>();
	ll anss = 1;
	for (int i = 0; i < m; ++i) {
		int x, y;
		cin >> x >> y;
		if (x == 0) {
			anss = anss * (y + 1) % TT;
			continue;
		}
		if (y) {
			ll t1 = 1LL * x * (y + 1);
			if (t1 <= n)
				a.push_back(t1);
		}
		b.push_back(x);
	}
	auto f = vector<int>(n + 1, 0);
	solve(a, f, -1, n);
	solve(b, f, 1, n);
	auto ans = polyExp(f);
	for (int i = 1; i <= n; ++i)
		cout << ans[i] << '\n';
    return 0;
}

BPMThor(BPM136)

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
loj556 咱们去烧菜吧 polyExp+调和级数枚举

咋一看就是要算∏i=1n(1−xai(bi+1))∏i=1n(1−xai)\frac{\prod_{i=1}^n (1-x^{a_i(b_i+1)})}{\prod_{i=1}^n (1-x^{a_i})}∏i=1n(1−xai)∏i=1n(1−xai(bi+1))大力分治NTT？完蛋发现好像是O(nmlognlogm)O(nmlognlogm)O(nmlognlogm)的看到一堆...
复制链接

扫一扫

专栏目录