P4449 于神之怒加强版[莫比乌斯反演]

· · 题解

P4449 于神之怒加强版[莫比乌斯反演]

题意: 给定i,j,k

求

\sum_{i=1}^n\sum_{j=1}^m gcd^k(i,j)

solution:

一步一步推式子,

我们利用莫比乌斯反演,设

f(n) = n^k f(n) = \sum_{d|n} g(d) g(n) = \sum_{d|n} f(d) \mu(\frac{n}{d})

化简原式

\sum_{i=1}^n\sum_{j=1}^m gcd^k(i,j) =\sum_{i=1}^n\sum_{j=1}^m f(gcd(i,j)) =\sum_{i=1}^n\sum_{j=1}^m\sum_{d|gcd(i,j)} g(d)

由于

d|gcd(i,j) \Leftrightarrow d|i,d|j

所以d必须是i,j的约数

于是我们想到换个思路,枚举i,j的约数,然后考虑约数的贡献

所以我们改变枚举顺序,将d提到前面来(细节:由于 d 为 i,j 的公共约数,故d最大为min(n,m)

\sum_{d=1}^{min(n,m)}\sum_{d|i}^n\sum_{d|j}^m g(d) =\sum_{d=1}^{min(n,m)} \lfloor{\frac{n}{d}}\rfloor \lfloor{\frac{m}{d}}\rfloor g(d)

于是我们只需枚举d便可以统计出所有答案,复杂度为O(n)

考虑 \lfloor\frac{n}{d}\rfloor \lfloor \frac{m}{d}\rfloor 可以用数论分块处理,再预处理一下g函数的前缀和,复杂度降为O(\sqrt{n})

那么麻烦来了,我们如何快速计算出g函数

本题数据范围为5e6,O(nloglogn)是过不了的 (事实上,O(nloglogn)最大能承受的数据案范围大致为1e6)

所以我们的解决办法就是线性筛

由于 g = f * \mu , f 和 \mu都为积性函数,故g也为积性函数

考虑线性筛的过程

当 pri[j] 不为 i 的因子时, pri[j] 与 i 互质,可得 g[pri[j] * i] = g[pri[j]] * g[i]

当 pri[j] | i 时 , 根据线性筛的原理我们可以得出 , pri[j] 一定是 i 的最小质因子 . 我们分

两种情况讨论 , 分类讨论的原因下面会讲 ( 此处设 low[i] = p_1^{a_1} , p_1 为i的最小质因子 )

1. low[i] != i

我们把 i 的最小质因子全部移到 pri[j] 上去 , 这样 i 中剩下的就没有 pri[j] 的因子了,

所以两者互质 , 即

g[pri[j] * i] = g[\frac{i}{low[i]}] * g[low[i] * pri[j]] 2. low[i] == i

换句话说 , 此时的 i 为 pri[j] 的幂 . 对于这种情况 , 我们为什么不能像上面那样求呢?

如果按照上面的步骤 , 我们会发现一个很尴尬的式子 g[p^k] = g[p^k] * g[1]

所以我们不能按照上面的方法求 对于质数幂的积性函数 , 手推一下可以发现下面这个式子(只有质数的幂才可以这样推!) $$g[p^k] = g[p^{k-1}] * f[p] + f[1] * \mu[p^k] (g = \mu * f)$$ 回到题目中来,式子即变成 $$g[pri[j] * i] = g[i] * f[pri[j]] + f[1] * \mu[pri[j] * i]$$ 然后我们的问题就圆满解决了 细节 : 注意每次 $ i $ 筛到一个质数时 , $ g[i] $ 都要更新为$ f[1] * \mu[i] + f[i] * mu[1] $ , 同时 $ low[i] = i i $ 与 $ pri[j] $ 互质的时候 , 更新 $ low $ , $ low[i * pri[j]] = pri[j] g $ 也要更新 , $ g[i * pri[j]] = g[i] * g[pri[j]]

下面上代码

// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define ll long long
#pragma GCC optimize(2)

using namespace std;

inline int read() {
    int res = 0 , flag = 1; char ch = getchar();
    while(!isdigit(ch)) { if(ch == '-') flag = -1; ch = getchar(); }
    while(isdigit(ch)) { res = (res << 3) + (res << 1) + ch - '0'; ch = getchar(); }
    return res * flag;
}

const int N=5e6;
const int mod=1000000007;

int mu[N+5],prime[N+5],ncnt;
ll k;
ll f[N+5],low[N+5];
ll g[N+5];
ll sum[N+5];
bool vis[N+5];

inline ll pow(ll a,ll b) {
    ll ret = 1;
    while(b) {
        if(b & 1) (ret *= a) %= mod;
        (a *= a) %= mod;
        b >>= 1;
    }
    return ret;
}

inline void Mobius() {
    mu[1] = 1;
    g[1] = 1;//积性函数g[1] == 1!!!!!!!!!!!!!

    for(register ll i = 1 ; i <= N ; ++ i) f[i] = pow(i,k);// , printf("%lld\n",f[i]);

    for(register int i = 2 ; i <= N ; ++ i) {
        if(!vis[i]) {
            prime[++ncnt] = i;
            low[i] = i;
            mu[i] = -1;
            g[i] = (f[i] - 1) % mod;
        }
        for(register int j = 1 ; j <= ncnt and i * prime[j] <= N ; ++ j) {
            vis[i * prime[j]] = 1;//筛掉//每次由最小的累积
            if(i % prime[j] == 0) {
//              mu[i * prime[j]] = 0;
                low[i * prime[j]] = low[i] * prime[j];

                if(low[i] == i) {//i为素数的幂形式
                    g[i * prime[j]] = (g[i] * f[prime[j]]) % mod;// + f[1] * mu[i * prime[j]]
                }
                else {
                    g[i * prime[j]] = (g[i / low[i]] * g[low[i] * prime[j]]) % mod;
                }
                break;
            }
            else {
                g[i * prime[j]] = (g[i] * g[prime[j]]) % mod;
                low[i * prime[j]] = prime[j];
                mu[i * prime[j]] = -mu[i];
            }
        }
    }

    for(register int i = 1 ; i <= N ; ++ i) sum[i] = (sum[i-1] + g[i]) % mod;
//  for(register int i = 1 ; i <= N ; ++ i) printf("%lld\n",sum[i]);
}

inline int min(int a , int b) {
    return a < b ? a : b;
}

inline ll work(int n,int m) {
    ll ans = 0;
    for(register int d = 1 , t; d <= min(n,m) ; d = t + 1) {
        t = min(n / (n/d) , m / (m/d));
        ans = (ans + 1ll * (n / d) * (m / d) % mod * ((sum[t] - sum[d - 1]) % mod + mod) % mod) % mod;
    }
    return ans;
}

int main() {
//  freopen("god.in","r",stdin);
//  freopen("god.out","w",stdout);

    int test = read(); k = read();
    Mobius();

    while(test --) {
        int n = read() , m = read();
        printf("%lld\n",work(n,m));
    }
    return 0;
}

end

莫比乌斯反演套路总结

用莫比乌斯反演优化算法的题目一般都涉及到 gcd

优化的原理和关键是 d | gcd(i , j) \Leftrightarrow d | i , d | j , 然后枚举约数 d 算贡献于是我们将 O(n^2) 枚举降为了

O(n) $ 枚举, 简化后式子中的 $ \lfloor \frac{n}{d} \rfloor * \lfloor \frac{m}{d} \rfloor $ 也可以用数论分块优化 , 复杂度再降为 $ O(\sqrt{n})

而运用莫比乌斯反演的关键是将题目中的式子化成如下形式

\sum_{i=1}^n\sum_{j=1}^m f(gcd(i,j))

然后就可以运用上述套路了

线性筛常见数论函数及自定义积性函数也是莫比乌斯反演可以考的一个板块

然而我还不会筛数论函数

以后有时间再写一篇这方面的博客