【BZOJ4407】于神之怒加强版

Description

给下N,M,K.求



Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。

Output

如题

Sample Input

1 2

3 3

Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

Source

命题人：成都七中张耀楠，鸣谢excited上传。

300年没写过反演的题,今天再做一道内心几乎是崩溃的(摔

感谢鸟神提供的答疑w

题目要求求∑i=1n∑j=1mgcd(i,j)k

即为∑d=1min(n,m)dk∑i=1n∑j=1m[gcd(i,j)=d]

后面那个东西很眼熟?似乎在梦中见过的样子w

∑i=1n∑j=1m[gcd(i,j)=d]

=∑i=1⌊nd⌋∑j=1⌊md⌋[gcd(i,j)=1]

显然可以反演咯w

此时的原式∑d=1min(n,m)dk∑i=1⌊nd⌋∑j=1⌊md⌋[gcd(i,j)=1]

=∑d=1min(n,m)dk∑x=1min(n,m)μ(x)⌊nxd⌋⌊mxd⌋

化到这里求解的复杂度是O(Tn),过不了.

怎么降复杂度? 丢掉d!

令y=xd,现在d是y的约数了.取代一下.

原式

∑d=1min(n,m)dk∑x=1min(n,m)μ(yd)⌊ny⌋⌊my⌋

=∑y=1min(n,m)⌊ny⌋⌊my⌋∑d|ymin(n,m)μ(yd)dk

里面那个是函数f(d)=dk与莫比乌斯函数的卷积.有前缀和就行了.

dk这个可以放线筛里求,质数的话暴力快速幂,非质数的话就可以用两个因子的值搞出来.这样是O(n+Tn√)的.

如果直接对每个数暴力求得话是O(nlogn+Tn√)的.(跑50+s的那些都是这种)

最后我竟然跪在不会线性的求这个前缀和上了(摔

蓝月亮跑rank1好厉害呀

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define MAXN 5000010
#define GET (ch>='0'&&ch<='9')
#define P 1000000007
#define LL long long
using namespace std;
int T,k,n[2010],m[2010],maxn;
bool not_prime[MAXN];
int prime[MAXN],top;
int mu[MAXN],fac[MAXN];
LL f[MAXN],F[MAXN];
void in(int &x)
{
char ch=getchar();x=0;
while (!GET)    ch=getchar();
while (GET) x=x*10+ch-'0',ch=getchar();
}
inline LL Pow(LL a,LL b)
{
a%=P;LL ret=1;
for (;b;b>>=1,a=a*a%P)  if (b&1)    ret*=a,ret%=P;
return (ret+P)%P;
}
void check()
{
mu[1]=1;F[1]=1;
for (int i=2;i<=maxn;i++)
{
if (!not_prime[i])  prime[++top]=i,fac[i]=i,mu[i]=-1,f[i]=Pow(i,k),F[i]=f[i]-1;
for (int j=1;j<=top&&i*prime[j]<=maxn;j++)
{
not_prime[i*prime[j]]=1;mu[i*prime[j]]=-mu[i];
if (i%prime[j]==0)
{
mu[i*prime[j]]=0;fac[i*prime[j]]=fac[i]*prime[j];
F[i*prime[j]]=fac[i]!=i?F[i/fac[i]]*F[fac[i]*prime[j]]%P:F[i]*f[prime[j]]%P;
break;
}
F[i*prime[j]]=F[i]*F[prime[j]]%P;fac[i*prime[j]]=prime[j];
}
}
for (int i=1;i<=maxn;i++)   mu[i]+=mu[i-1],F[i]+=F[i-1],F[i]%=P;
}
int main()
{
int i;
for (in(T),in(k),i=T;i;i--) in(n[i]),in(m[i]),maxn=max(maxn,max(n[i],m[i]));
for (check();T;T--)
{
int N=n[T],M=m[T];
LL ans=0;int t=min(N,M),last=1;
for (int i=1;i<=t;i=last+1)
{
last=min(N/(N/i),M/(M/i));
ans+=(F[last]-F[i-1]+P)*(N/i)%P*(M/i)%P;ans%=P;
}
printf("%lld\n",ans);
}
}