讀賈志鵬線性篩有感 (莫比烏斯函數的應用)

先拜大牛。感謝賈志鵬嚴謹的思惟。以及簡單清晰的論文描述。

必定要結合論文看。我只是提出我以爲關鍵的部分。論文在網上隨處可見。賈志鵬線性篩。

 

開頭兩種線性篩的比較。

一種是傳統的線性篩。時間複雜度爲N*log(log(N))。

另一種是優化了合數的篩法。文中稱做Euler線性篩。

其優化的地方。

舉個例子：合數6。 是2的倍數也是3的倍數。 當你用傳統的篩法的時候在遍歷2的倍數的時候會遍歷到6。遍歷3的倍數的時候一樣也會遍歷到6。

而另一種只會篩出6爲2的倍數。3就不會篩6了。

另外我的認爲篩法二有一個很重要的思想。當i爲合數的時候。其實腦海裏不認爲是合數。而是素數的乘積。這樣就能比較直觀地肯定這個算法的正確性了。

 

積性函數。

分爲徹底積性和條件積性。

咱們最喜歡的積性。大概就是互素積性了。由於知足互素積性的話。根據算術基本定理。就可以簡單作到推廣到任意實數。

f(1) = 1 。 這個在咱們高中數學題。抽象函數。就已經能簡單知道了。

 

歐拉函數。就再也不談了。包括其線性篩的那一步相當重要的證實。也在我其它博文提到過了。

其 歐拉定理和費馬小定理的做用。我得再多補充一點。

以及互質數和。 n的互質數和爲 n*φ(n)/2.

 

莫比烏斯函數和容斥定理的關係。

能夠發現莫比烏斯函數其實就是容斥定理的映射通常。

莫比烏斯函數 是咱們再熟悉不過的了。不熟悉能夠看這裏。

首先看 (-1)^r m = p1p2p3p4p5pr 其實就是在模擬容斥定理。

假如一但不是素數。那就爲0.

 

兩個函數的線性篩。這實際上是咱們處理問題的基本。這裏須要講的是。不必定只有積性函數才能夠用這種篩法。

只要你能找到f(kn) n整除k 和不整除的兩個時刻所對應的遞推式。這個在擴展問題中會出現。

問題一：求1~N對質數P的乘法逆元。

關於f(n)爲徹底積性函數。根據同餘定理能夠簡單得到。要證實的話。減法證同餘便可。

P = nt + k

n'  ≡ n*(t^2)*f(k)^2 (mod P)

這個證實過程很漂亮(很佩服這麼順暢,思惟這麼清晰)。也是根據同餘定理。還有逆元的性質。就能推理的。

這個問題的意義。能夠求N!的 mod P 的逆元了。逆元仍是頗有用的。由於畢竟除法並無特別好的同餘式。（依稀還記得那兩個。）

問題二：給T組N,M.依次求出的值.(N,M<=10^6,T<=10^3)

求解gcd(a,b).把gcd(a,b)當作n.再經過歐拉函數和式。推導過程以下。

第二個等式是用d來看待式子的方法來化簡和式的。

以後再窮舉d便可。

 

 

#include<stdio.h>
#include<string.h>
#define    N 100

bool mark[N+5];
int prime[N+5];
int num;
int euler[N+5];
int Min(int a,int b){return a>b?a:b;}
void Euler()
{
    int i,j;
    num = 0;
    memset(euler,0,sizeof(euler));
    memset(prime,0,sizeof(prime));
    memset(mark,0,sizeof(mark));
    euler[1] = 1; // multiply function
    for(i=2;i<=N;i++)
    {
        if(!mark[i])
        {    
            prime[num++] = i;
            euler[i] = i-1;
        }
        for(j=0;j<num;j++)
        {
            if(i*prime[j]>N){break;}
            mark[i*prime[j]] = 1;
            if(i%prime[j]==0)
            {
                euler[i*prime[j]] = euler[i]*prime[j];
                break;
            }
            else
            {
                euler[i*prime[j]] = euler[prime[j]]*euler[i];
            }
        }
    }
}

