SDOI2018 舊試題(莫比烏斯反演+三元環計數)

WAautomaton發表於2019-02-14

題目連結

題目大意

求:i=1Aj=1Bk=1Cd(ijk)\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^Cd(ijk)
其中d(i)d(i)ii的因數個數。A,B,C105A,B,C\le 10^5

題解

居然真的會有這麼美妙的結論qwq……
d(ijk)=aibjck[gcd(a,b)=1][gcd(a,c)=1][gcd(b,c)=1]d(ijk)=\sum_{a|i}\sum_{b|j}\sum_{c|k}[\gcd(a,b)=1][\gcd(a,c)=1][\gcd(b,c)=1]
感性理解一下就是,每個質因數分開考慮。假設它在i,j,ki,j,k中的指數分別為x,y,zx,y,z,那麼它總共會在上面的算式中出現恰好x+y+z+1x+y+z+1次。由於所有質數貢獻獨立,因此上面的算式就等於因數個數。
暴力帶入原式並使用莫比烏斯反演,可以得到:
i=1Aj=1Bk=1Caibjckxa,xbμ(x)ya,ycμ(y)zb,zcμ(z)=i=1Aj=1Bk=1Cxi,xjyi,ykzj,zkμ(x)μ(y)μ(z)d(ilcm(x,y))d(jlcm(x,z))d(klcm(y,z))=x=1Ay=1Az=1Bμ(x)μ(y)μ(z)D(Alcm(x,y))D(Blcm(x,z))D(Clcm(y,z))\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{a|i}\sum_{b|j}\sum_{c|k}\sum_{x|a,x|b}\mu(x)\sum_{y|a,y|c}\mu(y)\sum_{z|b,z|c}\mu(z) \\ =\sum_{i=1}^A\sum_{j=1}^B\sum_{k=1}^C\sum_{x|i,x|j}\sum_{y|i,y|k}\sum_{z|j,z|k}\mu(x)\mu(y)\mu(z)d\left(\frac{i}{lcm(x,y)}\right)d\left(\frac{j}{lcm(x,z)}\right)d\left(\frac{k}{lcm(y,z)}\right) \\ =\sum_{x=1}^A\sum_{y=1}^A\sum_{z=1}^B\mu(x)\mu(y)\mu(z)D\left(\frac{A}{lcm(x,y)}\right)D\left(\frac{B}{lcm(x,z)}\right)D\left(\frac{C}{lcm(y,z)}\right)
其中D(i)=j=1id(i)D(i)=\sum_{j=1}^id(i)
考慮x,y,zx,y,z全都不相等的情況。我們可以建無向圖,對於任意的x<y,lcm(x,y)Cx<y,lcm(x,y)\le C,連邊(x,y)(x,y)。然後找圖中的三元環,並任意輪換3!次加入答案。
如果x,y,zx,y,z中恰好有一對相等,則可以任意輪換3次計入答案。
如果全部相等,直接迴圈列舉即可。
三元環計數在本題中似乎暴力就能過(大霧)。

#include <bits/stdc++.h>
namespace IOStream {
	const int MAXR = 1 << 23;
	char _READ_[MAXR], _PRINT_[MAXR];
	int _READ_POS_, _PRINT_POS_, _READ_LEN_;
	inline char readc() {
	#ifndef ONLINE_JUDGE
		return getchar();
	#endif
		if (!_READ_POS_) _READ_LEN_ = fread(_READ_, 1, MAXR, stdin);
		char c = _READ_[_READ_POS_++];
		if (_READ_POS_ == MAXR) _READ_POS_ = 0;
		if (_READ_POS_ > _READ_LEN_) return 0;
		return c;
	}
	template<typename T> inline void read(T &x) {
		x = 0; register int flag = 1, c;
		while (((c = readc()) < '0' || c > '9') && c != '-');
		if (c == '-') flag = -1; else x = c - '0';
		while ((c = readc()) >= '0' && c <= '9') x = x * 10 + c - '0';
		x *= flag;
	}
	template<typename T1, typename ...T2> inline void read(T1 &a, T2 &...x) {
		read(a), read(x...);
	}
	inline int reads(char *s) {
		register int len = 0, c;
		while (isspace(c = readc()) || !c);
		s[len++] = c;
		while (!isspace(c = readc()) && c) s[len++] = c;
		s[len] = 0;
		return len;
	}
	inline void ioflush() {
		fwrite(_PRINT_, 1, _PRINT_POS_, stdout), _PRINT_POS_ = 0;
		fflush(stdout);
	}
	inline void printc(char c) {
		_PRINT_[_PRINT_POS_++] = c;
		if (_PRINT_POS_ == MAXR) ioflush();
	}
	inline void prints(char *s) {
		for (int i = 0; s[i]; i++) printc(s[i]);
	}
	template<typename T> inline void print(T x, char c = '\n') {
		if (x < 0) printc('-'), x = -x;
		if (x) {
			static char sta[20];
			register int tp = 0;
			for (; x; x /= 10) sta[tp++] = x % 10 + '0';
			while (tp > 0) printc(sta[--tp]);
		} else printc('0');
		printc(c);
	}
	template<typename T1, typename ...T2> inline void print(T1 x, T2... y) {
		print(x, ' '), print(y...);
	}
}
using namespace IOStream;
using namespace std;
typedef long long ll;
typedef pair<int, int> P;
#define cls(a) memset(a, 0, sizeof(a))

