luogu3911 最小公倍数之和

题目大意

给出一些数(A_1,A_2,cdots A_n)，求

[sum_{i=1}^{n}sum_{j=1}^{n}mathrm{lcm}(A_i,A_j) ]

(A_i,A_nleq 50000)

运用莫比乌斯反演思路

对于这种对多个数进行gcd、lcm统计的题，往往要用莫比乌斯反演。运用莫比乌斯反演的思路往往如下：

1. 我们要求的(g(x))是什么？
2. 比较容易求的(f(x))是什么？
3. 如果我们要求的(g(x))已知，则比较容易求的(f(x))应当如何表达？
4. 如果表达是以莫比乌斯反演公式的形式，则先求出(f(x))，然后反演出(g(x))即可。

我们要求的(g(x))是什么？

错误做法

根据我们的做题经验，(g(x))表示最大公约数是(x)的数的对数，即

[g(x)=sum_{i=1}^nsum_{j=1}^n [gcd(A_i,A_j)=k] ]

为什么可以利用它呢？因为

[原式=sum_{i=1}^{n}sum_{j=1}^{n}frac{A_i A_j}{gcd(A_i,A_j)} ag{1} ]

提取出(A_i,A_j)得

[原式=(sum_{i=1}^n A_i)^2sum_{i=1}^{n}sum_{j=1}^{n}frac{1}{gcd(A_i,A_j)} ]

OH,NO!这么化简是不对的。设(f(x),g(x))为任意函数，则

[sum_{i=1}^{n}sum_{j=1}^{n}f(i)g(j)=sum_{i=1}^n f(i)sum_{j=1}^n g(j) ]

此式成立，因为函数(f(i),g(j))的参数只关于一个变量。但是

[sum_{i=1}^nsum_{j=1}^n f(i,j)g(i,j) eq sum_{i=1}^nsum_{j=1}^n f(i,j)sum_{i=1}^nsum_{j=1}^n g(i,j) ]

这就很荒谬了。函数(f,g)是同时关于(i,j)的函数。两个函数相乘时，里面的((i,j))都应当是相等的，但化后的式子(f,g)内的(i,j)不相等时也乘起来了，这就错了。原式中，(f(i,j)=A_i A_j)，(g(i,j)=frac{1}{gcd(A_i,A_j)})。问题就出在这里。

正确做法

至少（1）式还是对的。因为(gcd(A_i,A_j))一定时，我们要求的是(A_i A_j)的和，所以

[g(x)=sum_{i=1}^nsum_{j=1}^n A_i A_j[gcd(A_i,A_j)=x] ]

(f(x))怎么求？

定义

[f(x)=sum_{i=1}^nsum_{j=1}^n A_i A_j[gcd(A_i,A_j)=kx] ag{2}$$$$=sum_{i=1}^nsum_{j=1}^n A_i A_j[x|A_i,x|A_j]$$$$=sum_{x|A_i}sum_{x|A_j}A_i A_j$$$$=(sum_{x|A_i}A_i)^2 ag{3} ]

(2)式即能体现莫比乌斯函数的性质。
运用（3）求(f(x))。

如何迅速地找到所有满足条件的(A_i)？

建立一个维护A个数的桶数组即可。

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

#define ll long long

const int MAX_N = 50010;
int Mu[MAX_N];
ll F[MAX_N], G[MAX_N], ExistCnt[MAX_N];
ll N, MaxA;

void GetMu(int *mu, int n)
{
	static int prime[MAX_N];
	static bool NotPrime[MAX_N];
	int primeCnt = 0;
	Mu[1] = 1;
	for (int i = 2; i <= n; i++)
	{
		if (!NotPrime[i])
		{
			prime[primeCnt++] = i;
			Mu[i] = -1;
		}
		for (int j = 0; j <= primeCnt; j++)
		{
			if (i*prime[j] > n)
				break;
			NotPrime[i*prime[j]] = true;
			if (i%prime[j] == 0)
			{
				mu[i*prime[j]] = 0;
				break;
			}
			else
				mu[i*prime[j]] = -mu[i];
		}
	}
}

void GetF()
{
	for (int cd = 1; cd <= MaxA; cd++)
	{
		ll sum = 0;
		for (int k = 1; k <= MaxA / cd; k++)
			sum += cd*k*ExistCnt[cd*k];
		F[cd] = sum*sum;
	}
}

void GetG()
{
	for (int gcd = 1; gcd <= MaxA; gcd++)
		for (int k = 1; k <= MaxA / gcd; k++)
			G[gcd] += F[gcd*k] * Mu[k];
}

ll Solve()
{
	ll ans = 0;
	for (int gcd = 1; gcd <= MaxA; gcd++)
		ans += G[gcd] / gcd;
	return ans;
}

int main()
{
	ll a;
	scanf("%lld", &N);
	for (int i = 1; i <= N; i++)
	{
		scanf("%lld", &a);
		ExistCnt[a]++;
		MaxA = max(MaxA, a);
	}
	GetMu(Mu, MaxA);
	GetF();
	GetG();
	cout << Solve() << endl;
	return 0;
}