int main()
{
    int i;
    int M1,M2;
    Euler();
    for(i=0;i<num;i++)printf("%d ",prime[i]);
    printf("\n");
    for(i=1;i<=N;i++)printf("%d ",euler[i]);
    printf("\n");
    M1 = 2;
    M2 = 3;
    int sum = 0;
    int min = Min(M1,M2);
    for(i=1;i<=min;i++)
    {
        sum += euler[i]*(M1/i)*(M2/i);
    }
    printf("%d\n",sum);
}

the second problem test

 

 

 

問題三：給T組N,M.依次求出的值.(N,M<=10^6,T<=10^3)

在證實以前，先證實如下式子。

 

 

問題的解決推導。

第一個等式。lcm(a,b) = ab/gcd(a,b).

第二個等式。令d=gcd(a,b)。

第三個等式。轉化爲d的視角。（這個手法常常有）。

第四個等式。轉化爲莫比烏斯函數。

第五個等式。利用上述的等式來轉化。注意d和d'

第六個等式。論文中提到的有趣的化簡性質。

第七個等式。實際上是d = dd'換元。不過有點老奸巨猾啊。幹嗎不設個T = dd'。這個我糾結了半天。

以後就是如論文中介紹的。g(d) 爲積性函數。線性篩之。

整體上算法仍是N的。

問題四：給T組N,M.依次求出的值.(N,M<=10^6,T<=10^3)

實質上就是求　其中x和y互素。的對數。

咱們是時候須要有和式化成的思想了。[gcd(a,b)=1]真是漂亮的莫比烏斯函數的和式的結果。

第一個等式：莫比烏斯函數擴寫

第二個等式：gcd(a,b)=p -> gcd(a/p,b/p)=1問題轉換。

第三個等式：一個和式的處理手段。

第四個等式：很常見的。

 

#include<stdio.h>
#include<string.h>
#define    N 100

bool mark[N+5];
int prime[N+5];
int num;
int mobi[N+5];
int Min(int a,int b){return a>b?a:b;}
void Mobi()
{
    int i,j;
    num = 0;
    memset(mobi,0,sizeof(mobi));
    memset(prime,0,sizeof(prime));
    memset(mark,0,sizeof(mark));
    mobi[1] = 1; // multiply function
    for(i=2;i<=N;i++)
    {
        if(!mark[i])
        {    
            prime[num++] = i;
            mobi[i] = -1;
        }
        for(j=0;j<num;j++)
        {
            if(i*prime[j]>N){break;}
            mark[i*prime[j]] = 1;
            if(i%prime[j]==0)
            {
                mobi[i*prime[j]] = 0;
                break;
            }
            else
            {
                mobi[i*prime[j]] = mobi[prime[j]]*mobi[i];
            }
        }
    }
}

int main()
{
    int i;
    int M1,M2;
    Mobi();
    for(i=0;i<num;i++)printf("%d ",prime[i]);
    printf("\n");
    for(i=1;i<=N;i++)printf("%d ",mobi[i]);
    printf("\n");
    M1 = 2;
    M2 = 3;
    int sum = 0;
    int min = Min(M1,M2);
    for(i=1;i<=min;i++)
    {
        sum += mobi[i]*(M1/i)*(M2/i);
    }
    printf("%d\n",sum);
}

Test problem forth

 

 

 

 

問題五：給T組N.依次求出的值.(N<=10^6,T<=10^3)

 

其實根據問題三.能夠直接得到該化簡出來的式子的。

而後解法和問題三同樣。

可是論文上尋找積性f(n)直接篩出答案。

首先佩服一下其思惟的緊密。一個變量啊。就尋找積性函數。這個轉化真是清晰而又巧。

畫個圖就能知道 -n 是用來去重複的統計的。

