BZOJ2820 YY的GCD 莫比乌斯+系数前缀和

/**
题目：BZOJ2820 YY的GCD
链接：http://www.cogs.pro/cogs/problem/problem.php?pid=2165
题意：神犇YY虐完数论后给傻×kAc出了一题
给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对
kAc这种傻×必然不会了，于是向你来请教……
T = 10000
N, M <= 10000000
思路：
f(n)表示gcd==n的对数。
g(n)表示gcd的n的倍数的对数。
u(d/p) = mu[d/p];

ans = sigma[p是质数,p<=min(n,m)]sigma[p|d] (u(d/p)*g(d))

    = sigma[p是质数,p<=min(n,m)]sigma[p|d] (u(d/p)*(n/d)*(m/d))

    = sigma[1<=d<=min(n,m)](n/d)*(m/d)sigma[p是d的约数且p是素数](u(d/p));

参考:http://www.cnblogs.com/candy99/p/6209609.html

*/
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <set>
#include <iostream>
#include <vector>
#include <map>
using namespace std;
typedef long long LL;
#define ms(x,y) memset(x,y,sizeof x)
typedef pair<int, int> P;
const LL INF = 1e10;
const int mod = 1e9 + 7;
const int maxn = 1e7 + 100;
int prime[maxn], tot, not_prime[maxn];
int mu[maxn], sum[maxn];
inline int read()
{
    char  c = getchar();
    int x = 0;
    while(c<'0'||c>'9'){
        c = getchar();
    }
    while(c>='0'&&c<='9') x = x*10+c-48,c = getchar();
    return x;
}
//法2.
//线性筛
//g[i*p[j]]
//当p[j]|i时结果显然为miu(i)
//否则考虑mu(i*p[j]/pp),当p[j]=pp时为mu[i]，p[j]!=pp时的所有的和就是-g(i)，所以总的结果为mu(i)-g(i)
/*
void mobius()
{
    mu[1] = 1;
    tot = 0;
    for(int i = 2; i < maxn; i++){
        if(!not_prime[i]){
            mu[i] = -1;
            prime[++tot] = i;
            sum[i] = 1;
        }
        for(int j = 1; prime[j]*i<maxn; j++){
            not_prime[prime[j]*i] = 1;
            if(i%prime[j]==0){
                mu[prime[j]*i] = 0;
                sum[prime[j]*i] = mu[i];
                break;
            }
            mu[prime[j]*i] = -mu[i];
            sum[prime[j]*i] = mu[i]-sum[i];
        }
    }

    for(int i = 1; i < maxn; i++) sum[i] += sum[i-1];
}*/
//法1.
//只需要枚举每个素数，将他的倍数的g更新就可以了
//由于有1/1+1/2+1/3+...+1/n=O(logn)这个结论
//因此每个质数枚举时是均摊O(logn)的(*n后好想，是nlogn，但是质数只有n/logn个)
//而质数恰好有O(n/logn)个 因此暴力枚举就是O(n)的
void mobius()
{
    mu[1] = 1;
    tot = 0;
    for(int i = 2; i < maxn; i++){
        if(!not_prime[i]){
            mu[i] = -1;
            prime[++tot] = i;
        }
        for(int j = 1; prime[j]*i<maxn; j++){
            not_prime[prime[j]*i] = 1;
            if(i%prime[j]==0){
                mu[prime[j]*i] = 0;
                break;
            }
            mu[prime[j]*i] = -mu[i];
        }
    }

    for(int i = 1; i <= tot; i++){
        for(int j = prime[i]; j < maxn; j+=prime[i]){
            sum[j] += mu[j/prime[i]];
        }
    }

    for(int i = 1; i < maxn; i++) sum[i] += sum[i-1];

}
LL solve(int n,int m)
{
    if(n>m) swap(n,m);
    LL ans = 0;
    int last;
    for(int i = 1; i <= n; i = last+1){
        last = min(n/(n/i),m/(m/i));
        ans += (LL)(sum[last]-sum[i-1])*(n/i)*(m/i);
    }
    return ans;
}


int main()
{
    freopen("YYnoGCD.in","r",stdin);
    freopen("YYnoGCD.out","w",stdout);
    int n, m;
    int T;
    T = read();
    mobius();
    while(T--)
    {
        n = read();
        m = read();
        printf("%lld
",solve(n,m));
    }
    return 0;
}