筛法求莫比乌斯，欧拉函数

网上关于素数筛的资料很多，这里只是给出弱鸟整理的几个线性筛和应用。

最朴素的素数筛——埃拉托斯特尼筛法(Sieve of Eratosthenes) 
复杂度 Olognlogn

int primes[MAXN],tot=0;
bool isPrime[MAXN];

void getPrime()
{
    memset(isPrime,true,sizeof(isPrime));
    for(int i=2;i<MAXN;i++)
    {
        if(isPrime[i])
        {
            primes[++tot]=i;
            for(int j=i+i;j<MAXN;j+=i)
                isPrime[j]=false;
        }
    }
}


   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

欧拉素数筛 
复杂度 On 
其原理是每个合数都只会被一个质数筛去，因此为 On 线性筛法 
关键代码原理是 每个比 i 大的合数，必可以拆分为一个比 i 小的质数和另一个合数之积

int primes[MAXN],tot=0;
bool isPrime[MAXN];

void getPrime()
{
    memset(isPrime,true,sizeof(isPrime));
    for(int i=2;i<MAXN;i++)
    {
        if(isPrime[i])
            primes[++tot]=i;
        for(int j=1;j<=tot;j++)
        {
            if(i*primes[j]>=MAXN) break;
            isPrime[i*primes[j]]=false;
            if(i%primes[j]==0) break;//*****
        }
    }
}

   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

---

掌握线性筛之后，可以利用线性筛求积性函数

• 线性筛求欧拉函数
• 线性筛莫比乌斯函数

积性函数f(x)满足：f(ab)=f(a)f(b),a,b互质。

1.线性筛求欧拉函数

顺便复（预）习欧拉函数 
欧拉函数：wikipedia 
在数论中，对正整数n，欧拉函数 ϕ(n) 是小于或等于n的正整数中与n互质的数的数目。

定义式： 
欧拉函数ϕ(n)=n(1−1/p1)(1−1/p2)…(1−1/pk),p1…pk是n的k个不同的质因数。

观察到之前欧拉函数的筛法，在筛去一个合数（break）的同时能够得到它最小质因数 p 。假设 p 出现了 k 次，那么 f（p^k）和 f（n/p^k） 可以得到 f（n）

如果我们在筛去一个合数 i 的同时，记录它的每个质因数 pj 的次数，那么在转移至 i*pt 的同时，在筛去i·pt的同时得到它里面任意一个质因数的次数。

由定义式可得到以下结论：

ϕ(p)=p−1,p是质数。 
假设当前从ϕ(i),ϕ(pt)转移到ϕ(ipt)： 
1、如果pt是在ipt中第一次出现的话(也就是pt不整除i)，则ϕ(ipt)=ϕ(i)(pt−1) 
2、如果pt不是在ipt中第一次出现的话(也就是pt整除i)，则ϕ(ipt)=ϕ(i)pt

代码实现：

int tot=0;
int phi[maxn],prime[maxn];
bool isPrime[maxn];
void getphi(){
    mem(isPrime,true);//宏
    phi[1]=1;
    for(int i=2;i<=maxn;i++){
        if(isPrime[i]) prime[++tot]=i,phi[i]=i-1;//*
        for(int j=1;j<=tot;j++){
            if(i*prime[j]>maxn) break;
            isPrime[i*prime[j]]=false;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];//*
                break;
            }else phi[i*prime[j]]=phi[i]*(prime[j]-1);//*
        }
    }
}
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

---

2. 莫比乌斯函数

莫比乌斯函数μ(d)的定义如下 
(1)若d=1，那么μ(d)=1 
(2)若d=p1p2…pk(p1…pk均为互异质数)，那么μ(d)=(−1)^k 
(3)其他情况下，μ(d)=0

假设当前从μ(i),μ(pt)转移到μ(i·pt)，有以下几个式子便于我们筛出莫比乌斯函数： 
1、如果pt是在ipt中第一次出现的话(也就是pt不整除i)，则μ(i·pt)=−μ(i) 
2、如果pt不是在ipt中第一次出现的话(也就是pt整除i)，则μ(i·pt)=0

int tot=0;
int miu[maxn],prime[maxn];
bool isPrime[maxn];
void getmiu(){
    mem(isPrime,true);
    miu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(isPrime[i]) prime[++tot]=i,miu[i]=-1;
        for(int j=1;j<=tot;j++){
            if(i*prime[j]>maxn) break;
            isPrime[i*prime[j]]=false;
            if(i%prime[j]==0){
                miu[i*prime[j]]=0;//略有不同
                break;
            }else miu[i*prime[j]]=-1*miu[i];
        }
    }
}
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

3.莫比乌斯反演

F(n),f(n) 是定义在非负整数集合上的两个函数，并且满足条件 
 
那么有结论: 

有些题目会要求求出f(n)的值，但是往往F(n)比f(n)更好求。如果我们难以求出f(n)的值，但是能找出一个F(n)函数，并且可以很轻松地求出F(n)的值的话，就应该使用莫比乌斯反演。

来个例题吧！ 
给一个正整数n，其中n<=10^7，求使得gcd(x,y)为素数的的(x,y)的个数，1<=x,y<=n。

分析： 
对于本题，因为是使得gcd(x,y)为素数，所以必然要枚举小于等于n的素数，那么对于每一个素数pi，只需要求在区间 [1,n/pi] 中，满足有序对互质的对数。 
因为gcd(x,y)=pi –> gcd(x/pi,y/pi)=1 
设 f(pi)=(gcd(x,y)=pi)的种类数，F(pi) =（pi | gcd(x,y)）的种类数 
可得

F(i)=∑d|if(d)

又由F(i)的含义得

F(i)=floor(n/i)*floor(n/i)
   
   
1
   
   
1

则利用莫比乌斯反演，得 

f(i)=∑d|nμ(d)F(i/d)

再枚举 i ，范围为 [1,n/pi]

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef set<int> Set;
typedef vector<int> Vec;
#define fast ios_base::sync_with_stdio(false)
#define mem(s,t) memset(s,t,sizeof(s))
#define D(v) cout<<#v<<" "<<v<<endl
#define mod 1000000007
#define inf 0x3f3f3f3f
const int maxn =1e7+10;
ll tot=0,ans=0,n;
int miu[maxn],prime[1000005];
bool isPrime[maxn];
void getmiu(){
    mem(isPrime,true);
    miu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(isPrime[i]) prime[++tot]=i,miu[i]=-1;
        for(int j=1;j<=tot;j++){
            if(i*prime[j]>maxn) break;
            isPrime[i*prime[j]]=false;
            if(i%prime[j]==0){
                miu[i*prime[j]]=0;
                break;
            }else miu[i*prime[j]]=-1*miu[i];
        }
    }
}

void solve(){
    scanf("%d",&n);
    for(int k=1;prime[k]<=n;k++){
        int t=n/prime[k];
        for(int i=1;i<=t;i++){
            ans+=(ll)(t/i)*(t/i)*miu[i];
        }
    }
    printf("%lld\n",ans);
}

int main(){
#ifdef LOCAL
    freopen("in.txt","r",stdin);
    freopen("out.txt","w",stdout);
#endif
    mem(prime,1);
    getmiu();
    solve();
    return 0;
}

   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
   
   

    
    
   
   
   
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52