LuoguP2257 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

题解

以下均设\(n<=m\)

\[\large{ \begin{align*} &\sum_{i=1}^{n}\sum_{j=1}^{m}{[gcd(i,j)=p](p\in prime)}\\ &=\sum_{p\in prime}\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{p} \rfloor}{[gcd(i,j)=1]}\\ &=\sum_{p\in prime}\sum_{i=1}^{\lfloor \frac{n}{p} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{p} \rfloor}\sum_{d|gcd(i,j)}{\mu(d)}\\ &=\sum_{p\in prime}\sum_{d=1}^{\lfloor \frac{n}{p}\rfloor}{\mu(d)\lfloor \frac{n}{dp} \rfloor \lfloor \frac{m}{dp} \rfloor} \end{align*} } \]

这样复杂度是\(O(p\sqrt{n})\)的（p为素数个数），会超时，要继续推
然而我只会推到这里了，数学题真毒瘤.jpg
题解里的大爷都是神仙.jpg
新操作get：在式子化到最简的时候，我们可以考虑一下更换枚举项，把这个式子搞成可以预处理的，然后降低复杂度

\[\large { 设T=dp\\ \begin{align*} &\sum_{p\in prime}\sum_{d=1}^{\lfloor \frac{n}{p}\rfloor}{\mu(d)\lfloor \frac{n}{dp} \rfloor \lfloor \frac{m}{dp} \rfloor}\\ &=\sum_{p\in prime}\sum_{d=1}^{\lfloor \frac{n}{p}\rfloor}{\mu(d)\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}\\ &=\sum_{T=1}^{n}{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}\sum_{p|T,p\in prime}{\mu(\frac{T}{p})} \end{align*} } \]

于是预处理后面的那块就好，具体做法是对每个素数p，枚举它的倍数T，加上\(\mu(\frac{T}{p})\)

式子推出来代码就很容易码了qwq

#include <bits/stdc++.h>
using namespace std;

#define ll long long
#define N 10000020
int n, m, cnt = 0;
int mu[N], vis[N], p[N];
ll f[N], sum[N];

void init() {
    mu[1] = 1;
    for(int i = 2; i < N; ++i) {
        if(!vis[i]) {p[++cnt] = i; mu[i] = -1;}
        for(int j = 1; j <= cnt && p[j] * i < N; ++j) {
            vis[p[j] * i] = 1;
            if(i % p[j] == 0) break;
            mu[i * p[j]] -= mu[i];
        }
    }
    for(int i = 1; i <= cnt; ++i) 
    	for(int j = 1; j * p[i] < N; j ++)
    		f[j * p[i]] += mu[j];
    for(int i = 1; i < N; ++i) sum[i] = sum[i - 1] + f[i];
}


ll calc(int a, int b) {
    ll s = 0;
    for(int l = 1, r; l <= a; l = r + 1) {
    	r = min(a / (a / l), b / (b / l));
		s += 1ll * (sum[r] - sum[l - 1]) * (a / l) * (b / l);
	}
    return s;
}

int main() {
#ifndef ONLINE_JUDGE
freopen("1.in","r",stdin);
freopen("1.out","w",stdout);
#endif
    init();
    int T;
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d", &n, &m);
        if(n > m) swap(n, m);
        printf("%lld\n", calc(n, m));
    }
    return 0;
}