先拜大牛。感謝賈志鵬嚴謹的思維。以及簡單清晰的論文描述。
一定要結合論文看。我只是提出我覺得關鍵的部分。論文在網上隨處可見。賈志鵬線性篩。
開頭兩種線性篩的比較。
一種是傳統的線性篩。時間復雜度為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); }
問題三：給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); }
問題五：給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]); } }
再考慮一個問題：篩一個數的所擁有的不重復的素因子之和。比如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]); } }
擴展問題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)); } }
實際上求互質數和有 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]); }
歡迎點個贊~