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的分块证明，虽然丑，还是可以看得