bzoj3561 DZY Loves Math VI 莫比乌斯反演

Description

---

求$$\sum _{i=1}^n\sum _{j=1}^m{lcm(i,j)}^{gcd(i,j)}$$

Solution

---

化柿子最终可以得到
$$ans=\sum _{d=1}^n{d}^d\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\mu (x)x^{2d}\sum _{i=1}^{\lfloor \frac{n}{dx}\rfloor }{i}^d\sum _{j=1}^{\lfloor \frac{m}{dx}\rfloor }{j}^d$$

我们枚举d，然后枚举d的倍数x，直接暴力做是nlogn的

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
const int MOD=1000000007;
const int N=500005;

int prime[N]; LL sum[N],cm[N];

bool not_prime[N];

short mu[N];

void pre_work(int n) {
	mu[1]=1;
	rep(i,2,n) {
		if (!not_prime[i]) {
			prime[++prime[0]]=i;
			mu[i]=-1;
		}
		for (int j=1;j<=prime[0]&&i*prime[j]<=n;j++) {
			not_prime[i*prime[j]]=1;
			if (i%prime[j]==0) {
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
}

LL ksm(LL x,LL dep) {
	LL res=1;
	for (;dep;dep>>=1) {
		(dep&1)?(res=res*x%MOD):0;
		x=x*x%MOD;
	}
	return res;
}

int main(void) {
	pre_work(N-1);
	int n,m; scanf("%d%d",&n,&m);
	if (n>m) std:: swap(n,m);
	rep(i,1,m) cm[i]=1;
	LL ans=0;
	rep(d,1,n) {
		for (int x=1;d*x<=m;++x) {
			cm[x]=cm[x]*x%MOD;
			sum[x]=(sum[x-1]+cm[x])%MOD;
		}
		LL wjp=ksm(d,d);
		for (int x=1;d*x<=n;++x) {
			ans=(ans+1LL*mu[x]*cm[x]%MOD*cm[x]%MOD*wjp%MOD*sum[n/d/x]%MOD*sum[m/d/x]%MOD+MOD)%MOD;
		}
	}
	printf("%lld\n", ans);
	return 0;
}