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

flag：会写数学公式之后一定好好写一发题解

2.16————————————————-

Ans=i=1nj=1mgcd(i,j)k

f(d)$f(d)$gcd(x,y)=d$gcd(x,y)=d$(x,y)$(x,y)$数对的数量
g(d)=i=1ndf(id)=ndmd

f(d)=i=1ndu(i)g(id)=i=1ndu(i)mdindi

Ans=d=1ndki=1ndu(i)mdindi

Ans=T=1nnTmTd|Tdku(Td)

h(T)=d|Tdku(Td)

#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;
}

©️2019 CSDN 皮肤主题: 大白 设计师: CSDN官方博客