「算法笔记」莫比乌斯反演进阶

YY 的 GCD

题目大意

求下列式子的值:

∑ i = 1 n ∑ j = 1 m [ gcd ⁡ { i , j } ∈ P ]   ( P  为素数集合 ) \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd\{i, j\} \in P]\ (P\ \text{为素数集合}) i=1nj=1m[gcd{i,j}P] (P 为素数集合)

1 0 4 10^4 104 组询问, n , m ≤ 1 0 7 n, m \le 10^7 n,m107

思路分析

将原式变形:


∑ p ∈ P ∑ i = 1 ⌊ n p ⌋ ∑ j = 1 ⌊ m p ⌋ [ gcd ⁡ { i , j } = 1 ] \sum_{p \in P} \sum_{i = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{p} \rfloor} [\gcd\{i, j\} = 1] pPi=1pnj=1pm[gcd{i,j}=1]

∑ p ∈ P ∑ i = 1 ⌊ n p ⌋ ∑ j = 1 ⌊ m p ⌋ ∑ d ∣ i , d ∣ j μ ( d ) \sum_{p \in P} \sum_{i = 1}^{\lfloor \frac{n}{p} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{p} \rfloor} \sum_{d | i, d | j} \mu(d) pPi=1pnj=1pmdi,djμ(d)

∑ p ∈ P ∑ d = 1 min ⁡ { ⌊ n p ⌋ , ⌊ n p ⌋ } μ ( d ) ∑ i = 1 ⌊ n p d ⌋ ∑ j = 1 ⌊ m p d ⌋ 1 \sum_{p \in P} \sum_{d = 1}^{\min\{\lfloor \frac{n}{p} \rfloor, \lfloor \frac{n}{p} \rfloor\}} \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{pd} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{pd} \rfloor} 1 pPd=1min{pn,pn}μ(d)i=1pdnj=1pdm1

∑ p ∈ P ∑ d = 1 min ⁡ { ⌊ n p ⌋ , ⌊ n p ⌋ } μ ( d ) ⌊ n p d ⌋ ⌊ m p d ⌋ \sum_{p \in P} \sum_{d = 1}^{\min\{\lfloor \frac{n}{p} \rfloor, \lfloor \frac{n}{p} \rfloor\}} \mu(d) \lfloor \frac{n}{pd} \rfloor \lfloor \frac{m}{pd} \rfloor pPd=1min{pn,pn}μ(d)pdnpdm

∑ T = 1 min ⁡ { n , m } ∑ p ∣ T , p ∈ P μ ( T p ) ⌊ n T ⌋ ⌊ m T ⌋ \sum_{T = 1}^{\min\{n, m\}} \sum_{p | T, p \in P} \mu(\frac{T}{p}) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor T=1min{n,m}pT,pPμ(pT)TnTm

∑ T = 1 min ⁡ { n , m } ⌊ n T ⌋ ⌊ m T ⌋ ∑ p ∣ T , p ∈ P μ ( T p ) \sum_{T = 1}^{\min\{n, m\}} \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor \color{red}{\sum_{p | T, p \in P} \mu(\frac{T} {p})} T=1min{n,m}TnTmpT,pPμ(pT)


红色部分可以使用线性筛算出 μ \mu μ,然后枚举每个素数并更新它的倍数算出,时间复杂度 O ( n log ⁡ log ⁡ n ) O(n \log \log n) O(nloglogn)。黑色部分可以使用数论分块的技巧求出。

代码实现

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long llong;
const int maxn = 1e7;
int T, n, m, p[maxn + 3];
int mrk[maxn + 3], miu[maxn + 3], sum[maxn + 3];

int main() {
	miu[1] = 1;
	for (int i = 2; i <= maxn; i++) {
		if (!mrk[i]) p[++m] = i, miu[i] = -1;
		for (int j = 1, t; j <= m && (t = i * p[j]) <= maxn; j++) {
			mrk[t] = 1;
			if (i % p[j] == 0) break;
			miu[t] = -miu[i];
		}
	}
	for (int i = 1; i <= m; i++) {
		for (int j = 1, t; (t = p[i] * j) <= maxn; j++) {
			sum[t] += miu[j];
		}
	}
	for (int i = 1; i <= maxn; i++) {
		sum[i] += sum[i - 1];
	}
	scanf("%d", &T);
	while (T--) {
		scanf("%d %d", &n, &m);
		if (n > m) swap(n, m);
		llong ans = 0;
		for (int i = 1, j; i <= n; i = j + 1) {
			j = min(n / (n / i), m / (m / i));
			ans += 1ll * (n / i) * (m / i) * (sum[j] - sum[i - 1]);
		}
		printf("%lld\n", ans);
	}
	return 0;
}

