BZOJ(本校) 2525 公约数 - 莫比乌斯反演

题目描述
给出一个数N，求1 <= x, y <= N，且gcd(x, y)为素数的数对x, y的数量。
输入
一行一个数N。
输出
一个数表示答案
样例输入
4
样例输出
4
提示
【数据范围】
N ≤ 10000000
来源
HZOI

gcd(x,y)为质数,1<=x,y<=n
<=> gcd(x,y)=1 1<=x,y<=n/prime[i] prime[i]需要枚举。
新问题就很熟悉了， 假设 x《y ,想象：枚举y,求[1,y]中与y互质的x的个数 就是欧拉函数
<=>枚举prime[i],求pre_sum(phi(n/prime[i]))

1.绝对正确

#include<cstdio>
#include<cstring>
#define MAXN 10000000
#define MAXM 700000

int n,phi[MAXN+10],prime[MAXM+10],cntpr;
long long ans,pre[MAXN+10];
bool isprime[MAXN+10];
void GetPrime()
{
    memset(phi,0,sizeof phi);
    memset(isprime,0,sizeof isprime);
    cntpr=0;
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++cntpr]=i;
            phi[i]=i-1;
        }
        for(int j=1;prime[j]*i<=n&&j<=cntpr;j++){
            isprime[prime[j]*i]=true;
            if(i%prime[j]==0){
                phi[prime[j]*i]=phi[i]*prime[j];
                break;
            }
            phi[prime[j]*i]=phi[i]*(prime[j]-1);
        }
    }
}
int main()
{
    scanf("%d",&n);
    GetPrime();
    for(int i=1;i<=n;i++)
        pre[i]=pre[i-1]+phi[i];
    for(int i=1;i<=cntpr;i++){
        //<=>[1,n/prime[i]]中(x,y)互质的数对个数
        ans+=(pre[n/prime[i]]-1)*2+1;
    }
    printf("%lld\n",ans);
}

2.学习了一下标程：竟然还能这样做！

#include<complex>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<vector>
#include<cstring>
#include<iostream>
using namespace std;
const int MAXN = 10000000;
const int INF = 2000000000;
typedef long long LL;
int mu[MAXN +1000],prime[MAXN +1000],sum[MAXN +1000];
bool vis[MAXN +1000];
int P,a,b,d,c,k;
void init()
{
    int i,j;
    mu[1] = 1;
    for(i = 2;i <= MAXN;i++)
    {
        if(!vis[i])
        {
            prime[++P] = i;
            mu[i] = -1;
        }
        for(j = 1;1LL * i * prime[j] <= MAXN && j <= P;j++)
        {
            vis[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            mu[i * prime[j]] = -mu[i];
        }
    }
    for(i = 1;i <= MAXN;i++)
        sum[i] = sum[i - 1] + mu[i];
}
LL deal(int n,int m,int k)
{
    int i,last;
    n /= k;
    m /= k;
    if(n > m) swap(n,m);
    LL ans = 0;
    for(i = 1;i <= n;i = last + 1)
    {
        last = min(n/(n/i),m/(m/i));
        ans += 1LL * (n/i) * (m/i) * (sum[last] - sum[i-1]);
    }
    return ans;
}
int main()
{
    int T,n,i;
    init();
    scanf("%d",&n);
    LL ans = 0;
    for(i = 2;i <= n;i++)
        if(!vis[i])
            ans += deal(n,n,i);
    cout << ans << endl;
}

或者，这样：(by LiuJunhao 见他的blog)

#include<cstdio>
#include<algorithm>
using namespace std;
#define MAXN 10000000
int sum[MAXN+10],mu[MAXN+10],p[MAXN+10],pcnt,m,n,T;
long long ans;
bool f[MAXN+10];
void Read(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
void prepare(){
    int i,j;
    for(i=2;i<=MAXN;i++){
        if(!f[i])
            p[++pcnt]=i,mu[i]=-1,sum[i]=1;
        for(j=1;p[j]*i<=MAXN;j++){
            f[p[j]*i]=1;
            if(i%p[j]==0){
                sum[i*p[j]]=mu[i];
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
            sum[i*p[j]]=mu[i]-sum[i];
        }
        sum[i]+=sum[i-1];
    }
}
void solve(){
    int i,last,t=min(m,n);
    ans=0;
    for(i=1;i<=t;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans+=1ll*(sum[last]-sum[i-1])*(n/i)*(m/i);
    }
}
int main()
{
    Read(T);
    prepare();
    while(T--){
        Read(n),Read(m);
        solve();
        printf("%lld\n",ans);
    }
}

莫比乌斯反演做法的推导：
参考
参考

推到最后的表达式后，并且第一篇blog中用sum()代入表达式中。
sum()可以在筛法算Mobius时算出来，但O(n)的复杂度太高。
利用分块，降低复杂度至O(2*( sqrt(n) + sqrt(m) ))
具体是因为floor(n/k)很大一部分是一样的，具体看一下Liu Junhao的另一篇blog的分块证明，虽然丑，还是可以看得

转载于:https://www.cnblogs.com/katarinayuan/p/6572848.html