专栏文章

一种不太快速计算多个少项式集合幂级数异或卷积的方法

算法·理论参与者 3已保存评论 3

文章操作

快速查看文章及其快照的属性,并进行相关操作。

当前评论
3 条
当前快照
1 份
快照标识符
@mjbozefv
此快照首次捕获于
2025/12/19 01:06
2 个月前
此快照最后确认于
2026/02/19 01:28
16 小时前
查看原文
首先求 (aix+bixSi){\prod}_\oplus(a_ix^{\emptyset}+b_ix^{S_i}),其中二项式个数为 nnSi{1,2,,m}S_i\subseteq \{1,2,\cdots,m\}
考虑 FWT,要将 FWT(aix+bixSi)\operatorname{FWT}(a_ix^{\emptyset}+b_ix^{S_i}) 对位乘得到 g(x)g(x) 之后 IFWT。
[xT]g(x)=[xT]FWT(aix+bixSi)=(ai+(1)SiTbi)=exp(ln(ai+(1)SiTbi))=exp([xT]FWT(ln(ai+bi)+ln(aibi)2x+ln(ai+bi)ln(aibi)2xSi))=exp([xT]FWT((ln(ai+bi)+ln(aibi)2x+ln(ai+bi)ln(aibi)2xSi)))\begin{aligned} & [x^T]g(x) \\ = & \prod [x^T]\operatorname{FWT}(a_i x^{\emptyset}+b_ix^{S_i}) \\ = & \prod (a_i+(-1)^{|S_i\cap T|}b_i) \\ = & \exp(\sum \ln (a_i+(-1)^{|S_i\cap T|}b_i)) \\ = & \exp(\sum [x^T] \operatorname{FWT}(\dfrac {\ln(a_i+b_i)+\ln(a_i-b_i)}{2} x^{\emptyset} + \dfrac {\ln(a_i+b_i)-\ln(a_i-b_i)}{2} x^{S_i})) \\ = & \exp([x^T]\operatorname{FWT}(\sum(\dfrac {\ln(a_i+b_i)+\ln(a_i-b_i)}{2} x^{\emptyset} + \dfrac {\ln(a_i+b_i)-\ln(a_i-b_i)}{2} x^{S_i})))\end{aligned}
只要对数能求就能直接做。
如果有较多项呢?对于 ai,jxSi,j\sum a_{i,j}x^{S_{i,j}},对于每个 ii,找一组 F2mF_2^m 中的基底,它们张成一个秩为 2ki2^{k_i} 的子空间。把 Si,jS_{i,j} 变换到 Si,jS'_{i,j},在子空间里做 FWT,取对数,做 IFWT,再把点值加到原空间的对应位置(具体来说,对于一个子空间的向量 vv,取出所有 vj=1v_j=1Si,jS_{i,j},将其 xor\operatorname{xor} 起来),做 FWT,取指数,做 IFWT。
复杂度 O(m2m+iki2ki+L(i2ki))O(m2^m+\sum_i k_i2^{k_i}+L(\sum_i 2^{k_i})),其中 L(x)L(x) 是计算 xx 个数的对数的复杂度,kik_i 是所有 Si,jS_{i,j} 张成子空间的大小。
不能取对数的时候,考虑额外记录 00 的指数。对于 x0x\neq 0,记 logx=(lnx,0)\log x=(\ln x,0),记 ln0=(0,1)\ln 0=(0,1) 即可。做完原空间上的 FWT 后,若某个向量是 (x,0)(x,0),将其替换为数值 exp(x)\exp(x),否则替换为数值 00,再 IFWT 即可。
以下是示例代码。由于没有高消,其累加复杂度是 O(2mi)O(2^{m_i})mim_i 是第 ii 个少项式的项数。可以在全局常量里调 mm 和模数,也可以改成输入模数。使用了较为朴素的方法求原根,在示例代码中这可能成为复杂度瓶颈。
CPP
#include<iostream>
#include<cstdio>
#include<fstream>
#include<queue>
#include<algorithm>
#include<cstring>
#include<random>
#include<map>
#include<cassert>
#include<stack>
#define fopen(x, y) freopen(x".in", "r", stdin); freopen(y".out", "w", stdout);
#define int long long
#define pii pair<int, int>
#define fi first
#define se second
#ifdef int
#define inf 0x3f3f3f3f3f3f3f3fll
#else
#define inf 0x3f3f3f3f
#endif

#define maxB 65536
#define maxm 1048576 

using namespace std;

//int m, M, mod, g, b, B, Mod;
const int m = 20, M = 1ll << m, mod = 998244353, g = 3, b = 16, B = 65536, Mod = 998244352ll * M;

