BZOJ2671 Calc 【莫比乌斯反演】

题目链接

BZOJ2671

题解

令(d = (a,b))，(a = dx,b = dy)
那么有

[egin{aligned} d(x + y) | d^2xy \ (x + y) | dxy end{aligned} ]

由于(x perp y)，所以

[(x + y) | d ]

令(d = k(x + y))，那么较大的(b = yk(x + 1))要满足(b <= n)
那么对于一组互质的((x,y))，合法的(d)的个数有

[lfloor frac{n}{y(x + y)} floor ]

那么就有式子

[egin{aligned} ans &= sumlimits_{x = 1}^{sqrt{n}} sumlimits_{y = x + 1}^{sqrt{n}}[(x,y) = 1] lfloor frac{n}{y(x + y)} floor \ &= sumlimits_{x = 1}^{sqrt{n}} sumlimits_{y = x + 1}^{sqrt{n}}lfloor frac{n}{y(x + y)} floor sumlimits_{d | (x,y)} mu(d) \ &= sumlimits_{d = 1}^{sqrt{n}} mu(d) sumlimits_{x = 1}^{lfloor frac{n}{d} floor} sumlimits_{y = x + 1}^{lfloor frac{n}{d} floor} lfloor frac{n}{y(x + y)} floor\ &= sumlimits_{d = 1}^{sqrt{n}} mu(d) sumlimits_{y = 1}^{lfloor frac{n}{d} floor}sumlimits_{x = 1}^{y - 1} lfloor frac{n}{y(x + y)} floor\ &= sumlimits_{d = 1}^{sqrt{n}} mu(d) sumlimits_{y = 1}^{lfloor frac{n}{d} floor}sumlimits_{x = y + 1}^{min(2y - 1,lfloor frac{n}{d} floor)} lfloor frac{n}{d^2 yx} floor end{aligned} ]

里边可以整除分块，复杂度(O(?n^{frac{3}{4}}))

#include<algorithm>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<cstdio>
#include<vector>
#include<queue>
#include<cmath>
#include<map>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define cls(s,v) memset(s,v,sizeof(s))
#define mp(a,b) make_pair<int,int>(a,b)
#define cp pair<int,int>
using namespace std;
const int maxn = 100005,maxm = 100005,INF = 0x3f3f3f3f;
inline int read(){
	int out = 0,flag = 1; char c = getchar();
	while (c < 48 || c > 57){if (c == '-') flag = 0; c = getchar();}
	while (c >= 48 && c <= 57){out = (out << 1) + (out << 3) + c - 48; c = getchar();}
	return flag ? out : -out;
}
LL ans;
int n,m;
int p[maxn],isn[maxn],mu[maxn],pi,N = 100000;
void init(){
	mu[1] = 1;
	for (int i = 2; i <= N; i++){
		if (!isn[i]) p[++pi] = i,mu[i] = -1;
		for (int j = 1; j <= pi && i * p[j] <= N; j++){
			isn[i * p[j]] = true;
			if (i % p[j] == 0) break;
			mu[i * p[j]] = -mu[i];
		}
	}
}
LL cal(int n,int m){
	LL re = 0;
	for (int y = 1; y <= m; y++){
		int t = n / y;
		for (int x = y + 1,nxt; x < (y << 1) && x <= t; x = nxt + 1){
			nxt = min(t / (t / x),(y << 1) - 1);
			re += 1ll * (nxt - x + 1) * (t / x);
		}
	}
	return re;
}
int main(){
	n = read(); m = (int)sqrt(n);
	init();
	for (int d = 1; d <= m; d++) if (mu[d]) ans += mu[d] * cal(n / d / d,m / d);
	printf("%lld
",ans);
	return 0;
}