const int MAXN = 100005, MOD = 1000000007;
int mu[MAXN], prime[MAXN], vis[MAXN], A, B, C = 100000, T;
ll d[MAXN], da[MAXN], db[MAXN], dc[MAXN];
vector<P> g[MAXN];
int main() {
	mu[1] = 1;
	for (int i = 2, cnt = 0; i <= C; i++) {
		if (!vis[i]) mu[prime[++cnt] = i] = -1;
		for (int j = 1; j <= cnt; j++) {
			int p = prime[j], q = i * p;
			if (q > C) break;
			vis[q] = 1;
			if (i % p == 0) break;
			mu[q] = -mu[i];
		}
	}
	for (int i = 1; i <= C; i++) {
		for (int j = i; j <= C; j += i) ++d[j];
		d[i] += d[i - 1];
	}
	for (scanf("%d", &T); T--;) {
		scanf("%d%d%d", &A, &B, &C);
		cls(vis);
		if (A > B) swap(A, B);
		if (A > C) swap(A, C);
		if (B > C) swap(B, C);
		for (int i = 1; i <= C; i++) {
			da[i] = d[A / i];
			db[i] = d[B / i];
			dc[i] = d[C / i];
			g[i].clear();
		}
		ll ans = 0;
		for (int i = 1; i <= C; i++) {
			for (int j = 1; j * i <= C; j++) {
				for (int k = j + 1; (ll)i * j * k <= C; k++) {
					int a = i * j, b = i * k, c = i * j * k;
					if (!mu[a] || !mu[b] || __gcd(j, k) > 1) continue;
					ans += mu[a] * mu[a] * mu[b] * (
						da[c] * db[c] * dc[a] +
						da[c] * db[a] * dc[c] +
						da[a] * db[c] * dc[c] );
					ans += mu[a] * mu[b] * mu[b] * (
						da[b] * db[c] * dc[c] +
						da[c] * db[b] * dc[c] +
						da[c] * db[c] * dc[b] );
					g[a].push_back(P(b, c));
				}
			}
			ans += mu[i] * mu[i] * mu[i] * da[i] * db[i] * dc[i];
		}
		for (int i = 1; i <= C; i++) if (mu[i]) {
			for (P j : g[i]) vis[j.first] = i, prime[j.first] = j.second;
			for (P j : g[i]) if (j.first > i && mu[j.first]) {
				int a = j.second;
				for (P k : g[j.first]) if (k.first > j.first && vis[k.first] == i && mu[k.first]) {
					int b = k.second, c = prime[k.first];
					ans += mu[i] * mu[j.first] * mu[k.first] * (
						da[a] * db[b] * dc[c] +
						da[a] * db[c] * dc[b] +
						da[b] * db[a] * dc[c] +
						da[b] * db[c] * dc[a] +
						da[c] * db[a] * dc[b] +
						da[c] * db[b] * dc[a] );
				}
			}
		}
		printf("%lld\n", (ans % MOD + MOD) % MOD);
	}
	return 0;
}

相關文章