inline int qpow(int x, int y) { int xum = 1; while(y) { if(y & 1) (xum *= x) %= mod; (x *= x) %= mod; y >>= 1; } return xum; }

int n, a, s[maxm], x[105], w[105], hfi[maxm], hse[maxm], ffi[maxm], fse[maxm];
int l3, l_1;

//void init_0();
void init_1();
int gxp(int x);
int gln(int x);

void fwt(int *f, int m, const int mod) {
	const int M = 1ll << m;
	for(int j = 0; j < m; j++) {
		int L = 1ll << (j + 1), bL = 1ll << j;
		for(int t = 0; t < M; t += L) {
			int *x = &f[t], *y = x + bL;
			for(int i = 0; i < bL; i++, x++, y++) {
				(*x += *y) %= mod;
				*y = (*x - (*y << 1)) % mod;
			}
		}
	}
}

void work() {
//	mod = 998244353; init_0();
	
	init_1();
	
	cin >> n;
	
	for(int i = 1; i <= n; i++) {
		int k;
		cin >> k;
		const int K = 1ll << k;
		for(int j = 0; j < K; j++) hfi[j] = hse[j] = 0;
		for(int j = 0; j < k; j++) {
			cin >> x[j] >> w[j];
			hfi[1ll << j] = w[j];
		}
		s[0] = 0;
		for(int j = 1; j < K; j++) {
			s[j] = s[j ^ (j & -j)] ^ x[__builtin_ctzll(j & -j)];
		}
		fwt(hfi, k, mod);
		for(int j = 0; j < K; j++) {
			(hfi[j] += mod) %= mod;
			if(hfi[j] == 0) hse[j] = 1ll << (m - k);
			else hfi[j] = (gln(hfi[j]) << (m - k)) % Mod;
		}
		fwt(hfi, k, Mod);
		fwt(hse, k, Mod);
		for(int j = 0; j < K; j++) {
			(ffi[s[j]] += hfi[j]) %= Mod;
			(fse[s[j]] += hse[j]) %= Mod;
		}
	}
	
	fwt(ffi, m, Mod);
	fwt(fse, m, Mod);
	
	int sum = 0, iv = qpow(M, mod - 2);
	for(int i = 0; i < M; i++) {
		if(fse[i]) {
			ffi[i] = 0;
			continue;
		}
		(ffi[i] >>= m) %= (mod - 1);
		if(ffi[i] < 0) ffi[i] += mod - 1;
		ffi[i] = gxp(ffi[i]);
	}
	fwt(ffi, m, mod);
	for(int i = 0; i < M; i++) {
		(ffi[i] *= iv) %= mod;
		cout << (ffi[i] + mod) % mod << ' ';
	}
	cout << '\n';
}

signed main() {
	ios::sync_with_stdio(false); cin.tie(0);
	int _ = 1;
	while(_--) work();
}

/*
void init_0() {
	Mod = (mod - 1) * maxn;
	vector<int> t;
	int s = mod - 1;
	B = 1, b = 0;
	while(B * B < mod) B <<= 1, b++;
	for(int i = 2; i * i <= s; i++) {
		if(s % i == 0) {
			while(s % i == 0) s /= i;
			t.push_back(i);
		}
	}
	for(int i = 2; i < mod; i++) {
		int flg = 1;
		for(int d : t) {
			if(qpow(i, (mod - 1) / d) == 1) {
				flg = 0;
				break;
			}
		}
		if(flg) {
			g = i;
			break;
		}
	}
}
//*/

int gb, gk[maxB], gbk[maxB];
pii gt[maxB + 1];
int find(int x) {
	pii *p = lower_bound(gt, gt + B, make_pair(x, 0ll));
	if(p->fi == x) return p->se;
	return -1;
}
int gln(int x) {
	int cnt = 0, p = find(x);
	while(p == -1) {
		(x *= gb) %= mod;
		cnt++;
		p = find(x);
	}
	return ((p - cnt * B) % (mod - 1) + (mod - 1)) % (mod - 1);
}
void init_1() {
	gb = 1;
	gk[0] = 1;
	for(int i = 1; i < B; i++) gk[i] = gk[i - 1] * g % mod;
	gb = gk[B - 1] * g % mod;
	gbk[0] = 1;
	for(int i = 1; i < B; i++) gbk[i] = gbk[i - 1] * gb % mod;
	
	for(int i = 0; i < B; i++) gt[i] = make_pair(gk[i], i);
	sort(gt, gt + B); 
	gt[B] = make_pair(mod, -1);
}
int gxp(int x) {
	return gk[x & (B - 1)] * gbk[x >> b] % mod;
}

评论

3 条评论,欢迎与作者交流。

正在加载评论...