bzoj 4407: 于神之怒加强版 （反演+线性筛）

题目描述

传送门

题目大意：
$\sum_{i=1}^n\sum_{j=1}^m gcd(i,j)^k \space mod \space10^9+7$

题解

一道非常不错的莫比乌斯反演，想了一晚啊。。。。
首先引入反演的两个公式
（1）如果$F(n)=\sum\limits_{d|n} f(n)$, 那么$f(n)=\sum\limits_{d|n}\mu(d)F({n\over d})$
（2）如果$F(n)=\sum\limits_{n|d} f(d)$,那么$f(n)=\sum\limits_{n|d}\mu({d\over n})*F(d)$
那么反演有什么用呢？可以实现F(n),f(n)之间的转换求解，使不易求的式子转化成易求的形式。
最常见的应用就是$[n=1]=\sum\limits_{d|n}\mu(d)$

那么对于这道题来说我们先对式子进行改动,设$n<=m$

$$\sum\limits_{d=1}^n d^k \sum_{i=1}^n\sum_{j=1}^m [gcd(i,j)=d]$$
设 $f(d)= \sum_{i=1}^n\sum_{j=1}^m [gcd(i,j)=d]$ 表示满足 $gcd(i,j)=d$且 $1<=i<=n,1<=j<=m$的 $(i,j)$的对数
$F(d)=\sum_{d|n} f(n)$ 表示满足 $d|gcd(i,j)$且 $1<=i<=n,1<=j<=m$的 $(i,j)$的对数
显然 $F(d)={\lfloor {n\over d}\rfloor}*{\lfloor {m\over d}\rfloor}$
利用反演公式(2)，可以得到
$$\sum\limits_{d=1}^n d^k \sum_{i=1}^n \mu(i)*F(d*i)$$
注意 $F(d)=\sum_{d|n} f(n)$ 中n的范围是无穷，但是如果n超过给出的min(n,m)，那么F(d)就会是0
所以对于枚举的d来说，最大的倍数只能是min(n,m).
$$\sum\limits_{d=1}^n d^k \sum_{i=1}^n \mu(i)*{\lfloor {n\over d*i}\rfloor}*{\lfloor {m\over d*i}\rfloor}$$
设 $T=d*i$,然后我们调整枚举的顺序
$$\sum_{T=1}^n {\lfloor {n\over T}\rfloor}*{\lfloor {m\over T}\rfloor} \sum\limits_{d|T} \mu({T\over d})*d^k$$
${\lfloor {n\over T}\rfloor}*{\lfloor {m\over T}\rfloor}$ 可以在 $O(\sqrt n+\sqrt m)$的时间内求解，那么问题的关键就是如果快速的计算 $h(d)=\sum\limits_{d|T} \mu({T\over d})*d^k$ 发现函数 $h$是狄利克雷卷积的形式。$\mu$与 $d^k$都是积性函数(a,b互质, $f(ab)=f(a)*f(b)$),所以h也是积性函数。
所以当a,b互质的时候，我们可以利用积性函数的性质 $h(ab)=h(a)*h(b)$
关键是a,b不互质的时候怎么求解。
$h(pi^{xi})=\mu(1)*(pi^{xi})^k+\mu(pi)*(pi^{xi-1})^k$ 剩下的项 $\mu$都是0
那么对于一个质因数来书, $h(pi^{xi})$ 变成 $h(pi^{xi+1})$的影响就是在 $h(pi^{xi})$的基础上乘 $pi^k$，其他质因数与其互质，直接乘起来就好啦。
那么对于函数h我们就可用线性筛进行预处理，利用的时候直接 $O(1)$

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 5000000
#define p 1000000007
#define LL long long
using namespace std;
LL f[N+3],g[N+3];
int n,m,k,T,pd[N+3],prime[N+3];
LL quickpow(LL num,int x)
{
    LL ans=1; LL base=num%p;
    while (x) {
        if (x&1) ans=base*ans%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
void init()
{
    f[1]=1;
    for (int i=2;i<=N;i++) {
        if (!pd[i]) {
            prime[++prime[0]]=i;
            g[prime[0]]=quickpow(i,k);
            f[i]=(g[prime[0]]-1+p)%p;
        }
        for (int j=1;j<=prime[0];j++) {
            if (i*prime[j]>N) break;
            pd[i*prime[j]]=1;
            if (i%prime[j]==0) {
                f[i*prime[j]]=f[i]*g[j]%p;
                break;
            }
            else f[i*prime[j]]=f[i]*f[prime[j]]%p;
        }
    } 
    for (int i=1;i<=N;i++) f[i]=(f[i-1]+f[i])%p;
    //for (int i=1;i<=10;i++) cout<<f[i]<<" ";
    //cout<<endl;
}
int main()
{
    freopen("a.in","r",stdin);
    scanf("%d%d",&T,&k);
    init();
    while (T--) {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
        int j=0; LL ans=0;
        for (int i=1;i<=n;i=j+1) {
            j=min(n/(n/i),m/(m/i));
            j=min(j,n);
            LL t=(LL)(n/i)*(m/i)%p;
            ans=(ans+(LL)t*(f[j]-f[i-1]+p)%p)%p;
        }
        printf("%I64d\n",(ans%p+p)%p);
    }
}