luogu2257 YY的gcd【莫比乌斯反演】

题目描述

神犇YY虐完数论后给傻×kAc出了一题

给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对

kAc这种傻×必然不会了，于是向你来请教……

多组输入

输入输出格式

输入格式：

第一行一个整数T 表述数据组数

接下来T行，每行两个正整数，表示N, M

输出格式：

T行，每行一个整数表示第i组数据的结果

输入输出样例

输入样例#1： 

2
10 10
100 100

输出样例#1： 

30
2791

说明

T = 10000

N, M <= 10000000

         蒟蒻自己码的第一道紫题，尽管是在老师讲了之后，A了以后开心地就来写题解。

为了方便推公式，我们用p表示小于m和n的质数，题目要求可以看做是求

那么回忆一下莫比乌斯反演是什么套路

那么上式就可以变形为

(注：打不了向下取整就用C++默认整除了）

然后我们就可以愉快地反演

然后我们可以改变一下顺序，把枚举d前置，然后再枚举i,j，由于d是gcd(i,j)的约数，所以它一定是i,j的约数，把它们的上界都再除以d，得到的就是约数个数，就相当于由枚举约数变成了枚举倍数，于是我们就可以和中间的两个Σ说再见了。

推到现在感觉差不多了？并不是，它显然会T。所以我们还可以进行进一步优化。

我们设T = p*d，则有

然后按照老套路，我们枚举T，由枚举倍数回到枚举约数，就可以得到

我们发现，后面那个Σ是可以预处理的，我们用一个 f[i]数组表示当T == i 时后面那一堆的值，只要在筛素数筛莫比乌斯函数的时候，对于每一个质数p，都在它的倍数T的 f (T)上加上μ（T / p），就可以求出f数组的值。

对于前面那一个，整除分块就可以解决，可以参照这个blog：http://www.cnblogs.com/peng-ym/p/8661118.html。

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define LL long long
using namespace std;
int prime[10000010], f[10000010], v[10000010], sum[10000010], tot, T, n, m;
LL miu[10000010];
void mb(){
    miu[1] = 1;
    for(int i = 2; i <= 10000001; ++ i){
        if(!v[i]){
            prime[++tot] = i;
            miu[i] = -1;
        }
        for(int j = 1; j <= tot && prime[j] * i <= 10000001; ++ j){
            v[i * prime[j]] = 1；
            if(i % prime[j] == 0) break;
            miu[i * prime[j]] = - miu[i];
 	    }
 		
    }//筛素数+求莫比乌斯函数
    for(int i = 1; i <= tot; ++ i)
        for(int j = 1; j * prime[i] <= 10000001; ++ j)
            f[j * prime[i]] += miu[j]; //预处理f数组
    for(int i = 1; i <= 10000001; ++ i)
        sum[i] = sum[i - 1] + f[i];
}
inline int read(){
    int x = 0, f = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){ if(ch == '-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') x = x * 10 + ( ch ^ 48 ), ch = getchar();
    return x * f;
}
inline LL calc(int a, int b){
    LL ans = 0;
    if(a > b) swap(a, b);
    for(int l = 1, r = 0; l <= a; l = r + 1){
        r = min( a / (a / l) , b / (b / l) ); //整除分块
        ans += (LL)( sum[r] - sum[l - 1] ) * (LL)( a / l ) * (LL)( b / l );
    }
    return ans;
}
int main()
{
    T = read();
    mb();
    while(T--){
        n = read(), m = read();
        if(n > m) swap(n, m);
        printf("%lld\n", calc(n, m));
    }
    return 0;
}