[BZOJ4407]于神之怒加强版（莫比乌斯反演）

题目描述

传送门

题解

感觉这题非常强
假设$n<m$

$$\sum\limits_{i=1}^n\sum\limits_{j=1}^m(i,j)^k$$
$$=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d=1}^n[(i,j)=d]d^k$$
$$=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^{n\over d}\sum\limits_{j=1}^{m\over d}[(i,j)=1]$$
用一下反演公式 $[n=1]=\sum\limits_{d|n}\mu(d)$
$$=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^{n\over d}\sum\limits_{j=1}^{m\over d}\sum\limits_{t|(i,j)}\mu(t)$$
$$=\sum\limits_{d=1}^nd^k\sum\limits_{t=1}^{n\over d}\sum\limits_{i=1}^{n\over d}[t|i]\sum\limits_{j=1}^{m\over d}[t|j]\mu(t)$$
$$=\sum\limits_{d=1}^nd^k\sum\limits_{t=1}^{n\over d}{\lfloor{n\over dt}\rfloor}{\lfloor{m\over dt}\rfloor}\mu(t)$$
令 $T=dt$，此时 $d$为$T$的约数
$$=\sum\limits_{T=1}^n{\lfloor{n\over T}\rfloor}{\lfloor{m\over T}\rfloor}\sum\limits_{d|T}\mu({T\over d})d^k$$
发现最后的 $\sum\limits_{d|T}\mu({T\over d})d^k$就是狄利克雷卷积的形式
令 $g(n)=n^k,f(n)=\sum\limits_{d|n}\mu({n\over d})g(d)$
$\mu$和 $g$都是积性函数，所以$f=\mu\times g$也是一个积性函数
那么可以用线性筛求 $f$
当$a$为质数时显然 $f(a)=a^k-1$
当 $(a,b)=1$时 $f(ab)=f(a)f(b)$
当 $ab$不互质时，也就是用 $f(p_i^{k_i})$和 $p_i$计算 $f(p_i^{k_i+1})$时，利用莫比乌斯函数的性质，发现只有 $\mu(1)$和 $\mu(p_i)$是有价值的
$f(p_i^{k_i})=\mu(1)*(p_i^{k_i})^k+\mu(p_i)*(p_i^{k_i-1})^k$
$f(p_i^{k_i+1})=\mu(1)*(p_i^{k_i+1})^k+\mu(p_i)*(p_i^{k_i})^k$
也就是说 $f(p_i^{k_i+1})=f(p_i^{k_i})*p_i^k$
那么就可以筛筛筛了
计算出来 $f(n)$了之后原式变成了
$$\sum\limits_{T=1}^n{\lfloor{n\over T}\rfloor}{\lfloor{m\over T}\rfloor}f(T)$$
预处理出 $f(n)$的前缀和，只有 $O(\sqrt n+\sqrt m)$种取值
分块求解
时间 $O(n+T(\sqrt n+\sqrt m))$

代码

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define Mod 1000000007
#define LL long long
#define N 5000005

int T,n,m,k;
int p[N],prime[N];
LL ans,f[N];

LL fast_pow(LL a,int p)
{
    LL ans=1;
    for (;p;p>>=1,a=a*a%Mod)
        if (p&1)
            ans=ans*a%Mod;
    return ans;
}
void get(int n)
{
    f[1]=1;
    for (int i=2;i<=n;++i)
    {
        if (!p[i])
        {
            prime[++prime[0]]=i;
            f[i]=(fast_pow(i,k)-1+Mod)%Mod;
        }
        for (int j=1;j<=prime[0]&&i*prime[j]<=n;++j)
        {
            p[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                f[i*prime[j]]=f[i]*fast_pow(prime[j],k)%Mod;
                break;
            }
            else f[i*prime[j]]=f[i]*f[prime[j]]%Mod;
        }
    }
    for (int i=2;i<=n;++i) f[i]=(f[i]+f[i-1])%Mod;
}
int main()
{
    scanf("%d%d",&T,&k);
    get(5000000);
    while (T--)
    {
        scanf("%d%d",&n,&m);
        if (n>m) swap(n,m);
        ans=0;
        for (int i=1,j=0;i<=n;i=j+1)
        {
            j=min(n/(n/i),m/(m/i));
            ans+=(f[j]-f[i-1])*(n/i)%Mod*(m/i)%Mod;
            ans=(ans%Mod+Mod)%Mod;
        }
        printf("%lld\n",ans);
    }
}