NOI 2016 循环之美

题目大意

求下列式子的值:

∑ i = 1 n ∑ j = 1 m [ gcd ⁡ { i , j } = 1 ] [ gcd ⁡ { j , k } = 1 ] \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd\{i, j\} = 1] [\gcd\{j, k\} = 1] i=1nj=1m[gcd{i,j}=1][gcd{j,k}=1]

n , m ≤ 1 0 9 , k ≤ 2000 n, m \le 10^9, k \le 2000 n,m109,k2000

思路分析

将原式变形:


∑ i = 1 m [ gcd ⁡ { i , k } = 1 ] ∑ j = 1 n [ gcd ⁡ { i , j } = 1 ] \sum_{i = 1}^{m} [\gcd\{i, k\} = 1] \sum_{j = 1}^{n} [\gcd\{i, j\} = 1] i=1m[gcd{i,k}=1]j=1n[gcd{i,j}=1]

∑ i = 1 m [ gcd ⁡ { i , k } = 1 ] ∑ j = 1 n ∑ d ∣ i , j μ ( d ) \sum_{i = 1}^{m} [\gcd\{i, k\} = 1] \sum_{j = 1}^{n} \sum_{d | i, j} \mu(d) i=1m[gcd{i,k}=1]j=1ndi,jμ(d)

∑ d = 1 min ⁡ { n , m } μ ( d ) ∑ i = 1 ⌊ m d ⌋ [ gcd ⁡ { i ⋅ d , k } = 1 ] ∑ j = 1 ⌊ n d ⌋ 1 \sum_{d = 1}^{\min\{n, m\}} \mu(d) \sum_{i = 1}^{\lfloor \frac{m}{d} \rfloor}[\gcd\{i \cdot d, k\} = 1] \sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor} 1 d=1min{n,m}μ(d)i=1dm[gcd{id,k}=1]j=1dn1

∑ d = 1 min ⁡ { n , m } μ ( d ) ⋅ [ gcd ⁡ { d , k } = 1 ] ⋅ ⌊ n d ⌋ ∑ i = 1 ⌊ m d ⌋ [ gcd ⁡ { i , k } = 1 ] \sum_{d = 1}^{\min\{n, m\}} \color{blue}{\mu(d) \cdot [\gcd\{d, k\} = 1]} \cdot \lfloor \frac{n}{d} \rfloor\sum_{i = 1}^{\lfloor \frac{m}{d} \rfloor}\color{red}{[\gcd\{i, k\} = 1]} d=1min{n,m}μ(d)[gcd{d,k}=1]dni=1dm[gcd{i,k}=1]


现在瓶颈就变成了计算红色部分和蓝色部分的前缀和。

注意到 k ≤ 2000 k \le 2000 k2000,于是红色部分直接可以递推计算 ≤ k \le k k 的部分,然后在乘上计数的次数即可。

对于蓝色部分,令 S ( n , k ) = ∑ i = 1 n [ gcd ⁡ ( i , k ) = 1 ] ⋅ μ ( i ) S(n, k) = \sum_{i = 1}^{n} [\gcd(i, k) = 1] \cdot \mu(i) S(n,k)=i=1n[gcd(i,k)=1]μ(i)。对原式进行变形,得到原式 = ∑ d ∣ k μ ( d ) ∑ i = 1 ⌊ n d ⌋ μ ( i ⋅ d ) = \sum_{d | k} \mu(d) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(i \cdot d) =dkμ(d)i=1dnμ(id)。发现如果 gcd ⁡ ( i , d ) ≠ 1 \gcd(i, d) \neq 1 gcd(i,d)̸=1,那么 μ ( i ⋅ d ) = 0 \mu(i \cdot d) = 0 μ(id)=0,对答案没有影响。于是可以继续化简,得到 S ( n , k ) = ∑ d ∣ k μ ( d ) 2 S ( ⌊ n d ⌋ , d ) S(n, k) = \sum_{d | k} \mu(d) ^ 2 S(\lfloor \frac{n}{d} \rfloor, d) S(n,k)=dkμ(d)2S(dn,d)。直接记忆化搜索 + + + 杜教筛即可。

代码实现

#include <cstdio>
#include <map>
#include <algorithm>
using namespace std;

