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

题目描述

传送门

题目大意：
∑ni=1∑mj=1gcd(i,j)k mod 109+7

题解

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

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

∑d=1ndk∑i=1n∑j=1m[gcd(i,j)=d]

设f(d)=∑ni=1∑mj=1[gcd(i,j)=d] 表示满足gcd(i,j)=d且1<=i<=n,1<=j<=m的(i,j)的对数
F(d)=∑d|nf(n) 表示满足d|gcd(i,j)且1<=i<=n,1<=j<=m的(i,j)的对数
显然F(d)=⌊nd⌋∗⌊md⌋
利用反演公式(2)，可以得到

∑d=1ndk∑i=1nμ(i)∗F(d∗i)

注意F(d)=∑d|nf(n) 中n的范围是无穷，但是如果n超过给出的min(n,m)，那么F(d)就会是0
所以对于枚举的d来说，最大的倍数只能是min(n,m).

∑d=1ndk∑i=1nμ(i)∗⌊nd∗i⌋∗⌊md∗i⌋

设T=d∗i,然后我们调整枚举的顺序

∑T=1n⌊nT⌋∗⌊mT⌋∑d|Tμ(Td)∗dk

⌊nT⌋∗⌊mT⌋ 可以在O(n√+m−−√)的时间内求解，那么问题的关键就是如果快速的计算h(d)=∑d|Tμ(Td)∗dk 发现函数h是狄利克雷卷积的形式。μ与dk都是积性函数(a,b互质,f(ab)=f(a)∗f(b)),所以h也是积性函数。
所以当a,b互质的时候，我们可以利用积性函数的性质h(ab)=h(a)∗h(b)
关键是a,b不互质的时候怎么求解。
h(pixi)=μ(1)∗(pixi)k+μ(pi)∗(pixi−1)k 剩下的项μ都是0
那么对于一个质因数来书,h(pixi) 变成h(pixi+1)的影响就是在h(pixi)的基础上乘pik，其他质因数与其互质，直接乘起来就好啦。
那么对于函数h我们就可用线性筛进行预处理，利用的时候直接O(1)

代码

#include<iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#define N 5000000#define p 1000000007#define LL long longusing 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);    }}