BZOJ 3309 莫比乌斯反演

题目链接

---

题意：多组数据，已知$a,b(<=1e7)$，定义$f(n)$为$n$所含质因子的最大幂指数。
求$$\sum _{i=1}^a\sum _{j=1}^{b}f(gcd(i,j))$$

此前写过一个$f(n)$为$n$的约数和的题解：我是链接= =

感觉此类题的推导都是一个套路，定义最终答案为$Ans$，首先转换枚举变量为$gcd$，并设$g(i):gcd==i$的方案数，则：
$$Ans=\sum _{i=1}^{min(a,b)}f(i)g(i)$$

然后利用莫比乌斯反演，设：
$G(i):gcd为i的倍数的方案数$，则：
$$G(i)=\lfloor \frac{a}{i}\rfloor \lfloor \frac{b}{i}\rfloor$$
$$G(i)=\sum _{i|d}g(d)\Rightarrow g(i)=\sum _{i|d}\mu (\frac{d}{i})G(d)$$

将$g(i)$代入$Ans$表达式，得出：
$$Ans=\sum _{i=1}^{min(a,b)}f(i)\sum _{i|d}\mu (\frac{d}{i})\lfloor \frac{a}{d}\rfloor \lfloor \frac{b}{d}\rfloor$$

转化枚举变量为$d$，得出：
$$Ans=\sum _{d=1}^{min(a,b)}\lfloor \frac{a}{d}\rfloor \lfloor \frac{b}{d}\rfloor \sum _{i|d}f(i)\mu (\frac{d}{i})$$

---

设$g(x)$(此处的$g(x)$与上面的$g(x)$的定义完全不同)：
$$g(x)=\sum _{i|x}f(i)\mu (\frac{x}{i})$$
现在的难点便在于如何利用线性筛快速得出$g(x)$然后再求前缀和了，如果能实现的话我们就可以在$O(n)$复杂度下预处理，$O(T\sqrt{n})$的复杂度下快速求解了。

此处参考了PoPoQQQ的博客
讲得挺好的，但确实很难理解，自己在纸上模拟了几个例子才想通。
数论真是玄学啊qwq

---

代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long ll;

const int A = 1e7 + 10;
bool vis[A];
int pri[A],cnt[A],tot;
ll g[A],val[A]; //假设p为i的最小质因子，则cnt[i]:p的次数 val[i]:p^cnt[i]的值

void init(){
    tot = 0;
    for(int i=2 ;i<A ;i++){
        if(!vis[i]){pri[++tot] = i;cnt[i]=1;val[i] = i;g[i] = 1;}
        for(int j=1 ;j<=tot && i*pri[j]<A ;j++){
            vis[i*pri[j]] = 1;
            if(i%pri[j] == 0){
                cnt[i*pri[j]] = cnt[i] + 1;
                val[i*pri[j]] = val[i]*pri[j];
                int now = i/val[i];
                if(now == 1) g[i*pri[j]] = 1;
                else         g[i*pri[j]] = (cnt[now] == cnt[i*pri[j]]?-g[now]:0);
                break;
            }
            cnt[i*pri[j]] = 1;
            val[i*pri[j]] = pri[j];
            g[i*pri[j]] = (cnt[i] == cnt[i*pri[j]]?-g[i]:0);
        }
    }
    for(int i=1 ;i<A ;i++) g[i] += g[i-1];
}

int main(){
    init();
    int T;
    scanf("%d",&T);
    while(T--){
        int a,b,last;
        scanf("%d%d",&a,&b);
        if(a>b) swap(a,b);
        ll ans = 0;
        for(int i=1 ;i<=a ;i=last+1){
            last = min(a/(a/i),b/(b/i));
            ans += 1LL*(a/i)*(b/i)*(g[last]-g[i-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
}