bzoj4407: 于神之怒加强版

链接

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

题解

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

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

$$ans=\sum_{x=1}^nx^k\sum_{i=1}^{\lfloor n/x\rfloor}\sum_{j=1}^{\lfloor m/x\rfloor}[gcd(i,j)=1]$$
　　根据经验，我们得知
　　 $$ans=\sum_{x=1}^nx^k\sum_{i=1}^{\lfloor n/x\rfloor}\mu(i)\lfloor \frac{n}{ix}\rfloor\lfloor \frac{m}{ix}\rfloor$$
　　然后我就写了一个SB暴力 $O(Tn^{\frac{3}{4}})$交上去T了一发。
　　根据题解，我们学会
　　 $$ans=\sum_{x=1}^nx^k\sum_{x|d}^n\mu(\frac{d}{x})\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\\ 　　=\sum_{d=1}^{n}\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\sum_{x|d}x^k\mu(\frac{d}{x})$$
　　令 $g(x)=x^k$，显然 $g(x)$是完全积性函数。
　　因此右面的那一坨就是一个狄利克雷卷积。
　　令 $f(d)=(g*\mu)(d)$
　　根据狄利克雷卷积的运算性质， $f(d)$也是积性函数。
　　所以 $f(d)$就可以线性筛，然后处理出前缀和之后，每次询问就变成 $O(\sqrt N)$的了。
　　那么问题来了，怎么线性筛？框架请参考国家集训队论文，任之洲的《积性函数求和的几种方法》。这里只介绍一些不好实现的细节。
　　对于一个质数，可以直接算出其 $f(x)=x^k-1$，这个是根据 $dirchlet$卷积的定义。对于能拆分成不同质数的数，直接搞即可。
　　对于一个只有一种素数因子的数，以由4推出8为例。
　　 $f(4)=g(1)\mu(4)+g(2)\mu(2)+g(4)\mu(1)$
　　 $f(8)=g(1)\mu(8)+g(2)\mu(4)+g(4)\mu(2)+g(8)\mu(1)$
　　由于 $g$是完全积性函数，因此不必考虑互质的情况，直接给$f(4)$乘 $g(2)$，得到
　　 $f(4)g(2)=g(2)\mu(4)+g(4)\mu(2)+g(8)\mu(1)$
　　离成功好像很接近了，还缺个 $g(1)f(8)$这个怎么办？显然 $f(8)=0$，所以这一项等于0。因此乘上 $g(2)$之后我们就得到 $f(8)$了。
　　关于复杂度：显然只需要用到素数的快速幂，而素数的个数是 $O(\frac{n}{logn})$的，所以算法的复杂度是 $O(N)$。
　　总的复杂度 $O(N+T\sqrt N)$

代码

//莫比乌斯反演
#include <cstdio>
#include <algorithm>
#define maxn 5000010
#define mod 1000000007
#define ll long long
using 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;
}