[luogu] P2257 YY的GCD(莫比乌斯反演)

最新推荐文章于 2022-04-28 20:05:12 发布

原创最新推荐文章于 2022-04-28 20:05:12 发布 · 290 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#GCD #莫比乌斯反演 #欧拉函数

数论-各种数论函数同时被 2 个专栏收录

9 篇文章

订阅专栏

math

8 篇文章

订阅专栏

本文详细解析了一道关于求解特定质数对数量的问题，运用莫比乌斯反演和欧拉函数两种方法，介绍了算法实现细节及优化技巧，适合对数学算法和数论感兴趣的读者。

题目

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

数据范围: $T(T\leq 1000)$ 组询问, $N,M \leq 10^7$

题解

这是一道练习莫比乌斯反演的模板题,下面先给出莫比乌斯反演的做法,最后再一种欧拉函数计算的思路

对于这种 $gcd(x,y)$ 二元组的题,我们设 $F[n],f[n]$ 分别为 $gcd(x,y)$ 为 $n$ 倍数与正好为 $n$ 的个数(范围内)

即: $F[n]=\sum_{n|gcd(x,y)} 1 = [\frac{N}{n}][\frac{M}{n}],f[n]=\sum_{gcd(x,y)==d}1$ ,两者关系有 $F[n] = \sum_{n|d} f[d]$ ,对其进行莫比乌斯反演

$f[n]=\sum_{n|d} \mu ({\frac{d}{n}})F[d]$ ,则原题的 $ans=\sum_{p \in prime}f[p]=\sum_{p \in prime} \sum_{p|d}\mu(\frac{d}{p})F[d]$

继续化解得 $ans=\sum_{p \in prime} \sum_{i=1}^{min(\frac{N}{p},\frac{M}{p})} \mu [i] F [ip]=\sum_{p \in prime} \sum_{i=1}^{min(\frac{N}{p},\frac{M}{p})} \mu [i] [\frac{N}{ip}][\frac{M}{ip}]$ ,这里有个 $trick$ 如果是单次询问则可以直接处理 $\mu []$ , $p[]$ (最外层)前缀和再利用整除分块(对内层求和上标分类)进行计算,时间复杂度大概单次 $O(min(N,M))$ ,不过由于询问过多还需要进一步化解(如有需求可以在单次询问情况写一写)

if(N > M) swap(N, M);
    int l, r;
    for(l = 2; l <= N; l = r + 1){
      r = min(N / (N / l), M / (M / l));
      ans += solve(N / l, M / l) * (sum[r] - sum[l - 1]);
      //solve()对最内层求和,sum[]维护最外层前缀和
    }

再思考一下,莫比乌斯反演求和的关键是整除分块,那么我们考虑将最内层的 $[\frac{N}{ip}][\frac{M}{ip}]$ 移至外层进行进一步化解

$\\ \sum_{p \in prime} \sum_{i=1}^{min(\frac{N}{p},\frac{M}{p})} \mu [i] [\frac{N}{ip}][\frac{M}{ip}]\\ =\sum_{p \in prime}\sum_{p|j}^{min(N,M)} [\frac{N}{j}][\frac{M}{j}] \mu [\frac{j}{p}]\\ =\sum_{j=1}^{min(N,M)} [\frac{N}{j}][\frac{M}{j}] \sum_{p \in prime, p|j}\mu [\frac{j}{p}]$ ,通过预处理 $\mu [\frac{j}{p}]$ 前缀和(大概 $O(cnt\_prime*log(MAXvalue))$ )后针对每次询问即可以在 $O(\sqrt {min(N,M)})$ 内处理得到结果

PS:这里注意不要乱用 $long\ long$ (莫名 $TLE$ )

一开始我写的

LL N,M;
LL solve(LL, LL){
    LL l,r;
    ...
}

调试了近半小时才,改成 $int$ 才发现竟然快了不止一倍（╯#-_-)╯╧═╧

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;

using LL = long long;
const int MAXN = 1e7 + 5;
int T, N, M;
bool noprime[MAXN];
int prime[MAXN], cnt_p, mu[MAXN];
void Euler_sieve(int top);
int g[MAXN], pre_g[MAXN];
LL solve(int, int);

int main(){
  Euler_sieve(1e7);
  int i, j;
  for(i = 1; i <= cnt_p; i++)
    for(j = prime[i]; j <= 1e7; j += prime[i])
      g[j] += mu[j / prime[i]];
  for(i = 2; i <= 1e7; i++) pre_g[i] = pre_g[i - 1] + g[i];
  cin >> T;
  while(T--){
    cin >> N >> M;
    LL ans = 0;
    ans = solve(N, M);
    cout << ans << endl;
  }
  return 0;
}

LL solve(int N, int M){
  if(N > M) swap(N, M);
  int l, r; LL res = 0;
  for(l = 2; l <= N; l = r + 1){
    r = min(N / (N / l), M / (M / l));
    res += (LL)(N / l) * (M / l) * (pre_g[r] - pre_g[l - 1]);
  }
  return res;
}

void Euler_sieve(int top){
  int i, j;
  mu[1] = 1;
  for(i = 2; i <= top; i++){
    if(!noprime[i])
      prime[++cnt_p] = i,  mu[i] = -1;
    for(j = 1; j <= cnt_p && prime[j] * i <= top; j++){
      noprime[prime[j] * i] = true;
      if(i % prime[j] == 0)
        break;
      mu[i * prime[j]] = -mu[i];
    }
  }
}

最后这里再给出欧拉函数的思路(这里限制 $N=M$ GCD)

主要思路可以看gcd求和中第二第三部分

这里介绍第三部分思路(主要原文另一种写的比较清楚)

$g[p]=\sum_{i=2}^{[\frac{N}{p}]} \varphi (i)$ , $g[p]=\sum_{gcd(x,y)=p,x < y}1$ ,则每个素数 $p$ 贡献为 $g[p](x<y)+1(x=y=p)+g[p](x>y)=2g[p]+1$

即可以计算得出GCD

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;

using LL = long long;
const int MAXN = 1e7 + 5;
LL N, ans;
bool noprime[MAXN];
int prime[MAXN], cnt_p, phi[MAXN];
void Euler_sieve(int top);
LL f[MAXN];

int main(){
  cin >> N;
  Euler_sieve(N);
  int i, j;
  for(i = 2; i <= N; i++)
    if(!noprime[i])
      for(j = 2; j * i <= N; j++)
        f[i] += phi[j];
  for(i = 2; i <= N; i++)
    if(!noprime[i]){
      ans += 2 * f[i];
      ans++;
    }
  cout << ans << endl;
  return 0;
}

void Euler_sieve(int top){
  int i, j;
  phi[1] = 1;
  for(i = 2; i <= top; i++){
    if(!noprime[i])
      prime[++cnt_p] = i,  phi[i] = i - 1;
    for(j = 1; j <= cnt_p && prime[j] * i <= top; j++){
      noprime[prime[j] * i] = true;
      if(i % prime[j] == 0){
        phi[i * prime[j]] = phi[i] * prime[j];
        break;
      }else{
        phi[i * prime[j]] = phi[i] * (prime[j] - 1);
      }
    }
  }
}