bzoj 3561: DZY Loves Math VI 莫比乌斯反演

题目大意：
给出$n$，$m$，求
$$\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j{)}^{gcd(i,j)}$$

分析：
直接推式子，
$$\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}(ijd{)}^d[gcd(i,j)==1]$$
$$\sum _{d=1}^n{d}^d\sum _{k=1}^{n/d}\mu (k)\ast \sum _{i=1}^{n/d/k}(ik{)}^d\ast \sum _{j=1}^{m/d/k}(jk{)}^d$$
$$\sum _{d=1}^n{d}^d\sum _{k=1}^{n/d}\mu (k)\ast k^{2d}\ast \sum _{i=1}^{n/d/k}{i}^d\ast \sum _{j=1}^{m/d/k}{j}^d$$

直接暴力算，自然数幂和以及幂次可以预处理，根据调和级数这个是$O(nlogn)$的。

代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const LL mod=1000000007;
const int maxn=5e5+7;

using namespace std;

int n,m,cnt;
LL ans,f[maxn],sum[maxn];
int not_prime[maxn],prime[maxn],mu[maxn];

LL power(LL x,LL y)
{
	if (y==1) return x;
	LL c=power(x,y/2);
	c=(c*c)%mod;
	if (y&1) c=(c*x)%mod;
	return c;
}

void getmu(int n)
{
	mu[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!not_prime[i])
		{
			prime[++cnt]=i;
			mu[i]=mod-1;
		}
		for (int j=1;j<=cnt;j++)
		{
			if (i*prime[j]>n) break;
			not_prime[i*prime[j]]=1;
			if (i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			else mu[i*prime[j]]=mod-mu[i];
		}
	}
}

int main()
{
	getmu(5e5);
	scanf("%d%d",&n,&m);
	if (n>m) swap(n,m);
	for (int i=1;i<=m;i++) f[i]=1;
	for (int i=1;i<=n;i++)
	{
		LL s=power(i,i);
		for (int j=1;j<=m/i;j++)
		{
			f[j]=f[j]*j%mod;
			sum[j]=(sum[j-1]+f[j])%mod;
		}
		for (int j=1;j<=n/i;j++)
		{
			ans=(ans+s*(LL)mu[j]%mod*f[j]%mod*f[j]%mod*sum[n/i/j]%mod*sum[m/i/j]%mod)%mod;
		}
	}
	printf("%lld\n",ans);
}