[BZOJ3309]DZY Loves Math

题面戳我
题意：多组数据，给出n，m，求

\[\sum_{i=1}^{n}\sum_{j=1}^{m}f(\gcd(i,j)) \]

其中\(f(i)\)表示\(i\)所含质因子的最大幂指数。
例如\(f(1960)=f(2^3*5^1*7^2)=3, f(10007)=1, f(1)=0\)。
\(T\le10^5,n,m\le10^7\)

sol

首先到了这一步$$ans=\sum_{T=1}^{n}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{d|T}f(d)\mu(\frac Td)$$
所以说\(h(T)=\sum_{d|T}f(d)\mu(\frac Td)\)这个函数的前缀和怎么求？
暴力筛？\(n\le10^7!\)
所以这个时候就需要研究这个函数的性质。
我们假设\(T\)被唯一分解为\(p_1^{a_1}p_2^{a_2}...p_k^{a_k}\)
为了保证有贡献我们要让\(\mu(\frac Td)\)非零，那么\(\frac Td\)的最大指数幂不能超过一。也就是说枚举的\(d=p_1^{b_1}p_2^{b_2}...p_k^{b_k}\)必须满足，对于所有\(i\in[1,k]\)有\(0\le a_i-b_i\le 1\)。那么这么来说对于每个\(T\)，设其一共含有\(k\)种不同的质因数，那么会给它贡献的就只有\(2^k\)个数。枚举的每个\(d\)，就有\(f(d)=max\{b_i\}\)
分两种情况讨论：
1、存在一组\(i,j\)使\(a_i< a_j\)。那么此时\(b_i\)就永远不可能成为\(\max\{b\}\)，换言之，无论你取不取这个\(i\)，\(f(d)\)的值算出来都是一样的。又因为取这个\(i\)的时候的\(\mu(\frac Td)\)值与不取这个\(i\)的时候的\(\mu(\frac Td)\)值恰好相反，所以加起来就会有\(h(T)=0\)。
2、所有\(a_i\)都相等。还是装作上一种情况算，不同之处在于当所有\(b_i\)取\(a_i-1\)的时候\(f(d)=b_i=a_i-1\)，那么这时\(h(T)\)的值就只与这个-1前面的\(\mu\)值有关（因为其他的都抵消掉了）。易知这个-1前面的系数是\(\mu(p_1p_2...p_k)\)，所以说\(h(T)=-(-1)^k=(-1)^{k+1}\)
线性筛：记录每个数的最小质因子的幂指数和最小质因子的指数次幂，判除掉最小质因子后最小值因子的幂指数是否相等（好像写的很复杂的样子。。。）

code

#include<cstdio>
#include<algorithm>
using namespace std;
const int N = 10000000;
int gi()
{
	int x=0,w=1;char ch=getchar();
	while ((ch<'0'||ch>'9')&&ch!='-') ch=getchar();
	if (ch=='-') w=0,ch=getchar();
	while (ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+ch-'0',ch=getchar();
	return w?x:-x;
}
int pri[N>>2],tot,zhi[N+5],a[N+5],low[N+5],h[N+5];
void Sieve()
{
	zhi[1]=1;
	for (int i=2;i<=N;i++)
	{
		if (!zhi[i]) low[i]=pri[++tot]=i,a[i]=h[i]=1;
		for (int j=1;j<=tot&&i*pri[j]<=N;j++)
		{
			zhi[i*pri[j]]=1;
			if (i%pri[j]==0)
			{
				a[i*pri[j]]=a[i]+1;
				low[i*pri[j]]=low[i]*pri[j];
				if (i==low[i])
					h[i*pri[j]]=1;
				else h[i*pri[j]]=(a[i/low[i]]==a[i*pri[j]])?-h[i/low[i]]:0;
				break;
			}
			a[i*pri[j]]=1;
			low[i*pri[j]]=pri[j];
			h[i*pri[j]]=(a[i]==1)?-h[i]:0;
		}
	}
	for (int i=1;i<=N;i++) h[i]+=h[i-1];
}
int main()
{
	Sieve();
	int T=gi();
	while (T--)
	{
		int n=gi(),m=gi();
		if (n>m) swap(n,m);
		int i=1;long long ans=0;
		while (i<=n)
		{
			int j=min(n/(n/i),m/(m/i));
			ans+=1ll*(n/i)*(m/i)*(h[j]-h[i-1]);
			i=j+1;
		}
		printf("%lld\n",ans);
	}
	return 0;
}