洛谷 P4449 于神之怒加强版（莫比乌斯反演+积性函数线性筛）

传送门

---

题目大意

给定$n,m,k$，计算${\sum }_{i=1}^n{\sum }_{j=1}^{m}gcd(i,j{)}^k$

反演部分

方便起见设$n\le m,gcd(i,j)\iff (i,j)$

常用的转化思路，枚举$gcd$：
${\sum }_{i=1}^n{\sum }_{j=1}^{m}gcd(i,j{)}^k={\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{d=1}^n{d}^k[gcd(i,j)=d]$

然后将后面部分提到前面：

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

根据$gcd(i,j)=d\iff gcd(\lfloor \frac{i}{d}\rfloor ,\lfloor \frac{i}{d}\rfloor )=1$，并转化到两个$\sum$上，得：

${\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]$

然后根据$[gcd(i,j)=1]\iff [\frac{1}{n}]\iff {\sum }_{d|n}\mu (d)$，得：

${\sum }_{d=1}^n{d}^k{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }{\sum }_{t|(i,j)}\mu (t)$

然后是关键的一步，改为枚举$t$，然后联系整除分块得：

${\sum }_{d=1}^n{d}^k{\sum }_{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t){\sum }_{i=1,t|i}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1,t|j}^{\lfloor \frac{m}{d}\rfloor }1={\sum }_{d=1}^n{d}^k{\sum }_{t=1}^{\lfloor \frac{n}{d}\rfloor }\mu (t)\lfloor \frac{n}{dt}\rfloor \lfloor \frac{m}{dt}\rfloor$

令$T=dt$，然后改为枚举$T$，最后化简为${\sum }_{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor {\sum }_{d|T}{d}^k\mu (\frac{T}{d})$

线性筛部分

对于上述式子的最后一部分：$f(n)={\sum }_{d|T}{d}^k\mu (\frac{T}{d})$，因为幂函数$i{d}_k(n)$和莫比乌斯函数$\mu (n)$都是积性函数，根据狄利克雷卷积，那么二者的乘积$f=i{d}_k\ast \mu$也是一个积性函数，积性函数可以考虑欧拉筛埃氏筛，对于$5e6$的数据范围，$nlogn$的筛法会$TLE$，只能想办法线性筛

蒟蒻第一次筛奇怪神奇的积性函数，参考了这篇博客

• 当$n=p,p\in primes$，$f(p)=1^k\ast \mu (p)+p^k\ast \mu (1)=p^k-1$
• 当$n=p^x,p\in primes$，注意到莫比乌斯函数只要含有平方项以上的素数就为0，那么$f(p)=(p^k{)}^{x-1}\ast \mu (p)+(p^k{)}^x\ast \mu (1)=(p^k{)}^{x-1}(p^k-1)$
  再考虑$f(p^{x+1})$，显然：$f(p^{x+1})=(p^k{)}^x(p^k-1)=f(p^x)p^k$

有了上面的结论，设$x,y;x\in Z_{+},y\in primes$，考虑如何线性筛：

• 当$y|x$时，取$x$中$y$的幂次$c$，由$f$的积性有$f(xy)=f(\frac{xy}{y^{c+1}})f(y^{c+1})$，又因为$f(p^{c+1})=f(p^c)p^k$，得$f(xy)=f(\frac{x}{y^c})f(y^c)y^k=f(x)y^k$
• 当$y\nmid x$时，由积性函数的性质得$f(xy)=f(x)f(y)$

代码

//
// Created by Happig on 2020/9/25
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>

using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
const int Mod = 1e9 + 7;
const int maxn = 5e6 + 10;

vector<int> prime;
bitset<maxn> vis;
ll f[maxn], sum[maxn];
int k;

ll qkp(ll x, ll n, ll p) {
    x %= p;
    ll ans = 1;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

void init() {
    memset(sum, 0, sizeof sum);
    f[1] = 1;
    vis.reset();
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) {
            prime.push_back(i);
            f[i] = (qkp(i, k, Mod) - 1 + Mod) % Mod;
        }
        for (int j = 0; j < prime.size() && i * prime[j] < maxn; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j]) {
                f[i * prime[j]] = f[i] * f[prime[j]] % Mod;
            } else {
                f[i * prime[j]] = f[i] * qkp(prime[j], k, Mod) % Mod;
                break;
            }
        }
    }
    for (int i = 1; i < maxn; i++) {
        sum[i] = (sum[i - 1] + f[i]) % Mod;
    }
}

int main() {
    //freopen("in.txt","r",stdin);
    //freopen("out.txt", "w", stdout);
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    int t, n, m;
    cin >> t >> k;
    init();
    while (t--) {
        cin >> n >> m;
        ll up = min(n, m), ans = 0;
        for (ll l = 1, r; l <= up; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            if (r > up) r = up;
            ans += (sum[r] - sum[l - 1] + Mod) % Mod * (n / l) % Mod * (m / l) % Mod;
            ans %= Mod;
        }
        cout << ans << ENDL;
    }
    return 0;
}