[SDOI2015]约数个数和

传送门

这个题……主要问题在于\(d(ij)\)应该怎么变形……容易想到改变成gcd的形式，不过不知道怎么改……

后来听大佬说有这么一个性质：

\[d(ij) = \sum_{x|i}\sum_{y|j}[gcd(x,y) = 1] \]

这个不知道怎么严格证明……不过可以感性理解一下，就是首先肯定是要枚举i，j的因子的，然后如果这两个枚举的数不互质的话，那就可以把他们的gcd挪到一边使其变为互质的，也就是说如果二者不互质还计算，就会导致重复计数。（好吧这确实只是感性理解……）

那我们要求的就是

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

这个式子的形式好像特别熟悉啊……把\(\sum\)之间的限制互换，再用莫比乌斯函数的性质，再替换一下，这些熟悉的操作之后，这个式子就成了这样：

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

然后发现外层的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 = 50005;
const int INF = 1000000009;
const double eps = 1e-7;

int read()
{
    int 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,m,T,p[M],mu[M];
ll pre[M],sum[M],tot;
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 <= tot;j++)
      {
	 np[i * p[j]] = 1;
	 if(!(i % p[j])) break;
	 mu[i * p[j]] = -mu[i];
      }
   }
   rep(i,1,M-2) pre[i] = pre[i-1] + mu[i];
   rep(k,1,M-2)
   {
      ll cur = 0;
      for(int i = 1,j;i <= k;i = j + 1)
      {
	 j = k / (k / i);
	 cur += (ll)(j - i + 1) * (ll)(k / i);
      }
      sum[k] = cur;
   }
}

ll solve(ll a,ll b)
{
   ll lim = min(a,b),cur = 0;
   for(int i = 1,j;i <= lim;i = j + 1)
   {
      j = min(a / (a / i),b / (b / i));
      cur += (pre[j] - pre[i-1]) * sum[a / i] * sum[b / i];
   }
   return cur;
}

int main()
{
   euler();
   T = read();
   while(T--) n = read(),m = read(),printf("%lld\n",solve(n,m));
   return 0;
}