題意明確，就是求∑i=ab∑j=cdgcd(i,j)==k
首先容易發現容斥定理，轉化為求
∑i=1n∑j=1mgcd(i,j)==k
我們設f(d)=∑i=1n∑j=1m(gcd(i,j)==d)
F(d)=∑i=1n∑j=1m(d|gcd(i,j))
那麼明顯有F(n)=∑n|df(d)
那麼反演會有f(i)=∑i|dμ(di)F(d)
又顯然F(d)=⌊nd⌋⌊md⌋
時間複雜度一次查詢O(n)
考慮到⌊nd⌋只有n‾‾√種取值，我們可以列舉取值，進行一下除法優化就可以了就優化為一次查詢O(n‾‾√)。

#include<iostream>
#include<cstdio>
 

#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
bool vis[50001];
int mu[50001],prime[50001],n,a,b,c,d,k;
ll sum[50001];
void get_mobious()
{int i,j;
    memset(vis,0,sizeof(vis));
    mu[1]=1;
    int tot=0;
     for (i=2;i<=50000;i++)
      {
          if
 
 (vis[i]==0)
          {
              prime[++tot]=i;
              mu[i]=-1;
          }
          for (j=1;j<=tot;j++)
          {
              if (i*prime[j]>50000) break;
              vis[i*prime[j]]=1;
              if (i%prime[j]==0)
              {
                  mu[i*prime[j]]=0;
                  break
 
;
            }
            else mu[i*prime[j]]=-mu[i];
          }
      }
    sum[0]=0;
    for (i=1;i<=50000;i++)
    sum[i]=sum[i-1]+mu[i];
}
ll cal(int x,int y)
{int r,i,pos;
ll s=0;
    r=min(x,y);
    for (i=1;i<=r;i=pos+1)
    {
        pos=min(x/(x/i),y/(y/i));
        s+=(sum[pos]-sum[i-1])*(x/i)*(y/i);
    }
    return s;
}
int main()
{
    cin>>n;
    get_mobious();
    while (n--)
    {
      scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
         ll ans1=cal(b/k,d/k);
         ll ans2=cal((a-1)/k,d/k);
         ll ans3=cal(b/k,(c-1)/k);
         ll ans4=cal((a-1)/k,(c-1)/k);
         printf("%lld\n",ans1-ans2-ans3+ans4);
    }
}