首先可以想到分为四个前缀区间进行加加减减,考虑
[1,a],[1,b]
这组:
∑i=1a∑j=1b[gcd(i,j)==k]
=∑i=1a/k∑j=1b/k[gcd(i,j)==1]
=∑i=1a/k∑j=1b/k∑t|i,t|jμ(t)
令 n=a/k , m=b/k ,则原式可改写为
∑t=1min(n,m)μ(t)[nt][mt]
这部分可以在 O(max{b,d}−−−−−−−−√) 分段求
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<set>
#include<map>
#include<queue>
#include<stack>
#include<vector>
#include<bitset>
using namespace std;
typedef long long LL;
int miu[51000];
int prim[51000],primm;
bool valid[51000];
void mobius(int N)
{
miu[1]=1;
for (int i=2;i<=N;i++)
{
if (!valid[i]) prim[++primm]=i,miu[i]=-1;
for (int j=1;j<=primm&&i*prim[j]<=N;j++)
{
valid[i*prim[j]]=1;
if (i%prim[j]==0)
{
miu[i*prim[j]]=0;
break;
}
else miu[i*prim[j]]=-miu[i];
}
}
}
int a,b,c,d,n,k;
LL sum[51000];
LL F(int n,int m)
{
n/=k,m/=k;
LL res=0;
for (int i=1;i<=min(n,m);i++)
{
int pos=i;
i=min(n/(n/i),m/(m/i));
LL a=n/i,b=m/i;
res+=(sum[i]-sum[pos-1])*a*b;
}
return res;
}
int main()
{
mobius(50000);
for (int i=1;i<=50000;i++)
sum[i]=sum[i-1]+(LL)miu[i];
scanf("%d",&n);
while(n--)
{
scanf("%d %d %d %d %d",&a,&b,&c,&d,&k);
LL ans=F(b,d)-F(a-1,d)-F(b,c-1)+F(a-1,c-1);
printf("%lld\n",ans);
}
return 0;
}