BZOJ 4407: 于神之怒加强版|莫比乌斯反演

不会搞数学公式很苦恼！！
flag：会写数学公式之后一定好好写一发题解
非常感谢龙爷（sd第一男选手！！可惜神犇都不写blog）提供线性筛做法
2.16————————————————-
一下均设$n<=m$

$$Ans=\sum_{i=1}^{n}\sum_{j=1}^{m}gcd(i,j)^k$$
设 $f(d)$为 $gcd(x,y)=d$的 $(x,y)$数对的数量
$$g(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i*d)={\lfloor\frac{n}{d}\rfloor}*\lfloor\frac{m}{d}\rfloor$$
$$f(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}u(i)*g(i*d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}u(i)*\lfloor\frac{m}{d*i}\rfloor*\lfloor\frac{n}{d*i}\rfloor$$
$$Ans=\sum_{d=1}^{n}d^k\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}u(i)*\lfloor\frac{m}{d*i}\rfloor*\lfloor\frac{n}{d*i}\rfloor$$
此时两次分块可以做到单次讯问 $O(n)$但是还不够
继续等价变换令 $T=d*i$
$$Ans=\sum_{T=1}^{n}\lfloor\frac{n}{T}\rfloor*\lfloor\frac{m}{T}\rfloor\sum_{d|T}d^k*u(\frac{T}{d})$$
$$h(T)=\sum_{d|T}d^k*u(\frac{T}{d})$$
只需要搞出 $h(T)$的前缀和再用分块处理就可以做到单次讯问 $\sqrt{n}$的复杂度
然后就是鬼畜的线性筛，还要感谢龙爷

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<queue>
#include<vector>
#include<set>
#include<map>
#include<iostream>
#include<algorithm>
#define N 5000005
#define ll long long
#define R 1000000007ll
using namespace std;
int sc()
{
    int i=0,f=1; char c=getchar();
    while(c>'9'||c<'0'){if(f=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9')i=i*10+c-'0',c=getchar();
    return i;
}
ll f[N],K;
int prime[N],low[N];
int a[N],n;
ll cal(ll x)
{
    ll ans=1,y=K;
    while(y)
    {
        if(y&1)ans=ans*x%R;
        x=x*x%R;
        y>>=1;
    }
    return (ans-1+R)%R;
}
void pre()
{
    int top=0; f[1]=low[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!a[i])
        {
            low[i]=prime[++top]=i;
            f[i]=cal(i);
        }
        for(int j=1;i*prime[j]<N;j++)
        {
            a[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                low[i*prime[j]]=low[i]*prime[j];
                if(low[i]==i) 
                    f[i*prime[j]]=f[i]*(f[prime[j]]+1)%R;
                else
                    f[i*prime[j]]=f[i/low[i]]*f[low[i]*prime[j]]%R;
                break;
            }
            low[i*prime[j]]=prime[j];
            f[i*prime[j]]=f[i]*f[prime[j]]%R;
        }
    }
    for(int i=2;i<N;i++)f[i]=(f[i]+f[i-1])%R;
}
int main()
{
    n=sc(),K=sc();
    pre();
    while(n--)
    {
        ll ans=0;
        int n=sc(),m=sc();
        if(n>m)swap(n,m);
        for(int i=1,last;i<=n;i=last+1)
        {
            last=min(n/(n/i),m/(m/i));
            ans=(ans+((ll)(n/i)*(m/i)%R)*(R+f[last]-f[i-1])%R)%R;
        }
        printf("%lld\n",ans);
    }
    return 0;
}