[BZOJ 3561] DZY Loves Math VI

BZOJ传送门

Description

给定正整数$n,m$。求

Input

一行两个整数$n,m$。

Output

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

Sample Input

5 4

Sample Output

424

HINT

数据规模：

$1\le n,m\le 500000$，共有3组数据。

解题分析

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j{)}^{gcd(i,j)} \\ =\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m(\frac{ij}{d}{)}^d[gcd(i,j)=d] \\ =\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }(dij{)}^d[gcd(i,j)=1] \\ =\sum _{d=1}^n{d}^d\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)\sum _{i=1}^{\lfloor \frac{n}{dp}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{dp}\rfloor }(ij{p}^2{)}^d \\ =\sum _{d=1}^n{d}^d\sum _{p=1}^{\lfloor \frac{n}{d}\rfloor }\mu (p)p^{2d}\sum _{i=1}^{\lfloor \frac{n}{dp}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{dp}\rfloor }(ij{)}^d \end{gathered}$$

看起来不是很能再化简了？

其实我们可以发现， 枚举$d$再枚举$p$的复杂度是$O(Nln(N))$， 完全可以过这道题， 那么我们只需要再维护一个$i^d$以及它的前缀和就好了。

总复杂度$O(Nln(N))$。

代码如下：

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cstdlib>
#include <cctype>
#include <cstring>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define MX 500050
#define MOD 1000000007ll
#define ll long long
int pcnt;
int miu[MX], pri[MX];
ll mul[MX], sum[MX];
bool npr[MX];
void get()
{
	miu[1] = 1; R int i, j, tar;
	for (i = 2; i <= 5e5; ++i)
	{
		if (!npr[i]) npr[i] = true, pri[++pcnt] = i, miu[i] = -1;
		for (j = 1; j <= pcnt; ++j)
		{
			tar = i * pri[j];
			if (tar > 5e5) break;
			npr[tar] = true;
			if (!(i % pri[j])) {miu[tar] = 0; break;}
			miu[tar] = -miu[i];
		}
	}
}
IN ll fpow(ll base, R int tim)
{
	ll ret = 1;
	W (tim)
	{
		if (tim & 1) ret = ret * base % MOD;
		base = base * base % MOD, tim >>= 1;
	}
	return ret;
}
int main(void)
{
	int n, m, bdn, bdm; ll ans = 0, mt, tot;
	scanf("%d%d", &n, &m); get();
	if (n > m) std::swap(n, m);
	for (R int i = 1; i <= m; ++i) mul[i] = 1;
	for (R int i = 1; i <= n; ++i)
	{
		bdn = n / i, bdm = m / i, tot = 0;
		for (R int j = 1; j <= bdm; ++j)
		mul[j] = mul[j] * j % MOD, sum[j] = (sum[j - 1] + mul[j]) % MOD;
		for (R int j = 1; j <= bdn; ++j)
		{
			mt = miu[j] % MOD * mul[j] % MOD * mul[j] % MOD;
			(tot += mt * sum[bdn / j] % MOD * sum[bdm / j] % MOD) %= MOD;
		}
		ans = (ans + fpow(i, i) * tot % MOD) % MOD;
	}
	printf("%lld", (ans + MOD) % MOD);
}