[P5348]密码解锁

Description

给一个长度为 (n) 的数组 (a[1dots n]) ，满足 (sum_{m|x}a[x] = mu(m))，求 (a[m])。

(nle 10^{18}, mle 10^9, frac{n}{m}le10^9,ngeq m)

Solution

由另一种形式的莫比乌斯反演：

[egin{aligned} a[m] &= sum_{m|x}mu(frac{x}{m})mu(x)\ &=sum_{i=1}^{frac{n}{m}}mu(i)mu(im)\ &=mu(m)sum_{i=1}^{lfloorfrac{n}{m} floor}mu(i)^2[gcd(i, m) = 1]\ end{aligned} ]

后面那个 (sum) 就是在求 (1dots frac{n}{m}) 中与 (m) 互质且不能写成完全平方数的倍数的个数。

类似于 [中山市选2011]完全平方数，可以容斥求：

令 (N = lfloorfrac{n}{m} floor)，

[egin{aligned} a[m] &=mu(m)sum_{i=1}^{frac{n}{m}}mu(i)^2[gcd(i, m) = 1]\ &=mu(m)sum_{i=1}^{sqrt{N}}mu(i)sum_{j=1}^{lfloorfrac{N}{i^2} floor}[gcd(i^2j,m)=1]\ &=mu(m)sum_{i=1}^{sqrt{N}}mu(i)[gcd(i,m)=1]sum_{j=1}^{lfloorfrac{N}{i^2} floor}[gcd(j,m)=1]\ &=mu(m)sum_{i=1}^{sqrt{N}}mu(i)[gcd(i,m)=1]sum_{j=1}^{lfloorfrac{N}{i^2} floor}sum_{d|gcd(j,m)}mu(d)\ &=mu(m)sum_{i=1}^{sqrt{N}}mu(i)[gcd(i,m)=1]sum_{d|m}mu(d)lfloorfrac{lfloorfrac{N}{i^2} floor}{d} floor end{aligned} ]

然后就可以把 (m) 的所有约数处理出来，暴力算（复杂度上界为 (O(Tsqrt frac{n}{m} sqrt m))，实际后面的 (sqrt m) 跑不满）。

注意(mu)要筛到(sqrt N)复杂度才是对的（不然多一个根号）。

code

#include <bits/stdc++.h>

typedef long long LL;
typedef unsigned long long uLL;

#define SZ(x) ((int)x.size())
#define ALL(x) (x).begin(), (x).end()
#define MP(x, y) std::make_pair(x, y)
#define DEBUG(...) fprintf(stderr, __VA_ARGS__)
#define GO cerr << "GO" << endl;

using namespace std;

inline void proc_status()
{
	ifstream t("/proc/self/status");
	cerr << string(istreambuf_iterator<char>(t), istreambuf_iterator<char>()) << endl;
}

template<class T> inline T read() 
{
	register T x(0);
	register char c;
	register int f(1);
	while (!isdigit(c = getchar())) if (c == '-') f = -1;
	while (x = (x << 1) + (x << 3) + (c xor 48), isdigit(c = getchar()));
	return x * f;
}

template<typename T> inline bool chkmin(T &a, T b) { return a > b ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }

const int maxN = 1e5;

bool vis[maxN + 1];
vector<int> prime;
int mu[maxN + 1];
LL n, m;

void Init()
{
	mu[1] = 1;
	for (int i = 2; i <= maxN; ++i)
	{
		if (!vis[i])
		{
			prime.push_back(i);
			mu[i] = -1;
		}
		for (int j = 0; j < SZ(prime) and prime[j] * i <= maxN; ++j)
		{
			vis[prime[j] * i] = 1;
			if (i % prime[j] == 0)
				break;
			else 
				mu[i * prime[j]] = -mu[i];
		}
	}
}

int Mu(LL M)
{
	if (M <= maxN) return mu[M];
	//GO;
	int cnt = 0;
	for (LL i = 2; i * i <= M; ++i)
		if (M % i == 0)
		{
			M /= i;
			cnt++;
			if (M % i == 0)
				return 0;
		}
	if (M != 1) cnt++;
	return (cnt & 1) ? -1 : 1;
}

vector<pair<LL, int> > p;

LL calc(LL i)
{
	LL ans(0);
	for (int l = 0; l < SZ(p); ++l)
		ans += (LL)p[l].second * (((n / m) / (i * i)) / p[l].first);
	return ans;
}

void GetP(LL m)
{
	p.clear();
	for (LL d = 1; d * d <= m; ++d)
		if (m % d == 0)
		{
			p.push_back(MP(d, Mu(d)));
			if (m / d == d) continue;
			p.push_back(MP(m / d, Mu(m / d)));
		}
}

void Solve()
{
	int T = read<int>();
	while (T--)
	{
		n = read<LL>(), m = read<LL>();

		GetP(m);

		LL ans(0);
		for (LL i = 1; i * i <= (n / m); ++i)
			if (__gcd((LL)i, m) == 1)
				ans += Mu(i) * calc(i);
		ans *= Mu(m);

		printf("%lld
", ans);
	}
}

int main() 
{
#ifndef ONLINE_JUDGE
	freopen("xhc.in", "r", stdin);
	freopen("xhc.out", "w", stdout);
#endif
	Init();
	Solve();
	return 0;
}