SPOJ-PGCD Primes in GCD Table (Mobius反演 好题)

题意：

给你a，b，求gcd(i, j) = k (k为素数)的个数，1 <= i <= a, 1 <= j <= b。

解题思路：

Mobius反演有两种形式：

设 F(n) = sigma( f(d) ), 其中d | n (也就是n%d==0) 。则 f(n) = sigma( F(d)*mu(n/d) )

设 F(d) = sigma( f(n) )，其中d | n 。则f(d) = sigma( F(n)*mu(n/d) )。

                     1          (x = 1)

其中mu(x) =   (-1)^k   (x含有k个不同素因子且每个素因子有且仅有一个)

                        0          (其他情况)

其实Mobius反演就是偏序集上的容斥原理，mu函数就是容斥系数。

要求gcd(i, j) = k，设F(x)为gcd(i, j)为x的倍数的个数，f(x)为gcd(i,j)为x的个数。则F(x) = (a/x)*(b/x)

根据公式可以得到 f(k) = sigma( F(n)*mu(n/k))，其中k | n，n <= a && n <= b。

假如我只需要求一个gcd(i, j) = k的个数，可以对于a,b都除k，则转化成求gcd(i, j) = 1 (1 <= i <= a , 1<= j <= b)。

则答案即为f(1) = sigma( F(x)*mu(x) ) ( x <= min(a, b) )，然后可以求出mu函数，1~min(a,b)枚举得到答案

这样的复杂度为O(n)的，有点高，F(x) = (a/x)*(b/x)，很容易可以看出F(x)是递减的且为若干段值相等的段，即有些连续部分值都相等，所以可以分段处理，这样子的复杂度是O( sqrt(n) )的。

好，现在求单独gcd(i ,j) = k的已经可以处理了，不过题目要求k是所有的素数，很显然不能枚举所有的素数。

重新观察下上面的式子 f(k) = sigma( F(n)*mu(n/k) )，对于所有的f(k)都加起来为

sigma( F(n) * ( mu(n/k1) + mu(n/k2) +mu(n/k3) +... ) )其中n >= 2，ki为不大于n的素因子。

设cnt(n) = F(n) * ( mu(n/k1) + mu(n/k2) +mu(n/k3) +... )

对于mu( n/k1 )如果这个值为0的话就可以忽略了。所以可以枚举所有mu(x)不为0的数，然后 cnt(x * k) += mu(x)，就可以处理出所有的cnt(n)，然后用上面的分段的方法就可以解决了。具体见代码

/* **********************************************Author      : JayYeCreated Time: 2013-10-11 16:04:00File Name   : JayYe.cpp*********************************************** */#include <stdio.h>#include <string.h>#include <algorithm>using namespace std;typedef long long ll;const int maxn = 10000000 + 5;bool vis[maxn];int pri[666666], pnum, mu[maxn], cnt[maxn];void Mobius(int n) {    mu[1] = vis[1] = 1;    pnum = 0;    // 预处理出mu函数    for(int i = 2;i <= n; i++) {        if(!vis[i]) {            pri[pnum++] = i;            mu[i] = -1;        }        for(int j = 0;j < pnum; j++) {            if(i*pri[j] > n)    break;            vis[i * pri[j]] = 1;            if(i % pri[j] == 0) {                mu[i * pri[j]] = 0;                break;            }            mu[i * pri[j]] = -mu[i];        }    }    // 计算cnt(n)     for(int i = 1;i <= n; i++) if(mu[i]) {        for(int j = 0;j < pnum; j++) {            if(i * pri[j] > n)  break;            cnt[i * pri[j]] += mu[i];        }    }    // 处理出 cnt(n)的前缀和，因为要分段处理    for(int i = 1;i <= n; i++)  cnt[i] += cnt[i-1];}int main() {    Mobius(10000000);    int t, a, b;    scanf("%d", &t);    while(t--) {        scanf("%d%d", &a, &b);        if(a > b)   swap(a, b);        ll ans = 0;        for(int i = 1;i <= a; i++) {            int tmp1 = a/i, tmp2 = b/i;            int next = min(a/tmp1, b/tmp2); // 分段处理，找到下一段            ans += (ll)tmp1*tmp2*(cnt[next] - cnt[i-1]);            i = next;        }        printf("%lld\n", ans);    }    return 0;}