[Update] 我好像现在都看不懂我当时在写什么了=-=
(Description)
求(sum_{i=a}^bsum_{j=c}^d[(i,j)=k])
(Solution)
首先是把下界作为1.可以化为求
[sum_{i=1}^{lfloorfrac{N}{k}
floor}sum_{j=1}^{lfloorfrac{M}{k}
floor}[(i,j)=1]
]
说明:大概就我不能直接看出来了。。
首先要求([1,N])中有多少(i,i|k),再求[1,j]中有多少(j,j|k且(i,j)=1),显然这个(i,j)的上界就分别是(lfloorfrac{N}{k}
floor,lfloorfrac{M}{k}
floor),答案就是((i,j)=1)的((i,j))数对个数。
现在考虑如何求上面的式子。
由莫比乌斯反演,有
[F(d)=sum_{d|n}f(n)Leftrightarrow f(d)=sum_{d|n}mu(frac{n}{d})F(n)
]
设(f[i])为满足(gcd(x,y)=i)的((x,y))对数,其中(1leq xleq b,1leq yleq d);
(F[i])为满足(i|gcd(x,y))的((x,y))对数,其中(1leq xleq b,1leq yleq d)。
显然有(F[i]=sum_{i|n}f[n]Rightarrow f[i]=sum_{i|n}mu(frac{n}{i})F[n])
又显然有(F[i]=lfloorfrac{b}{i}
floorlfloorfrac{d}{i}
floor),那么
[f[i]=sum_{i|n,1leq nleq min(b,d)}mu(frac{n}{i})lfloorfrac{b}{n}
floorlfloorfrac{d}{n}
floor
]
令(k=frac{n}{i}),即(n=ki),令(b'=frac{b}{i},d'=frac{d}{i}),则
[f[i]=sum_{k=1}^{min(b',d')}mu(k)lfloorfrac{b'}{k}
floorlfloorfrac{d'}{k}
floor
]
(本题i就是1。)
上面这个式子还是(O(n^2))的。。还是要分块计算。
/*
最后求解要容斥:(a,c为开区间,另外其实并不分左右)
ans=f[a~b,c~d]=f[b,d]-f[a,d]-f[b,c]+f[a,c]
*/
#include<cstdio>
#include<cctype>
#include<algorithm>
#define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
const int N=5e4+3,MAXIN=1<<17;
int cnt,P[N+2],mu[N+2],sum[N+2];
bool Not_P[N+2];
char IN[MAXIN],*SS=IN,*TT=IN;
inline int read()
{
int now=0,f=1;register char c=gc();
for(;!isdigit(c);c=gc()) if(c=='-') f=-1;
for(;isdigit(c);now=now*10+c-'0',c=gc());
return now*f;
}
void Init()
{
mu[1]=1;
for(int i=2;i<N;++i)
{
if(!Not_P[i]) P[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*P[j]<N;++j)
{
Not_P[i*P[j]]=1;
if(!(i%P[j])) {mu[i*P[j]]=0; break;}
mu[i*P[j]]=-mu[i];
}
}
for(int i=1;i<N;++i) sum[i]=sum[i-1]+mu[i];
}
int calc(int n,int m)
{
int t=std::min(n,m),ans=0;//应该不需要longlong
// for(int k=1;k<=t;++k) ans+=mu[k]*(n/k)*(m/k);//TLE:O(n^2)
for(int las,i=1;i<=t;i=las+1)
{
las=std::min(n/(n/i),m/(m/i));
ans+=(sum[las]-sum[i-1])*(n/i)*(m/i);
}
return ans;
}
int main()
{
Init();
int t=read(),a,b,c,d,k;
while(t--)
a=read()-1,b=read(),c=read()-1,d=read(),k=read(),
a/=k,b/=k,c/=k,d/=k,printf("%d
",calc(b,d)-calc(a,d)-calc(b,c)+calc(a,c));
return 0;
}