【LG-P4449】于神之怒加强版

P4449 于神之怒加强版

给定 $n,m,k$，计算

${\sum }_{i=1}^n{\sum }_{j=1}^m\gcd (i,j{)}^k$ 对 $10^9+7$ 取模的结果。

---

推柿子：

假设 $n\le m$

${\sum }_{i=1}^n{\sum }_{j=1}^m\gcd (i,j{)}^k$

$={\sum }_{d=1}^n{d}^k{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }[\gcd (i,j)=1]$

$={\sum }_{d=1}^n{d}^k{\sum }_{u=1}^{\lfloor \frac{n}{d}\rfloor }\mu (u)\ast \lfloor \frac{n}{du}\rfloor \ast \lfloor \frac{m}{du}\rfloor$

设 $T=du$，

$=\boxed{\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \ast \lfloor \frac{m}{T}\rfloor \sum _{d|T}^T{d}^k\ast \mu (\frac{T}{d})}$

以上都是基本操作。

---

易知，框起来的柿子中，前半部分可以直接用整除分块做，所以重点就是求后半部分。

设 $f(x)={\sum }_{d|x}^x{d}^k\ast \mu (\frac{x}{d})$。

显然，当 $k=1$ 时它就是个卷积。易得：$f=I{d}_k\ast \mu$。

根据积性函数的性质可得：在 $I{d}_k$ 和 $\mu$ 都是积性函数的情况下，$f$ 也必然是个积性函数。

所以：$f(x)=\prod f(p_i^{a_i})$

即可以线性求 $f(x)$。

---

考虑求 $f(p_i^{a_i})$。

$f(p_i^{a_i})={\sum }_{d|p_i^{a_i}}^x{d}^k\ast \mu (\frac{p_i^{a_i}}{d})$

根据莫比乌斯函数的性质，有：当且仅当 $d=p_i^{a_i-1}$ 或 $d=p_i^{a_i}$ 时 $\mu (\frac{p_i^{a_i}}{d})$ 的值不为 0。

所以：$f(p_i^{a_i})=p_{i}{k* a_i} - p_{i}^{k* (a_i-1)} = (p_i^k-1)* p_{i}^{k* (a_i-1)} $

• 当 $a_i>1$ 时，也有：$f(p_i^{{a_i-1})=p_{i}}{k* (a_i-1)} - p_{i}^{k* (a_i-2)} = (p_i^k-1)* p_{i}^{k* (a_i-2)} $

  综上，得：$\boxed{f(p_i^{a_i})=f(p_i^{a_i-1})\ast p_i^k(a_i>1)}$

• 当 $a_i=1$ 时，$\boxed{f(p_i^{a_i})=p_j^k-1}$

---

推到这里之后就可以根据积性函数的性质线性筛了。

CODE：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;

#define int long long
#define maxm 5000005
#define maxn 5000000
#define mod 1000000007

int t, n, m, k;
int prm[maxm], tot, f[maxm], sum[maxm], low[maxn];
bool vis[maxm];

inline int power(int x, int p)
{
	int res = 1;
	while(p > 0)
	{
		if(p & 1) res = res * x % mod;
		p /= 2, x = x * x % mod;
	}
	return res;
}

inline void pre()
{
	vis[1] = f[1] = sum[1] = 1;
	for(int i = 2; i <= maxn; ++i)
	{
		if(!vis[i])
			f[i] = (-1 + power(i, k) + mod) % mod,
			prm[++tot] = i;
		for(int j = 1; j <= tot and i * prm[j] <= maxn; ++j)
		{
			vis[i * prm[j]] = 1;
			if(i % prm[j] == 0)
			{
				f[i * prm[j]] = f[i] * power(prm[j], k) % mod;
				break;
			}
			f[i * prm[j]] = f[i] * f[prm[j]] % mod;
		}
		sum[i] = sum[i - 1] + f[i];
	}
}

inline int slv(int a, int b)
{
	int ans = 0;
	for(int l = 1, r; l <= a; l = r + 1)
	{
		r = min(a / (a / l), b / (b / l));
		ans = (ans + (sum[r] - sum[l - 1] < 0 ? sum[r] - sum[l - 1] + mod : sum[r] - sum[l - 1]) * (a / l) % mod * (b / l) % mod) % mod;	
	}
	return ans;
}

signed main()
{
	scanf("%lld %lld", &t, &k);
	pre();
	while(t--)
	{
		scanf("%lld %lld", &n, &m);
		if(n > m) swap(n, m);
		printf("%lld\n", slv(n, m));
	}
	return 0;
}

---

——$End$——