Primes in GCD Table SPOJ - PGCD(莫比乌斯反演+分段)

Primes in GCD Table SPOJ - PGCD

Johnny has created a table which encodes the results of some operation -- a function of two arguments. But instead of a boring multiplication table of the sort you learn by heart at prep-school, he has created a GCD (greatest common divisor) table! So he now has a table (of height a and width b), indexed from (1,1) to (a,b), and with the value of field (i,j) equal to gcd(i,j). He wants to know how many times he has used prime numbers when writing the table.
Input

First, t ≤ 10, the number of test cases. Each test case consists of two integers, 1 ≤ a,b < 107.
Output

For each test case write one number - the number of prime numbers Johnny wrote in that test case.
Example

Input:
2
10 10
100 100

Output:
30
2791

题意：

给出一个数$N,M，求1 \le x \le N,1 \le y \le M，且gcd(x, y)为素数的数对x, y$的数量。

分析:

假设$p$为素数，那么这道题就转化成了求

$$gcd(x,y) = p$$ 的数对有多少个。

这样就和我们之前做过的题求$gcd(x,y) = k$的数对有多少个是一样的了

---

回想以下之前的求解步骤

1）要求$gcd(x,y) = k$的数对多少个相当于求$gcd(\frac{x}{k},\frac{y}{k}) = 1$的数对有多少个

2）利用莫比乌斯反演：
设$f(d)表示gcd(x,y) = d的个数$

设$F(d)表示gcd(x,y)=d以及d的倍数的个数$

那么根据莫比乌斯反演公式

$$f(n) = \sum_{n|d}\mu(\frac{d}{n})F(d)$$

而且我们易求$F(d) = \frac{n}{d} \cdot \frac{m}{d}$

所以答案实际上是求

$$f(1) = \sum_{1|d}\mu(d)F(d)$$

此时的$F(d) = \frac{\frac{n}{k}}{d}\cdot \frac{\frac{m}{k}}{d}$

---

对于这道题目，求$gcd(x,y) = p的个数$是不是就是和上面解法一样了，求$gcd(\frac{x}{p},\frac{y}{p}) = 1$的个数，只不过p是不知道的，需要我们枚举素数p

因此我们得到公式

$$ans = \sum_{p}^{min(n,m)}\sum_{1|d}^{min(n,m)}\mu(\frac{d}{1})F(d)$$

$$=\sum_{p}^{min(n,m)}\sum_{1|d}^{min(n,m)}\mu(d)\frac{\frac{n}{p}}{d}\cdot\frac{\frac{m}{p}}{d}$$

$$\sum_{p}^{min(n,m)}\sum_{1|d}^{min(n,m)}\mu(d)\frac{n}{pd}\cdot \frac{m}{pd}$$

但是直接这样枚举的话会超时，那么我们需要考虑优化

令$t = pd，即d = \frac{t}{p}$

通过变量代换得到

$$ans = \sum_{t=1}^{min(m,n)}\sum_{p|t}\mu(\frac{t}{p})\frac{n}{t}\cdot \frac{m}{t}$$

$$= \sum_{t=1}^{min(m,n)}\frac{n}{t}\cdot \frac{m}{t}\sum_{p|t}\mu(\frac{t}{p})$$

设$a[t] = \sum_{p|t}\mu(\frac{t}{p})$

$$= \sum_{t=1}^{min(m,n)}\frac{n}{t}\cdot \frac{m}{t}a[t]$$
因此我们可以预处理出来 $a[t]$

同时我们还发现$\frac{n}{t}\cdot \frac{m}{t}$因为是整数除法向下取整，因此中间必定有连续的一段（比如$t-1,t,t+1,t+2$）这是$\frac{n}{t}\cdot \frac{m}{t}$的值是一样的只不过$a[t]$不同，所以我们可以进行分段，预处理出$a[t]$的前缀和，这样只要我们知道了这一段的范围，就能通过前缀和O(1)查询这一段的$a[t]$的和了然后直接乘上$\frac{n}{t}\cdot \frac{m}{t}$的值，最后求和

那么这个每一段的范围怎么求呢？

我们可以求出$d_1 = \frac{n}{t},d_2 = \frac{m}{t}$,因此每个这一段的最大范围是$r_1 = \frac{n}{d_1},r_2 = \frac{m}{d_2}$,这样$min(r_1,r_2)$便是这一段的范围

证明可见下面这道题 Ice Rain（数论+分段）是分段思想的简单应用

code：

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
typedef long long ll;
const int maxn = 1e7 + 5;
const int N = 1e6 + 10;
bool vis[maxn];
int prime[N],pnum,mu[maxn],a[maxn];
void init(int n){
    mu[1] = vis[1] = 1;
    pnum = 0;
    //预处理处莫比乌斯函数
    for(int i = 2; i <= n; i++){
        if(!vis[i]){
            prime[pnum++] = i;
            mu[i] = -1;
        }
        for(int j = 0; j < pnum; j++){
            if(i * prime[j] > n) break;
            vis[i*prime[j]] = 1;
            if(i % prime[j] == 0){
                mu[i*prime[j]] = 0;
                break;
            }
            mu[i*prime[j]] = -mu[i];
        }
    }
    //求a[t]与前缀和
    for(int i = 1; i <= n; i++){
        if(mu[i]){
            for(int j = 0; j < pnum; j++){
                if(i * prime[j] > n) break;
                a[i*prime[j]] += mu[i];
            }
        }
    }
    for(int i = 1; i <= n; i++) a[i] += a[i-1];
}
int main(){
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    init(10000000);
    int t,n,m;
    scanf("%d",&t);
    while(t--){
        scanf("%d%d",&n,&m);
        if(n > m) swap(n,m);
        ll ans = 0;
        for(int i = 1; i <= n; i++){
            int tmp1 = n / i,tmp2 = m / i;
            int next = min(n/tmp1,m/tmp2);
            ans += (ll)tmp1 * tmp2 * (a[next] - a[i-1]);
            i = next;
        }
        printf("%lld\n",ans);
    }   
    return 0;
}