51nod 1192 莫比乌斯反演

很基础的反演题，练练手。

题目要求：
$$\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p]$$
枚举p：
$$\sum _{k=1}^{primes}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=prime[k]]$$
然后把后面单独拉出来
$$\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=d]$$
化简：（套路）
$$\sum _{i=1}^{n/d}\mu [i]\ast \frac{n}{i\ast d}\ast \frac{m}{i\ast d}$$
带回去：
$$\sum _{k=1}^{primes}\sum _{i=1}^{n/d}\mu [i]\ast \frac{n}{i\ast prime[k]}\ast \frac{m}{i\ast prime[k]}$$
一看看出来id套路T=id,化简得：
$$\sum _{T=1}^n\ast \frac{n}{T}\ast \frac{m}{T}\ast \sum _{d|T}\mu [T/d]$$
令
$$g(T)=\sum _{d|T}\mu [T/d]$$
原式为：
$$\sum _{T=1}^n\ast \frac{n}{T}\ast \frac{m}{T}\ast g(T)$$
对于这个式子，很显然n/Tm/T可以整除分块。这样多次查询+预处理的复杂度是
$$O(T\ast \sqrt{n}+n\ast ln(n))$$
代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<queue>
#include<stack>
#include<map>
#include<ctime>
#define up(i,x,y) for(int i = x;i <= y;i ++)
#define down(i,x,y) for(int i = x;i >= y;i --)
#define mem(a,b) memset((a),(b),sizeof(a))
#define mod(x) ((x)%MOD)
#define lson p<<1
#define rson p<<1|1
using namespace std;
typedef long long ll;
const int SIZE = 5000010;
const int INF = 2147483640;
const double eps = 1e-8;

inline void RD(int &x)
{
    x = 0;  char c; c = getchar();
    bool flag = 0;
    if(c == '-')    flag = 1;
    while(c < '0' || c > '9')   {if(c == '-')   {flag = 1;} c = getchar();}
    while(c >= '0' && c <= '9') x = (x << 1) + (x << 3) + c - '0',c = getchar();
}

int primes,prime[SIZE],mu[SIZE];
bool vis[SIZE];
int g[SIZE];

void Init()
{  
    mu[1]=1;
    primes=0;
    for (ll i=2;i<=5000000;i++)
    {  
        if (!vis[i])//素数
        {  
            prime[primes++]=i;
            mu[i]=-1; 
        }  
        for (ll j=0;j<primes && i*prime[j]<=5000000;j++)
        {  
            vis[i*prime[j]]=1;  
            if(i%prime[j])//i pj互素
            {
                mu[i*prime[j]]=-mu[i];
            }
            else
            {
                mu[i*prime[j]]=0;
                break;
            }  
        }  
    }
    for(int i = 0;i < primes;i ++)
    {
    //	printf("prime:%lld\n",prime[i]);
    	for(int j = prime[i];j <= 5000000;j += prime[i])
    	{
    		g[j] += mu[j/prime[i]];
    	}
    }
    for(int i = 1;i <= 5000000;i ++)	g[i] += g[i-1];
}

void solve(int a,int b)
{
	if(a > b)	swap(a,b);
	int last;
	ll ans = 0;
	for(int i = 1;i <= a;i = last+1)
	{
		last = min(a/(a/i),b/(b/i));
		ans += (ll)(a/i)*(ll)(b/i)*(ll)(g[last]-g[i-1]);
	}
	printf("%lld\n",ans);
}

int main(int argc, char const *argv[])
{
	Init();
	int t;
	scanf("%d",&t);
	while(t--)
	{
		int a,b;
		scanf("%d%d",&a,&b);
		solve(a,b);
	}
	return 0;
}