BZOJ3561: DZY Loves Math VI

Description

给定正整数n,m。求

Input

一行两个整数n,m。

Output

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

Sample Input

5 4

Sample Output

424

HINT

数据规模：

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

Source

By Jcvb

Solution

\[\begin{aligned} 原式=&\sum_{d=1}^nd^d\sum_{k=1}^{\frac{n}{d}}k^{2d}\mu(k)\sum_{i=1}^{\frac{n}{kd}}i^d\sum_{j=1}^{\frac{m}{kd}}j^d \end{aligned} \]

实际上就不用继续化简了，只需要可以均摊O(1)求出后面那个就行了。这样总复杂度就是\(O(nlogn)\)的。

考虑对于每个d，因为上界是\(\frac{n}{kd}\)，所以根据调和级数算一下枚举均摊是\(O(nlogn)\)的，所以暴力算出来就好了。

Dzy Loves Math和这题比不知道高到哪去了

#include <bits/stdc++.h>
using namespace std;

#define mod 1000000007
#define N 500010
#define ll long long
ll sum[N], c[N], ans = 0;
int cnt, n, p[N], vis[N], mu[N], m;

void init(int n) {
	mu[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!vis[i]) p[++cnt] = i, mu[i] = -1;
		for(int j = 1; j <= cnt && i * p[j] <= n; ++j) {
			vis[i * p[j]] = 1;
			if(i % p[j] == 0) break;
			mu[i * p[j]] = -mu[i];
		}
	}
}

ll power(ll a, ll b) { ll ans = 1;
	while(b) {
		if(b & 1) ans = ans * a % mod;
		a = a * a % mod; b >>= 1;
	}
	return ans;
}

int main() {
	scanf("%d%d", &n, &m);
	if(n > m) swap(n, m);
	init(m);
	for(int i = 1; i <= m; ++i) c[i] = 1ll;
	for(int d = 1; d <= n; ++d) {
		for(int j = 1; j <= m / d; ++j) {
			c[j] = 1ll * c[j] * j % mod;
			sum[j] = sum[j - 1] + c[j];
			sum[j] %= mod;
		}
		ll s = 0;
		for(int i = 1; i <= n / d; ++i) {
			s = (s + (c[i] * c[i] % mod * mu[i] * sum[n/d/i] % mod * sum[m/d/i]) % mod) % mod;
		}
		ans = (ans + s * power(d, d)) % mod;
	}
	printf("%lld\n", (ans%mod+mod)%mod);
	return 0;
}