「BZOJ2693」jzptab-莫比乌斯反演

博客围绕「BZOJ2693」jzptab问题,运用莫比乌斯反演求解。问题是求特定条件下的最小公倍数之和,通过设函数、推导关系得出求解公式,还提到可对部分式子线性筛,每次询问时数论分块求解。

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

2019.5.2

「BZOJ2693」jzptab-莫比乌斯反演

Description

∑i=1N∑j=1Mlcm(i,j)\sum_{i=1}^N\sum_{j=1}^M lcm(i,j)i=1Nj=1Mlcm(i,j)

n,m≤107n,m \leq 10^7n,m107,多组数据。

Solution

f(n)=∑i=1N∑j=1M[gcd(i,j)=1]ijf(n) = \sum_{i=1}^N\sum_{j=1}^M[gcd(i,j)=1]ijf(n)=i=1Nj=1M[gcd(i,j)=1]ij
F(n)=∑i=1N∑j=1M[n∣gcd(i,j)]ijF(n) = \sum_{i=1}^N\sum_{j=1}^M[n|gcd(i,j)]ijF(n)=i=1Nj=1M[ngcd(i,j)]ij

显然

F(n)=n2⌊Nn⌋2⌊Mn⌋2F(n) = n^2\frac{\lfloor \frac{N}{n} \rfloor}{2} \frac {\lfloor \frac{M}{n} \rfloor} {2}F(n)=n22nN2nM

又因为

F(n)=∑n∣df(d)F(n) = \sum_{n | d}f(d)F(n)=ndf(d)

所以

f(n)=∑n∣dμ(dn)F(d)f(n) = \sum_{n | d} \mu(\frac{d}{n})F(d)f(n)=ndμ(nd)F(d)

所以

∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mijgcd(i,j)=∑k=1n1kf(k)=∑k=1n1k∑k∣dμ(dk)F(d)=∑dd⌊Nd⌋2⌊Md⌋2∑k∣dμ(dk)dk=∑dd⌊Nd⌋2⌊Md⌋2∑k∣dμ(d)d\begin{matrix} & \sum_{i=1}^n\sum_{j=1}^m lcm(i,j) \\ = & \sum_{i=1}^n\sum_{j=1}^m \frac{ij}{gcd(i,j)} \\ = & \sum_{k=1}^n \frac{1}{k}f(k) \\ = & \sum_{k=1}^n \frac{1}{k} \sum_{k|d} \mu(\frac{d}{k})F(d) \\ = & \sum_d d\frac{\lfloor \frac{N}{d} \rfloor}{2} \frac {\lfloor \frac{M}{d} \rfloor} {2} \sum_{k|d} \mu(\frac{d}{k})\frac{d}{k} \\ = & \sum_d d\frac{\lfloor \frac{N}{d} \rfloor}{2} \frac {\lfloor \frac{M}{d} \rfloor} {2} \sum_{k|d} \mu(d)d \end{matrix}=====i=1nj=1mlcm(i,j)i=1nj=1mgcd(i,j)ijk=1nk1f(k)k=1nk1kdμ(kd)F(d)dd2dN2dMkdμ(kd)kddd2dN2dMkdμ(d)d

∑k∣dμ(d)d\sum_{k|d} \mu(d)dkdμ(d)d可以线性筛(可以发现每次添加一个已有的质因子时,新增加的约数的μ\muμ000,当添加一个没有的质因子ppp时,新增加的约数d⋅pd\cdot pdp的贡献为d⋅pμ(d⋅p)=dμ(d)⋅pμ(p)d \cdot p \mu(d \cdot p) = d \mu(d) \cdot p \mu(p)dpμ(dp)=dμ(d)pμ(p))。

每次询问时数论分块求就好了。

#include <bits/stdc++.h>
using namespace std;

inline int gi()
{
	char c = getchar();
	while(c < '0' || c > '9') c = getchar();
	int sum = 0;
	while('0' <= c && c <= '9') sum = sum * 10 + c - 48, c = getchar();
	return sum;
}

typedef long long ll;
const int maxn = 10000005, mod = 100000009;

int T, n, m, vis[maxn], p[maxn], tot, f[maxn], sum[maxn], s[maxn];

void pre()
{
	register int i, j;
	n = 1e7;
	f[1] = 1;
	for (i = 2; i <= n; ++i) {
		if (!vis[i]) p[++tot] = i, f[i] = mod + 1 - i;
		for (j = 1; j <= tot && i * p[j] <= n; ++j) {
			vis[i * p[j]] = 1;
			if (i % p[j]) f[i * p[j]] = (ll)f[i] * f[p[j]] % mod;
			else {f[i * p[j]] = f[i]; break;}
		}
	}
	for (i = 1; i <= n; ++i)
		sum[i] = (sum[i - 1] + (ll)i * f[i]) % mod, s[i] = s[i - 1] + i >= mod ? s[i - 1] + i - mod : s[i - 1] + i;
}

int main()
{
	freopen("jzptab.in", "r", stdin);
	freopen("jzptab.out", "w", stdout);

	int T = gi();
	register int i, j, ans;
	
	pre();
	
	while (T--) {
		n = gi(); m = gi();
		if (n > m) swap(n, m);
		ans = 0;
		for (i = 1; i <= n; i = j + 1) {
			j = min(n / (n / i), m / (m / i));
			ans = (ans + (ll)(sum[j] - sum[i - 1] + mod) * s[n / i] % mod * s[m / i]) % mod;
		}
		cout << ans << endl;
	}
	
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值