BZOJ3561 DZY Loves Math VI 【莫比乌斯反演】

题目

给定正整数n,m。求

输入格式

一行两个整数n,m。

输出格式

一个整数，为答案模1000000007后的值。

输入样例

5 4

输出样例

424

提示

数据规模：

1<=n,m<=500000，共有3组数据。

题解

套路反演：

\[\begin{aligned} ans &= \sum\limits_{d = 1}^{n} \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j = 1}^{\lfloor \frac{m}{d} \rfloor} (\frac{id * jd}{d})^{d} \qquad [gcd(i,j) == 1] \\ &= \sum\limits_{d = 1}^{n} d^{d} \sum\limits_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum\limits_{j = 1}^{\lfloor \frac{m}{d} \rfloor} i^{d}j^{d} \qquad [gcd(i,j) == 1] \\ &= \sum\limits_{d= 1}^{n} d^{d} \sum\limits_{k = 1}^{\lfloor \frac{n}{d} \rfloor} \mu(k) k^{2d} \sum\limits_{i = 1}^{\lfloor \frac{n}{kd} \rfloor} i^{d} \sum\limits_{j = 1}^{\lfloor \frac{m}{kd} \rfloor} j^{d} \end{aligned} \]

然后就慌了，好像搞不下去了
仔细分析一下复杂度，对于每个\(d\)，里面那堆玩意只需要\(O(\lfloor \frac{m}{d} \rfloor)\)就可以计算出来
所以复杂度为\(n\)乘一个调和级数
即\(O(nlogn)\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 500005,maxm = 100005,INF = 1000000000,P = 1e9 + 7;
int p[maxn],isn[maxn],pi,mu[maxn],n,m;
LL a[maxn],sum[maxn];
LL qpow(LL a,LL b){
	LL ans = 1;
	for (; b; b >>= 1, a = a * a % P)
		if (b & 1) ans = ans * a % P;
	return ans;
}
void init(){
	mu[1] = 1;
	for (int i = 2; i <= n; i++){
		if (!isn[i]) p[++pi] = i,mu[i] = -1;
		for (int j = 1; j <= pi && i * p[j] <= n; j++){
			isn[i * p[j]] = true;
			if (i % p[j] == 0){
				mu[i * p[j]] = 0;
				break;
			}
			mu[i * p[j]] = -mu[i];
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);
	if (n > m) swap(n,m);
	init();
	REP(i,m) a[i] = 1;
	LL ans = 0;
	for (int d = 1; d <= n; d++){
		for (LL i = 1; i <= m / d; i++){
			a[i] = a[i] * i % P;
			sum[i] = (sum[i - 1] + a[i]) % P;
		}
		LL tmp = 0;
		for (int k = 1; k <= n / d; k++){
			tmp = (tmp + (mu[k] * qpow(k,2 * d) % P + P) % P * sum[n / (k * d)] % P * sum[m / (k * d)] % P);
			tmp = (tmp + P) % P;
		}
		ans = (ans + qpow(d,d) * tmp % P) % P;
	}
	ans = (ans % P + P) % P;
	printf("%lld\n",ans);
	return 0;
}