Description：

$\begin{align}\sum_{i=1}^N\sum_{j=1}^Mlcm(i,j)^{gcd(i,j)}\end{align}$

$1\le N,M\le 500000$

Solution：

$\begin{align}\sum_{i=1}^N\sum_{j=1}^Mlcm(i,j)^{gcd(i,j)}\end{align}$

$\begin{align}=\sum_d^{min(N,M)}\sum_{i=1}^N\sum_{j=1}^M(\frac {ij} d)^d[gcd(i,j)==d]\end{align}$

$\begin{align}=\sum_d^{min(N,M)}d^d\sum_{i=1}^{\lfloor \frac N d\rfloor}i^d\sum_{j=1}^{\lfloor \frac M d\rfloor}j^d\sum_{k|i,k|j}\mu(k)\end{align}$

$\begin{align}=\sum_d^{min(N,M)}d^d\sum_k^{min(\lfloor\frac N d\rfloor,\lfloor\frac M d\rfloor)}\mu(k)\sum_{i=1,k|i}^{\lfloor \frac N d\rfloor}i^d\sum_{j=1,k|j}^{\lfloor \frac M d\rfloor}j^d\end{align}$

$\begin{align}=\sum_d^{min(N,M)}d^d\sum_k^{min(\lfloor\frac N d\rfloor,\lfloor\frac M d\rfloor)}\mu(k)\sum_{i=1}^{\lfloor \frac N {dk}\rfloor}(ik)^d\sum_{j=1}^{\lfloor \frac M {dk}\rfloor}(jk)^d\end{align}$

$\begin{align}=\sum_d^{min(N,M)}d^d\sum_{k|d}^{min(\lfloor\frac N d\rfloor,\lfloor\frac M d\rfloor)}\mu(k)k^{2d}\times\sum_{i=1}^{\lfloor \frac N {dk}\rfloor}i^d\times\sum_{j=1}^{\lfloor \frac M {dk}\rfloor}j^d\end{align}$

后面那个式子可以递推

Code：

#include

#include

#include

#include

#include

#include

using namespace std;

int n,m;

#define MAXN 500010

#define MOD 1000000007

int prime[MAXN],tot = 0,mu[MAXN];

bool isprime[MAXN];

void init(int n)

{

	for(int i = 2;i <= n;++i)isprime[i] = true;

	mu[1] = 1;

	for(int i = 2;i <= n;++i)

	{

		if(isprime[i])

		{

			prime[++tot] = i;

			mu[i] = -1;

		}

		for(int j = 1;j <= tot && i * prime[j] <= n;++j)

		{

			isprime[i * prime[j]] = false;

			if(i % prime[j] == 0)

			{

				mu[i * prime[j]] = 0;

				break;

			}

			else

			{

				mu[i * prime[j]] = -1 * mu[i];

			}

		}

	}

	return;

}

typedef long long ll;

ll powd[MAXN],sum[MAXN];

ll power(ll a,ll b)

{

	ll res = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a % MOD;

		a = a * a % MOD;

		b = b >> 1;

	}

	return res;

}

int main()

{

	scanf("%d%d",&n,&m);

	init(min(n,m));

	if(n > m)swap(n,m);

	for(int i = 1;i <= m;++i)powd[i] = 1;

	ll ans = 0;

	for(int d = 1;d <= n;++d)

	{

		ll res = 0;

		for(int i = 1;i <= m / d;++i)

		{

			powd[i] = powd[i] * i % MOD;sum[i] = (sum[i - 1] + powd[i]) % MOD;

		}

		for(int i = 1;i <= n / d;++i)

		{

			res = (res + mu[i] * powd[i] % MOD * powd[i] % MOD * sum[n / i / d] % MOD * sum[m / i / d] % MOD + MOD) % MOD;

		}

		ans = (ans + res * power(d,d) % MOD) % MOD;

	}

	cout << ans << endl;

	return 0;

}