【BZOJ】3561: DZY Loves Math VI

题意

求(sum_{i=1}^{n} sum_{j=1}^{m} lcm(i, j)^{gcd(i, j)})（(n, m<=500000)）

分析

很显然要死推莫比乌斯

题解

设(n le m)

[egin{aligned} ans & = sum_{i=1}^{n} sum_{j=1}^{m} lcm(i, j)^{gcd(i, j)} \ & = sum_{i=1}^{n} sum_{j=1}^{m} (frac{ij}{gcd(i, j)})^{gcd(i, j)} \ & = sum_{d=1}^{n} sum_{i=1}^{a} sum_{j=1}^{b} left( frac{ijdd}{d} ight)^{d} sum_{k|(i, j)} mu(k) left( a=left lfloor frac{n}{d} ight floor, b=left lfloor frac{m}{d} ight floor ight) \ & = sum_{d=1}^{n} d^d sum_{k=1}^{a} mu(k) sum_{k|i}^{a} i^d sum_{k|j}^{b} j^d \ & = sum_{d=1}^{n} d^d sum_{k=1}^{a} mu(k) k^{2d} sum_{i=1}^{left lfloor frac{a}{k} ight floor} i^d sum_{j=1}^{left lfloor frac{b}{k} ight floor} j^d \ & = sum_{d=1}^{n} d^d sum_{k=1}^{left lfloor frac{n}{d} ight floor} mu(k) k^{2d} sum_{i=1}^{left lfloor frac{n}{kd} ight floor} i^d sum_{j=1}^{left lfloor frac{m}{kd} ight floor} j^d \ end{aligned} ]

于是我们对于每一个(d)，暴力维护一下(mu(k) k^{2d})，暴力维护一下(displaystyle sum_{i=1}^{left lfloor frac{m}{kd} ight floor} j^d)，总复杂度(O(nlogn))

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mo=1000000007, N=500005;
int mu[N], p[N], pcnt, np[N], c[N], C[N], b[N];
int ipow(int a, int b) {
	int x=1;
	for(; b; b>>=1, a=(ll)a*a%mo) if(b&1) x=(ll)x*a%mo;
	return x;
}
void init(int n) {
	mu[1]=1;
	for(int i=2; i<=n; ++i) {
		if(!np[i]) {
			p[pcnt++]=i;
			mu[i]=-1;
		}
		for(int j=0; j<pcnt; ++j) {
			int t=p[j]*i;
			if(t>n) break;
			np[t]=1;
			if(i%p[j]==0) {
				mu[t]=0;
				break;
			}
			mu[t]=-mu[i];
		}
	}
}
int main() {
	int n, m, ans=0;
	scanf("%d%d", &n, &m);
	if(n>m) {
		swap(n, m);
	}
	init(n);
	for(int i=1; i<=m; ++i) {
		c[i]=1;
	}
	for(int d=1; d<=n; ++d) {
		int A=ipow(d, d);
		int nn=n/d, mm=m/d;
		for(int k=1; k<=mm; ++k) {
			c[k]=(ll)c[k]*k%mo;
			C[k]=C[k-1]+c[k];
			if(C[k]>=mo) {
				C[k]-=mo;
			}
		}
		int temp=0;
		for(int k=1; k<=nn; ++k) if(mu[k]) {
			temp+=(ll)c[k]*c[k]%mo*C[nn/k]%mo*C[mm/k]%mo*mu[k];
			if(temp>=mo) {
				temp-=mo;
			}
			if(temp<0) {
				temp+=mo;
			}
		}
		ans+=(ll)A*temp%mo;
		if(ans>=mo) {
			ans-=mo;
		}
	}
	printf("%d
", ans);
	return 0;
}