.

f(n)是積性的。具體證實如論文上解釋。

第一個等式:d = gcd(n,i)

第二個等式:k = i/d.且所有都除以d.gcd(a,b)=d轉化成求互素(gcd(a,b)=1)的問題。

第三個等式:令d=n/d。是對應的。  其實在第二個等式就能看出是歐拉函數求約數和問題了。

第四個等式:不解釋了吧。

第五個等式:手算一下容易得。

歐拉函數求互質數和的函數是積性函數

有一道題。就是利用這個。後會介紹。

見到積性函數咱們如今應該是very happy的。

擴展問題1: T組N.依次求出的值.(N<=10^6,T<=10^3)

借鑑了賈志鵬上面全部問題的證實。這個是我本身寫的擴展證實。不免有錯誤。見諒。還望留言提醒。

我以爲這樣的證實是很是輕快明瞭的。而後網上還有流行一種。用莫比烏斯反演的另一種表示式的。也是很是神奇的。

不過。那個反演我尚未證實過。不過仍是借鑑了其下半部分的設T。（也是這個設T點醒了我。賈志鵬第3個問題的證實的最後一步。）

這裏並不能由於p是素數。而n/p不必定是素數。因此並非對稱的。(若是看過具體數學就能很快明白了。)

 

處理。

分類對於 g(kx) .有

 g(kx)=μ(x)                 k|p

 g(kx)=μ(x)-g(x)          k!|p

 

結合莫比烏斯函數。能夠知道分類成立:

咱們能夠借這個 而且借用以前兩個積性函數的篩法 來篩 g(n)。

這是明顯可行的。也就是說。咱們不須要函數必須是積性的才能去篩。

咱們只須要找到g(kx)是由g(x)得到的。或者是在g(x)以前就篩掉的值得到的。eg:g(x-1) (篩法老是從小到大。)

更甚的是。咱們只須要得到大值和小值的關係！就能夠篩法。可是該篩法。是創建在素數表達式之上的。

 

這段闡述也許很混亂。可是我也只能描述個大概的我的體會。理解不理解不要緊。

 給你個例子。篩一個數的素因子之和。

對於上述的篩法.

讓F[n] 爲n的素因子之和。

F[i*prime[j]] = F[i]+prime[j]                            i\prime[j]時。

F[i*prime[j]] = F[i]+prime[j]                             i!\prime[j]時。 

兩種狀況是同樣的。緣由顯而易見。不過咱們仍是得判斷。由於i\prime[j]的時候咱們更新後能夠直接break;

再考慮一個問題:篩一個數的所擁有的素因子之和。 好比12 爲2*2*3 咱們只計算2+3.

那麼有：

F[i*prime[j]] = F[i]                                             i\prime[j]時。

F[i*prime[j]] = F[i]+prime[j]                             i!\prime[j]時。 

對於這個問題的code.

#include<stdio.h>
#include<string.h>
#define N 100

int num,prime[N+5],f[N+5];
bool mark[N+5];
void Init()
{
    int i,j;
    num = 0;
    f[1] = 0;
    for(i=2;i<=N;i++)
    {
        if(!mark[i])
        {
            prime[num++] = i;
            f[i] = i;
        }

        for(j=0;j<num;j++)
        {
            if(i*prime[j]>N){break;}
            mark[i*prime[j]] = 1;
            if(i%prime[j]==0)
            {
                f[i*prime[j]] = f[i];
                break;
            }
            else
            {
                f[i*prime[j]] = f[i]+prime[j];    
            }
        }

    }
}

int main()
{
    int i;
    Init();
    for(i=1;i<=N;i++)
    {
        printf("%d = %d \n",i,f[i]);
    }
}

Problem test

再考慮一個問題：篩一個數的所擁有的不重複的素因子之和。好比12 爲2*2*3 咱們只計算3

那麼有：

 i\prime[j]時。

