BZOJ2820 YY的GCD 【莫比乌斯反演】

2820: YY的GCD

Time Limit: 10 Sec Memory Limit: 512 MB
Submit: 2398 Solved: 1269
[Submit][Status][Discuss]
Description

神犇YY虐完数论后给傻×kAc出了一题给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对kAc这种
傻×必然不会了，于是向你来请教……多组输入
Input

第一行一个整数T 表述数据组数接下来T行，每行两个正整数，表示N, M
Output

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

2

10 10

100 100

Sample Output

30

2791

HINT

T = 10000

N, M <= 10000000

题解

还是经典的莫比乌斯反演
我们设f(n)表示gcd等于n的对数
我们设F(n)表示n|gcd的对数
则有
F(n)=⌊Nn⌋⌊Mn⌋

f(n)=∑n|dμ(dn)F(d)

=∑n|dμ(dn)⌊Nn⌋⌊Mn⌋

=∑Ni=1μ(i)⌊Nn⌋⌊Mn⌋

ans=∑Np∈prime∑Ni=1μ(i)⌊Nn⌋⌊Mn⌋

=∑NT=1⌊Nn⌋⌊Mn⌋∗∑Np|Tμ(Tp)

对于前半部分∑NT=1⌊Nn⌋⌊Mn⌋我们可以分块O(2N−−√)

后半部分我们可以预处理，枚举每个质数以及它的倍数
由于质数的个数大约是nlogn，而每个质数倍数的枚举均摊是logn的，加上线性筛，总的预处理复杂度O(n)

总复杂度O(T∗N−−√+N)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<bitset>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u]; k != -1; k = ed[k].nxt)
using namespace std;
const int maxn = 10000005,maxm = 100005,INF = 1000000000;
inline int RD(){
    int out = 0,flag = 1; char c = getchar();
    while (c < 48 || c > 57) {if (c == '-') flag = -1; c = getchar();}
    while (c >= 48 && c <= 57) {out = (out << 1) + (out << 3) + c - '0'; c = getchar();}
    return out * flag;
}
bool isn[maxn];
int prime[maxn],primei = 0,mu[maxn];
LL A[maxn];
void init(){
    mu[1] = 1;
    for (int i = 2; i <= 10000000; i++){
        if (!isn[i]) prime[++primei] = i,mu[i] = -1;
        for (int j = 1; j <= primei && i * prime[j] <= 10000000; j++){
            isn[i * prime[j]] = true;
            if (i % prime[j] == 0) {mu[i * prime[j]] = 0; break;}
            mu[i * prime[j]] = -mu[i];
        }
    }
    REP(i,primei) for (int j = 1; prime[i] * j <= 10000000; j++) A[prime[i] * j] += mu[j];
    REP(i,10000000) A[i] += A[i - 1];
}
int main(){
    init();
    int Tim = RD(),N,M;
    LL ans;
    while (Tim--){
        N = RD(); M = RD(); ans = 0; if (N > M) swap(N,M);
        for (int i = 1,nxt; i <= N; i = nxt + 1){
            nxt = min(N / (N / i),M / (M / i));
            ans += (A[nxt] - A[i - 1]) * (N / i) * (M / i);
        }
        printf("%lld\n",ans);
    }
    return 0;
}