bzoj4407: 于神之怒加强版

链接

　　http://www.lydsy.com/JudgeOnline/problem.php?id=4407

题解

夜里挑灯做题，梦回莫比乌斯。

　　所以这是一道莫比乌斯反演。
　　设n<m
　　根据题意，我们写出式子
　　

ans=∑x=1nxk∑i=1⌊n/x⌋∑j=1⌊m/x⌋[gcd(i,j)=1]

　　根据经验，我们得知
　　

ans=∑x=1nxk∑i=1⌊n/x⌋μ(i)⌊nix⌋⌊mix⌋

　　然后我就写了一个SB暴力O(Tn34)交上去T了一发。
　　根据题解，我们学会
　　

ans=∑x=1nxk∑x|dnμ(dx)⌊nd⌋⌊md⌋　　=∑d=1n⌊nd⌋⌊md⌋∑x|dxkμ(dx)

　　令g(x)=xk，显然g(x)是完全积性函数。
　　因此右面的那一坨就是一个狄利克雷卷积。
　　令f(d)=(g∗μ)(d)
　　根据狄利克雷卷积的运算性质，f(d)也是积性函数。
　　所以f(d)就可以线性筛，然后处理出前缀和之后，每次询问就变成O(N−−√)的了。
　　那么问题来了，怎么线性筛？框架请参考国家集训队论文，任之洲的《积性函数求和的几种方法》。这里只介绍一些不好实现的细节。
　　对于一个质数，可以直接算出其f(x)=xk−1，这个是根据dirchlet卷积的定义。对于能拆分成不同质数的数，直接搞即可。
　　对于一个只有一种素数因子的数，以由4推出8为例。
　　f(4)=g(1)μ(4)+g(2)μ(2)+g(4)μ(1)
　　f(8)=g(1)μ(8)+g(2)μ(4)+g(4)μ(2)+g(8)μ(1)
　　由于g是完全积性函数，因此不必考虑互质的情况，直接给f(4)乘g(2)，得到
　　f(4)g(2)=g(2)μ(4)+g(4)μ(2)+g(8)μ(1)
　　离成功好像很接近了，还缺个g(1)f(8)这个怎么办？显然f(8)=0，所以这一项等于0。因此乘上g(2)之后我们就得到f(8)了。
　　关于复杂度：显然只需要用到素数的快速幂，而素数的个数是O(nlogn)的，所以算法的复杂度是O(N)。
　　总的复杂度O(N+TN−−√)

代码

//莫比乌斯反演#include <cstdio>#include <algorithm>#define maxn 5000010#define mod 1000000007#define ll long longusing namespace std;int x[maxn], K, prime[maxn], mark[maxn], tmp[maxn], f[maxn];inline int fastpow(int a, int b, int p){    int ans=1, t=a;    for(;b;b>>=1,t=(ll)t*t%p)if(b&1)ans=(ll)ans*t%p;    return ans;}void init(){    int i, j, t;    f[1]=1;    for(i=2;i<maxn;i++)    {        if(!mark[i])        {            prime[++prime[0]]=i;            tmp[i]=i;            x[i]=fastpow(i,K,mod);            f[i]=(x[i]-1+mod)%mod;        }        for(j=1;j<=prime[0] and i*prime[j]<maxn;j++)        {            t=i*prime[j];            mark[t]=1;            if(i%prime[j]==0)            {                tmp[t]=tmp[i]*prime[j];                if(tmp[t]!=t)f[t]=(ll)f[t/tmp[t]]*f[tmp[t]]%mod;                else f[t]=(ll)f[i]*x[prime[j]]%mod;                break;            }            f[t]=(ll)f[i]*f[prime[j]]%mod;            tmp[i*prime[j]]=prime[j];        }    }    for(i=2;i<maxn;i++)f[i]=(f[i]+f[i-1])%mod;}inline void calc(int n, int m){    ll ans=0, lim=1e17;    int i, last;    if(n>m)swap(n,m);    for(i=1;i<=n;i=last+1)    {        last=min(n/(n/i),m/(m/i));        ans+=(ll)(n/i)*(m/i)%mod*(f[last]-f[i-1])%mod;        if(ans>lim)ans%=mod;    }    printf("%d\n",(int)(ans%mod+mod)%mod);}int main(){    int T, N, M;    scanf("%d%d",&T,&K);    for(init();T;T--)    {        scanf("%d%d",&N,&M);        calc(N,M);    }    return 0;}