对于每次询问 $a,b,c,d,n$ 答案就是 $f[n]=sum_{i=a}^{b}sum_{j=c}^{d}[gcd(i,j)==n]$
看到熟悉的 $gcd$ ,莫比乌斯反演
设 $F[n]=sum_{i=a}^{b}sum_{j=c}^{d}[n|gcd(i,j)]$,那么 $F[n]=sum_{n|d}f[d]$
并且易知 $F[n]=(left lfloor frac{b}{n} ight floor-left lfloor frac{a-1}{n} ight floor) (left lfloor frac{d}{n} ight floor-left lfloor frac{c-1}{n} ight floor)$
直接反演:$f[n]=sum_{n|d} mu (d/n) F[d]=sum_{n|d} mu (d/n) (left lfloor frac{b}{n} ight floor-left lfloor frac{a-1}{n} ight floor) (left lfloor frac{d}{n} ight floor-left lfloor frac{c-1}{n} ight floor)$
枚举 $n$ 的倍数 $k$:变成 $f[n]=sum_{k} mu (k) (left lfloor frac{b}{kn} ight floor-left lfloor frac{a-1}{kn} ight floor) (left lfloor frac{d}{kn} ight floor-left lfloor frac{c-1}{kn} ight floor)$
因为 $left lfloor frac{a-1}{kx} ight floor=left lfloor frac{left lfloor frac{a-1}{x} ight floor}{k} ight floor$
所以 $f[n]=sum_{k} mu (k) (left lfloor frac{ left lfloor frac{b}{n} ight floor }{k} ight floor-left lfloor frac{left lfloor frac{a}{n} ight floor}{k} ight floor) (left lfloor frac{left lfloor frac{d}{n} ight floor}{k} ight floor-left lfloor frac{left lfloor frac{c}{n} ight floor}{k} ight floor)$
然后就可以数论分块了
总复杂度 $O( n sqrt(n) )$
注意代码中的 $k$ 就是指前面的 $n$
#include<iostream> #include<cstdio> #include<algorithm> #include<cmath> #include<cstring> using namespace std; typedef long long ll; inline int read() { int x=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9') { if(ch=='-') f=-1; ch=getchar(); } while(ch>='0'&&ch<='9') { x=(x<<1)+(x<<3)+(ch^48); ch=getchar(); } return x*f; } const int N=2e6+7,M=5e4; int Q,a,b,c,d,k,pri[N],mu[N],tot; ll ans; bool not_pri[N]; void pre() { mu[1]=1; for(int i=2;i<=M;i++) { if(!not_pri[i]) pri[++tot]=i,mu[i]=-1; for(int j=1;j<=tot;j++) { ll t=1ll*i*pri[j]; if(t>M) break; not_pri[t]=1; if(!(i%pri[j])) break; mu[t]=-mu[i]; } } for(int i=2;i<=M;i++) mu[i]+=mu[i-1]; } int main() { pre(); Q=read(); while(Q--) { a=read(),b=read(),c=read(),d=read(),k=read(); ans=0; a--,c--;/*注意a--,c--之后a,c可能为0*/ a/=k,b/=k,c/=k,d/=k; int mi=min(b,d); for(int l=1,r;l<=mi;l=r+1) { r=min( b/(b/l) , d/(d/l) ); if(a/l) r=min(r,a/(a/l)); if(c/l) r=min(r,c/(c/l));//注意特判a,c等于0的情况 ans+=1ll*(mu[r]-mu[l-1])*(b/l-a/l)*(d/l-c/l); } printf("%lld ",ans); } return 0; }