typedef pair<int, int> pii;
const int maxm = 1e7, maxk = 2000;
bool mrk[maxm + 3];
int n, m, k, cnt, p[maxm + 3], miu[maxm + 3], sum[maxm + 3], f[maxk + 3];
map<pii, int> mp;

int gcd(int a, int b) {
	return b ? gcd(b, a % b) : a;
}

int F(int n) {
	return f[n % k] + n / k * f[k];
}

int S(int n, int k) {
	if ((k == 1 && n <= maxm) || !n) {
		return sum[n];
	}
	pii x = pii(n, k);
	if (mp.find(x) != mp.end()) {
		return mp[x];
	}
	int res;
	if (k == 1) {
		res = 1;
		for (int l = 2, r; l <= n; l = r + 1) {
			r = n / (n / l);
			res -= (r - l + 1) * S(n / l, 1);
		}
	} else {
		res = 0;
		for (int i = 1, j; i * i <= k; i++) {
			if (k % i == 0) {
				j = k / i;
				if (miu[i]) res += S(n / i, i);
				if (miu[j] && i * i != k) res += S(n / j, j);
			}
		}
	}
	return mp[x] = res;
}

int main() {
	#ifdef LOCAL
		freopen("input.md", "r", stdin);
		freopen("output.md", "w", stdout);
	#endif
	miu[1] = 1;
	for (int i = 2; i <= maxm; i++) {
		if (!mrk[i]) {
			p[++cnt] = i;
			miu[i] = -1;
		}
		for (int j = 1, t; j <= cnt && (t = i * p[j]) <= maxm; j++) {
			mrk[t] = 1;
			if (i % p[j] == 0) break;
			miu[t] = -miu[i];
		}
	}
	for (int i = 1; i <= maxm; i++) {
		sum[i] = sum[i - 1] + miu[i];
	}
	scanf("%d %d %d", &n, &m, &k);
	for (int i = 1; i <= k; i++) {
		f[i] = f[i - 1] + (gcd(i, k) == 1);
	}
	int x = min(n, m);
	long long ans = 0;
	for (int l = 1, r; l <= x; l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		ans += 1ll * (n / l) * F(m / l) * (S(r, k) - S(l - 1, k));
	}
	printf("%lld\n", ans);
	return 0;
}

UOJ 62 怎样跑的更快

题目大意

神仙题,参考官方题解

思路分析

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

代码实现

#include <cstdio>

const int maxn = 1e5, mod = 998244353;
int n, c, d, q, f_r[maxn + 3], f_z[maxn + 3], hx[maxn + 3];

int qpow(int a, int b) {
	int c = 1;
	for (; b; b >>= 1, a = 1ll * a * a % mod) {
		if (b & 1) c = 1ll * a * c % mod;
	}
	return c;
}

void add(int &a, int b) {
	a += b;
	a < mod ? NULL : a -= mod;
}

void sub(int &a, int b) {
	a -= b;
	a < 0 ? a += mod : NULL;
}

void solve(int n, int f[], int type) {
	if (!type) {
		for (int i = 1; i <= n; i++) {
			for (int j = 2 * i; j <= n; j += i) {
				sub(f[j], f[i]);
			}
		}
	} else {
		for (int i = n; i; i--) {
			for (int j = 2 * i; j <= n; j += i) {
				sub(f[i], f[j]);
			}
		}
	}
}

int main() {
	#ifdef LOCAL
		freopen("input.md", "r", stdin);
		freopen("output.md", "w", stdout);
	#endif
	scanf("%d %d %d %d", &n, &c, &d, &q);
	c %= mod - 1, d %= mod - 1;
	c = (c - d + mod - 1) % (mod - 1);
	for (int i = 1; i <= n; i++) {
		f_r[i] = qpow(i, c);
	}
	solve(n, f_r, 0);
	while (q--) {
		for (int i = 1; i <= n; i++) {
			scanf("%d", &f_z[i]);
			f_z[i] = 1ll * f_z[i] * qpow(i, mod - 1 - d) % mod;
		}
		solve(n, f_z, 0);
		bool flag = true;
		for (int i = 1; i <= n; i++) {
			if (f_r[i] > 0) {
				hx[i] = 1ll * f_z[i] * qpow(f_r[i], mod - 2) % mod;
			} else if (f_z[i] > 0) {
				flag = false;
				break;
			}
		}
		if (!flag) {
			puts("-1");
			continue;
		}
		solve(n, hx, 1);
		for (int i = 1; i <= n; i++) {
			printf("%d%c", 1ll * hx[i] * qpow(i, mod - 1 - d) % mod, " \n"[i == n]);
		}
	}
	return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值