【UR #5】怎样跑得更快莫比乌斯反演

最新推荐文章于 2022-01-18 23:39:00 发布

原创最新推荐文章于 2022-01-18 23:39:00 发布 · 394 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#莫比乌斯反演

xgc的做题记录同时被 2 个专栏收录

168 篇文章

订阅专栏

莫比乌斯反演

6 篇文章

订阅专栏

本文介绍了一种利用莫比乌斯反演和快速幂模运算解决特定形式线性方程组的方法，通过预处理得到莫比乌斯函数和逆元，将原方程转换并求解未知数序列。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description
给你一个b序列，求x序列，b和x满足下式。
$∑j=1ngcd(i,j)c⋅lcm(i,j)d⋅xj≡bi(modp)\sum_{j = 1}^{n} gcd(i, j)^c \cdot lcm(i, j)^d \cdot x_j \equiv b_i \pmod{p}$

Sample Input
3 1 0 2
1 0 0
1 2 3

Sample Output
499122179 998244352 499122176
998244352 1 1

式子首先可以化成：
$∑j=1njdgcd(i,j)c−d⋅xj≡biid(modp)\sum_{j = 1}^{n} j^dgcd(i, j)^{c-d} \cdot x_j \equiv \frac {b_i}{i^d} \pmod{p}$
我们考虑设 $g_i$ 为 $i^d$ ， $f_i$ 为 $i^{c-d}$ ， $w_i$ 为 $biid\frac {b_i}{i^d}$ 。
$∑j=1ngj⋅fgcd(i,j)⋅xj≡wi(modp)\sum_{j = 1}^{n} g_j \cdot f_{gcd(i, j)} \cdot x_j \equiv w_i \pmod{p}$
考虑莫反搞掉这个 $g c d$ ，设 $fri=∑d∣ifd⋅μidfr_i=\sum_{d|i}f_d \cdot \mu_{\frac id}$ 。
设 $ti=xi⋅git_i=x_i \cdot g_i$ ，那么式子变为：
$∑j=1n∑d∣i,d∣jfrd⋅tj≡wi(modp)\sum_{j = 1}^{n}\sum_{d|i,d|j} fr_d \cdot t_j \equiv w_i \pmod{p}$
$∑d∣ifrd∑d∣jtj≡wi(modp)\sum_{d|i} fr_d \sum_{d|j} t_j \equiv w_i \pmod{p}$
设 $qd=∑d∣jtjq_d=\sum_{d|j} t_j$ ，
$∑d∣ifrd⋅qd≡wi(modp)\sum_{d|i} fr_d\cdot q_d \equiv w_i \pmod{p}$
设 $Fd=frd⋅qdF_d=fr_d\cdot q_d$ ，
莫反一下，这样就搞得出 $F$ ，进而搞得出 $q$ 。
对于这个 $q$ 再莫反一下，就搞了出 $t$ 。

#include <cstdio>
#include <cstring>

using namespace std;
typedef long long LL;
const int N = 300001;
const LL mod = 998244353;
int read() {
	int s = 0, f = 1; char ch = getchar();
	while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
	return s * f;
}
void put(int x) {
	if(x >= 10) put(x / 10);
	putchar(x % 10 + '0');
}

LL w[N], g[N], f[N], fr[N], F[N], q[N], t[N], x[N];
int plen, p[N], mu[N];
bool v[N];

void get_p() {
	mu[1] = 1;
	for(int i = 2; i < N; i++) {
		if(!v[i]) p[++plen] = i, mu[i] = -1;
		for(int j = 1; j <= plen && (LL)p[j] * i < N; j++) {
			v[i * p[j]] = 1;
			if(i % p[j] == 0) {mu[i * p[j]] = 0; break;}
			mu[i * p[j]] = -mu[i];
		}
	}
}

LL pow_mod(LL a, int k) {
	LL ans = 1;
	k = (k + mod - 1) % (mod - 1);
	while(k) {
		if(k & 1) (ans *= a) %= mod;
		(a *= a) %= mod; k /= 2;
	} return ans;
}

int main() {
	int n = read(), c = read(), d = read(), tt = read();
	get_p();
	for(int i = 1; i <= n; i++) g[i] = pow_mod(i, d), f[i] = pow_mod(i, c - d);
	for(int i = 1; i <= n; i++) {
		for(int j = i, s = 1; j <= n; j += i, s++) {
			(fr[j] += mu[s] * f[i]) %= mod;
		} (fr[i] += mod) %= mod;
	}
	while(tt--) {
		bool bk = 0;
		memset(F, 0, sizeof(F));
		memset(t, 0, sizeof(t));
		for(int i = 1; i <= n; i++) w[i] = read(), w[i] = w[i] * pow_mod(g[i], mod - 2) % mod;
		for(int i = 1; i <= n; i++) {
			for(int j = i, s = 1; j <= n; j += i, s++) {
				(F[j] += mu[s] * w[i]) %= mod;
			} (F[i] += mod) %= mod;
		} for(int i = 1; i <= n; i++) {
			if(F[i] != 0 && fr[i] == 0) {bk = 1; break;}
			q[i] = F[i] * pow_mod(fr[i], mod - 2) % mod;
		} if(bk) {puts("-1"); continue;}
		for(int i = 1; i <= n; i++) {
			for(int j = i, s = 1; j <= n; j += i, s++) {
				(t[i] += mu[s] * q[j]) %= mod;
			} (t[i] += mod) %= mod;
		} for(int i = 1; i <= n; i++) {
			if(t[i] != 0 && g[i] == 0) {bk = 1; break;}
			x[i] = t[i] * pow_mod(g[i], mod - 2) % mod;
		} if(bk) {puts("-1"); continue;}
		for(int i = 1; i <= n; i++) put(x[i]), putchar(' ');
		putchar('\n');
	}
	return 0;
}