[POI2007]ZAP-Queries

传送门

题目要求：求出

\[\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j) = d] \]

我们先假设m<n.
这个还是老套路先把d除进去然后互相替代……

\[\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}[gcd(i,j) = 1] \]

之后我们用莫比乌斯函数的性质！再把互相之间的限制条件调换一下！就直接得到结果就是：

\[\sum_{p=1}^{\frac{n}{d}}\mu(p)\left\lfloor\frac{n}{dp}\right\rfloor\left\lfloor\frac{m}{dp}\right\rfloor \]

然后就做完啦，先把给定的n，m除以d，剩下的直接整除分块即可。

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<set>
#include<vector>
#include<map>
#include<queue>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('\n')
#define fr friend inline
#define y1 poj
#define mp make_pair
#define pr pair<int,int>
#define fi first
#define sc second
#define pb push_back
#define I puts("bug")

using namespace std;
typedef long long ll;
const int M = 100005;
const int INF = 1000000009;
const double eps = 1e-7;
const double pi = acos(-1);
const ll mod = 1e9+7;

ll read()
{
    ll ans = 0,op = 1;char ch = getchar();
    while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
    while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
    return ans * op;
}

int n,a,b,d,tot,p[M],sum[M],mu[M];
bool np[M];

void euler()
{
   np[1] = 1,mu[1] = 1;
   rep(i,2,M-2)
   {
      if(!np[i]) p[++tot] = i,mu[i] = -1;
      for(int j = 1;i * p[j] <= M-2;j++)
      {
     np[i*p[j]] = 1;
     if(!(i % p[j])) {mu[i * p[j]] = 0;break;}
     mu[i*p[j]] = -mu[i];
      }
   }
   //rep(i,1,10) printf("%d ",mu[i]);enter;
   rep(i,1,M-2) sum[i] = sum[i-1] + mu[i];
}

int solve(int k1,int k2)
{
   //printf("%d %d\n",k1,k2);
   int m = min(k1,k2),ans = 0;
   for(int i = 1,j;i <= m;i = j + 1)
   {
      j = min(k1 / (k1 / i),k2 / (k2 / i));
      //printf("%d\n",j);
      ans += (k1 / i) * (k2 / i) * (sum[j] - sum[i-1]);
   }
   return ans;
}

int main()
{
   euler();
   n = read();
   while(n--)
   {
      a = read(),b = read(),d = read();
      int k1 = a / d,k2 = b / d;
      printf("%d\n",solve(k1,k2));
   }
   return 0;
}