Min_25 篩這個東西，徹底理解花了我很長的時間，因此寫點東西來記錄一些本身的理解。php
它能作什麼
對於某個數論函數 f，若是知足如下幾個條件，那麼它就能夠用 Min_25 篩來快速求出這個函數的前綴和。ios
1.它是一個積性函數c++
2.對於一個質數 p ，f(p) 的表達式必須是一個項數比較小的多項式。即 
數組
3.對於一個質數 p ，f(pk的表達式必須能夠由 f(p) 快速獲得。微信
看起來這些條件比較緊，其實仔細想一想許多函數仍是知足這些性質的。接下來討論它的實現過程。函數
實現過程
Min_25 篩的核心思想就是對於 ≤n 的數中，除掉質數以後，每一個數的最小質因子大小不超過 √n 。即相似埃氏篩的方式，從每一個質因子篩掉一部分的貢獻。整個過程能夠分紅三部分，即質數的部分的貢獻之和與合數的部分的貢獻和，以及 1 的特判。spa
記號約定
這裏用 P 表示質數集合，Pi 表示第 i 小的質數。.net
質數部分
首先咱們不妨設 f(p)=pk 。對於更加通常的項數小的多項式形式，咱們能夠分別記每個項的貢獻，而後分別相加。注意到這裏 f(p) 是一個徹底積性函數。對於 p 爲合數的狀況，爲了方便，咱們先姑且一樣定義 f(p)=pk ，到最後你會明白這樣定義的緣由。3d
咱們再設g(n,i) 表示 ≤n 的全部數 x 中，x 的最小質因子 >Pi 或者 x 自己爲質數的全部 f(x) 之和，這樣 g(n,|P|) 就是全部質數的貢獻之和了。而 g(n,0) 的答案也十分顯然，即 code
。咱們只須要考慮，如何從 g(∗,i−1) 如何獲得g(∗,i) 。
考慮咱們從i−1 推到 i 的過程當中，會減小哪一個部分的貢獻，容易發現最小質因子剛好爲 Pi 的貢獻被減掉了。若是 Pi2>n ，那麼這一部分是沒有貢獻的，即 g(n,i)=g(n,i−1)。
對於 Pi2≤n 的狀況，咱們須要把全部 Pi 的倍數，且知足其最小質因子除去 Pi 以外 >Pi−1 的貢獻所有減掉。若是直接減去
你會發現貢獻是減多了的，緣由就是 g 函數的值還包含了全部質數的貢獻之和，仔細分析一下，能夠發現，對於
的貢獻在 j 這個地方已經減過了，因此還要加回來，即
式子彙總起來也就是
回頭來，咱們再來看最初咱們提出的一些問題。
咱們最開始爲何要定義，對於一個合數 p ，f(p)=pk 呢？無非就是方便計算，方便利用徹底積性函數的性質。你發現到最後，合數的貢獻仍是沒有算進去，因此對於合數的定義是無所謂的。
咱們到底在哪裏利用了 f(p) 是一個徹底積性函數的性質呢？這個性質是在上式中第二種狀況中利用到的。注意到咱們在減貢獻的時候，並無枚舉 Pi 的次數，這也就意味着 中可能包含了某些數仍然是 Pi 的倍數。因此這一步須要有一個對於全部質數 p, 的條件。因爲以前咱們假定了 f(p) 對於 p 爲合數的狀況也知足 f(p)=pk ，所以這個條件是獲得知足了的。
遞歸實現很顯然。考慮一下非遞歸的實現。首先對於這個步驟，咱們有用的 nn 只有全部不一樣的 ⌊ni⌋,(1≤i≤n)⌊ni⌋,(1≤i≤n)，因此能夠離散掉這些值，用一個大小爲 2n−−√2n 的數組存下來。因爲這個遞推是一層一層推的，因此實現上能夠不用開二維數組，能夠直接離散以後用一個一維數組存下來，須要注意枚舉順序。
合數部分
考慮咱們在求質數部分的時候，除了質數的答案之和以外，還求出了些什麼東西。事實上，對於全部的 ⌊n/i⌋,(1≤i≤n) ，咱們求出了 g([n/i⌋,|P|) 的值。
咱們設 S(n,i) 表示 ≤n 的全部數 x 中，x 的最小質因子≥Pi 的 f(x) 之和。注意這裏的 f(x) 在 x 爲合數的時候再也不沿用質數部分時的定義，也注意這裏的定義與 g 的區別。
接下來咱們先算上全部知足條件的質數貢獻之和，即。對於合數，咱們直接枚舉其最小質因子以及質因子的個數，直接遞歸計算。注意在這裏形如pk,(p∈P) 的貢獻並無算進去，因此還要單獨加一下。式子的形式以下：
直接對着式子看，應該仍是很好理解的。對這個式子，直接遞歸暴力計算便可，跑得會比第一部分快。
至此 Min_25 篩的全過程就講完了。
附Min25篩求質數和代碼
#include <math.h>#include <stdio.h>#include <assert.h>#define LINT long long
inline LINT V2IDX(LINT v, LINT N, LINT Ndr, LINT nv) { return v >= Ndr ? (N/v - 1) : (nv - v);}
LINT primesum(LINT N) { LINT *S; LINT *V;
LINT r = (LINT)sqrt(N); LINT Ndr = N/r;
assert(r*r <= N and (r+1)*(r+1) > N);
LINT nv = r + Ndr - 1;
V = new LINT[nv]; S = new LINT[nv];
for (LINT i=0; i<r; i++) { V[i] = N/(i+1); } for (LINT i=r; i<nv; i++) { V[i] = V[i-1] - 1; }
for (LINT i=0; i<nv; i++) { S[i] = V[i] * (V[i] + 1) / 2 - 1; }
for (LINT p=2; p<=r; p++) { if (S[nv-p] > S[nv-p+1]) { LINT sp = S[nv-p+1]; LINT p2 = p*p; for (LINT i=0; i<nv; i++) { if (V[i] >= p2) { S[i] -= p * (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp); } else { break; } } } }
return S[0];}
int main() { LINT N = 1000000000; printf("%lld\n", primesum(N));}
附加一道剛作過的題
http://acm.hdu.edu.cn/contests/contest_showproblem.php?pid=1002&cid=909
題目代碼以下：
#include <iostream>#include <string.h>#include <stdlib.h>#include <stdio.h>#include <math.h>#include <assert.h>#include <bits/stdc++.h>using namespace std;const int INF = 0x3f3f3f3f;const int N = 1000010;typedef long long ll;const int mod =1e9+7;ll n,k;
long long ksm(long long x,long long y,long long mod){ long long ans=1; while(y){ if(y&1) ans=ans*x%mod; y>>=1; x=x*x%mod; } return ans;}
bool ip(ll x){ for(int i=2;i<sqrt(x);i++) { if(x%i==0) return false; } return true;}
ll zsh(ll n){ ll ans=0; for(int i=2;i<=n;i++) { if(ip(i)) ans+=i; } return ans;}
inline ll jf(ll x, ll y, ll p, ll q) { return x >= p ? (y/x - 1) : (q - x);}
ll solution(ll N) { ll *S; ll *V; ll r = (ll)sqrt(N); ll qq = N/r; assert(r*r <= N and (r+1)*(r+1) > N); ll maxn = r + qq - 1; V = new ll[maxn]; S = new ll[maxn]; for (ll i=0; i<r; i++) V[i] = N/(i+1); for (ll i=r; i<maxn; i++) V[i] = V[i-1] - 1; for (ll i=0; i<maxn; i++) S[i] = V[i] * (V[i] + 1) / 2 - 1; for (ll p=2; p<=r; p++) { if (S[maxn-p] > S[maxn-p+1]) { ll sp = S[maxn-p+1]; ll p2 = p*p; for (ll i=0; i<maxn; i++) { if (V[i] >= p2) { S[i] -= p * (S[jf(V[i]/p, N, qq, maxn)] - sp); } else { break; } } } } return S[0];}
int main(){ int t; cin>>t; while(t--) { scanf("%lld%lld",&n,&k); ll sum=0; cout<<(solution(n+1)%k-4+(n+3)%k*n%k*ksm(2,k-2,k)%k+k)%k<<endl; } return 0;}
本文分享自微信公衆號 - WHICH工做室（which_cn）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。