ACM学习历程—HYSBZ 2818 Gcd(欧拉函数 || 莫比乌斯反演)

Description

给定整数N，求1<=x,y<=N且Gcd(x,y)为素数的
数对(x,y)有多少对.

Input

一个整数N

Output

如题

Sample Input

Sample Output

Hint

对于样例(2,2),(2,4),(3,3),(4,2)

1<=N<=10^7

这个题目可以用欧拉函数或者莫比乌斯反演。

第一种欧拉函数：

因为gcd(x, y) = p，所以gcd(x/p, y/p) = 1。

不妨设y较大，那么就是求所有比y/p小的数k，phi(k)的和。phi(k)是欧拉函数，表示小于等于k，与k互质的数的个数。

然后每种乘2就得到了组数，但是需要减1，因为(1, 1)这组被计算了两次。

预处理所有phi(k)的前缀和就OK了，就能加速运算。

然后枚举质数p就OK了。

这个OJ我揣测是单组测试的，像下面预处理打表会T，需要对每个n进行一次打表。

代码：（欧拉函数）

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <set>
#include <map>
#include <queue>
#include <string>
#include <algorithm>
#define LL long long

using namespace std;

const int maxN = 1e7+5;
int n;
bool isprim[maxN];
int phi[maxN], prim[maxN], top;
LL s[maxN];

//埃氏筛法求素数
void isPrim()
{
    top = 0;
    memset(isprim, true, sizeof(isprim));
    isprim[0] = isprim[1] = false;//初始化
    for (LL i = 2; i < maxN; ++i)//筛法
    {
        if (isprim[i])
        {
            for (LL j = i*i; j < maxN; j += i)//上界太大可能会爆int
            {
                isprim[j] = false;
            }
            prim[top++] = i;
        }
    }
}

//***筛法求欧拉函数
void phi_table()
{
    memset(phi, 0, sizeof(phi));
    phi[1] = 1;
    s[0] = 0;
    s[1] = 1;
    for (int i = 2; i < maxN; i++)
    {
        if (!phi[i])
            for (LL j = i; j < maxN; j += i)//i如果太大可能会爆int
            {
                if (!phi[j])
                    phi[j] = j;
                phi[j] = phi[j]/i*(i-1);
            }
        s[i] = s[i-1]+phi[i];
    }
}

void work()
{
    LL ans = 0;
    for (int i = 0; i < top && prim[i] <= n; ++i)
    {
        ans += 2*s[n/prim[i]];
        ans--;
    }
    printf("%lld
", ans);
}

int main()
{
    //freopen("test.in", "r", stdin);
    isPrim();
    phi_table();
    while (scanf("%d", &n) != EOF)
    {
        work();
    }
    return 0;
}

View Code

第二种莫比乌斯反演：

设F(d)表示d | gcd(x, y) (1 <= x <= N, 1 <= y <= M)

f(d)表示d = gcd(x, y) (1 <= x <= N, 1 <= y <= M)

那么自然F(d) = sum(f(n)) (d | n)

且F(d) = (N/d)*(M/d)

根据莫比乌斯反演：

F(d) = sum(f(n)) (d|n) => f(d) = sum(u(n/d)*F(n)) (d|n)

所以f(d) = f(d) = sum(u(n/d)*F(n)) (d|n)

然后d取遍质数

则ans = sum(f(p)) = sum( sum(u(n/p)*F(n)) (p|n) )

外层sum上界为min(M, N), 内层sum也为min(M, N)

这样复杂度是min(M^2, N^2)的，时间上过不去。

此外就是这两个sum的位置是可以替换的，

因为整个式子的结果就是所有p|n的情况都要加上，

原式是对于一个p把所有n的情况加起来，最后求和。

同样可以，对于每一个n把所有p的情况加起来，最后求和。

那么式子可以改写为：ans = sum( F(n)*sum(u(n/p)) ) (p|n)

这样，如果能在时限内预处理出每个n对应的sum(u(n/p))，那么就能在min(M, N)时间复杂度内求解。

这里可以先把所有sum初始化为0，然后类似于筛法的做法，对于一个p，把所有p的倍数n的sum，加上u(n/p)。

这样预处理的复杂度要nlogn左右。不过常数有点大，本地需要1S多才能跑出来。（网上有一种O(n)的做法）

同样的，下面的打表也会T，需要对每个n进行一次打表。

代码：（莫比乌斯）

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <set>
#include <queue>
#include <vector>
#define LL long long

using namespace std;

const int maxN = 1e7+10;
int u[maxN], cnt;
LL prime[maxN], sum[maxN];
bool vis[maxN];

void mobius()
{
    memset(vis, false,sizeof(vis));
    u[1] = 1;
    cnt = 0;
    for(LL i = 2; i < maxN; i++)
    {
        if(!vis[i])
        {
            prime[cnt++] = i;
            u[i] = -1;
        }
        for(int j = 0; j < cnt && i*prime[j] < maxN; j++)
        {
            vis[i*prime[j]] = true;
            if(i%prime[j])
                u[i*prime[j]] = -u[i];
            else
            {
                u[i*prime[j]] = 0;
                break;
            }
        }
    }
}

int n;
LL ans;

void init()
{
    mobius();
    ans = 0;
    memset(sum, 0, sizeof(sum));
    for (int i = 0; i < cnt; ++i)
    {
        for (LL j = prime[i]; j < maxN; j += prime[i])
            sum[j] += u[j/prime[i]];
    }
}

void work()
{
    for (LL i = 1; i <= n; ++i)
        ans += (n/i)*(n/i)*sum[i];
    printf("%lld
", ans);
}

int main()
{
    //freopen("test.in", "r", stdin);
    init();
    while (scanf("%d", &n) != EOF)
        work();
    return 0;
}

View Code