　　狀況1: （i/prime[j]）\prime[j]時.

　　F[i*prime[j]] = F[i]

　　狀況2： （i/prime[j]）！\prime[j]時.

　　FF[i*prime[j]] = F[i]-prime[j].

 

 i!\prime[j]時。 

F[i*prime[j]] = F[i]+prime[j]                           

 

#include<stdio.h>
#include<string.h>
#define N 100

int num,prime[N+5],f[N+5];
bool mark[N+5];
int Max(int a,int b)
{
    return a>b?a:b;
}
void Init()
{
    int i,j;
    num = 0;
    f[1] = 0;
    for(i=2;i<=N;i++)
    {
        if(!mark[i])
        {
            prime[num++] = i;
            f[i] = i;
        }

        for(j=0;j<num;j++)
        {
            if(i*prime[j]>N){break;}
            mark[i*prime[j]] = 1;
            if(i%prime[j]==0)
            {
                if((i/prime[j])%prime[j]==0)
                {
                f[i*prime[j]] = f[i];
                }
                else
                {
                f[i*prime[j]] = f[i] - prime[j];
                }

                break;
            }
            else
            {
                f[i*prime[j]] = f[i]+prime[j];    
            }
        }

    }
}

int main()
{
    int i;
    Init();
    for(i=1;i<=N;i++)
    {
        printf("%d = %d \n",i,f[i]);
    }
}

Problem test

 

 擴展問題1: T組N.求1~N範圍上與N互素的數的和。

 

值得一提的是推導到最後的。按照以往的手段彷佛沒有繼續下去的可能了。（可是若是你仔細觀察的話。能夠發現 n/k 不須要取底符號。那麼就能提取出一個n的因子）

Code:

#include<stdio.h>
#include<string.h>
#define N 100
int num;
int prime[N+5];
int mobius[N+5];
bool mark[N+5];

void Mobius()
{
    int i,j;
    num = 0;
    mobius[1] = 1;

    for(i=2;i<=N;i++)
    {
        if(!mark[i])
        {
            prime[num++] = i;
            mobius[i] = -1;
        }
        for(j=0;j<num;j++)
        {
            if(i*prime[j]>N){break;}
            mark[i*prime[j]] = 1;
            if(i%prime[j]==0)
            {
                mobius[i*prime[j]] = 0;
            }
            else
            {
                mobius[i*prime[j]] = -mobius[i];
            }
        }
    }
}
int solve(int n)
{
    int i,r;
    r = 0;
    for(i=1;i<=n;i++)
    {
        if(n%i==0)
        {
            r += mobius[i]*i*(n/i)*(n/i+1);
        }
    }
    r /= 2;
    return r;
}
int main()
{
    int i,n;
    Mobius();
    while(scanf("%d",&n)!=EOF)
    {
        printf("%d = %d\n",n,solve(n));
    }
}

Problem two

 

 實際上求互質數和有 n*φ(n)/2 。

用莫比烏斯函數表示

 上面公式得證。十分感謝yzq986的留言。告訴了我後續的解法！！！~~~

若是咱們直接用n*φ(n)/2。該函數咱們是能夠直接篩出來的。

 對於互質數咱們探討得較多了。個數（歐拉函數）。互素數和。就是以上的。

那麼對於約數呢？另外開一篇隨筆去探討這個問題。

 

論文上的一個優化:

論文上sqrt的優化具體原理論文已經給得很清楚了。

即存在 a/x = a/(x+k)  這個是取整除法

稍微講述一下代碼的構造。

咱們預處理出目標函數以後。再預處理出其前綴和用sum數組保存.經過如下代碼進行結果的處理。便可。

 for(int i=1,last;i<=n;i=last+1)  
 {  
    last = min(n/(n/i),m/(m/i));  
    ans += (n/i)*(m/i)*(sum[last]-sum[i-1]);  
  }  

優化

 

歡迎點個贊~