洲閣篩 & min_25篩學習筆記

洲閣篩

給定一個積性函數$F(n)$，求$\sum_{i = 1}^{n}F(n)$。而且$F(n)$知足在素數和素數次冪的時候易於計算。

顯然有：

$\sum_{i = 1}^{n} F(n) = \sum_{i = 1}^{\sqrt{n}}F(i) \left(\sum_{\sqrt{n} < p\leqslant n/i, p\ is\ a\ prime} F(p) \right) + \sum_{i = 1, i\ has\ no\ prime\ factor\ greater\ than\ \sqrt{n}}^{n} F(i)$。

Part 1

第一部分是要求出$\sum_{\sqrt{n}<p \leqslant n/i, p\ is\ a\ prime} F(p)$

假設小於等於$\sqrt{n}$的素數從小到大分別是$P_1, P_2, \cdots, P_{m}$。

設$f_{i, j}$表示$[1, j]$中和前$i$個素數互質的數的函數值和。

當$i = 0$時候想辦法計算。

考慮當$i > 0$的時候的轉移。考慮從$f_{i - 1, j}$中須要刪掉最小質因子等於$P_i$的數的貢獻。

因此有轉移式子$f_{i, j} = f_{i - 1, j} - F(P_i) f_{i - 1, \left \lfloor j / P_i \right \rfloor}$

若是$F(n)$自己不是徹底積性函數，須要把它拆成若干積性函數的和。注意這一部分只用保證素數處取值相同便可。

注意到第一維約有$O\left (\frac{\sqrt{n}}{\log n} \right)$種取值，第二維有$O\left (\sqrt{n}\right)$種取值。暴力計算的複雜度是$O\left (\frac{n}{\log n} \right )$

優化

考慮當$j < P_{i+ 1}$的時候，$f_{i, j} = 1$。考慮把條件放縮成$j < P_i$。

所以，當$P_i \leqslant j < P_{i}^2$的時候$f_{i, j} = f_{i - 1, j} - F(P_i)$。

因此作法是：

• 維護一個$F(P_i)$的前綴和
• 考慮轉移時的$f_{i - 1, \left \lfloor j / P_i \right \rfloor}$
  • 當$j \geqslant P_{i}^2$的時候暴力計算
  • 當$P_i \leqslant j < P_{i}^2$的時候，須要這樣的狀態的時候減去一段和就好了。
  • 當$j < P_i$的時候，因爲$f_{i, 0} = 0$，因此這一部分的值不會改變。

因此對於每一個$j$，有效的狀態只有$O\left (\frac{\sqrt{j}}{\log j} \right )$。

時間複雜度爲：$\sum_{i = 1}^{\sqrt{n}} O \left (\frac{\sqrt{i}}{\log i} \right ) + \sum_{i = 1}^{\sqrt{n}} O\left ( \frac{\sqrt{n/i}}{\log {n/i}}\right)$

因爲第一部分顯然小於第二部分，因此只用考慮第二部分的和。

因此有：$\sum_{i = 1}^{\sqrt{n}}O \left (\frac{\sqrt{n/i}}{\log n - \log i}\right)$

因爲$\log i \leqslant \log \sqrt{n} = \frac{1}{2}\log n$

因此：$\sum_{i = 1}^{\sqrt{n}}O \left (\frac{\sqrt{n/i}}{\frac{1}{2}\log n}\right) = \sum_{i = 1}^{\sqrt{n}}O\left (\frac{\sqrt{n / i}}{\log n}\right) = O\left ( \frac{n^{3/4}}{\log n}\right)$

Part 2

這一部分是要求出$\sum_{i = 1, i\ has\ no\ prime\ factor\ greater\ than\ \sqrt{n}}^{n} F(i)$

假設小於等於$\sqrt{n}$的素數從小到大分別是$P_1, P_2, \cdots, P_{m}$。

設$g_{i, j}$表示在$[1, j]$中，質因子只包含$P_{i}, P_{i + 1}, \cdots, P_{m}$的數的函數值的和。

轉移考慮枚舉當前質數$P_i$的指數。因此有：$g_{i, j} = g_{i - 1, j} + \sum_{e = 1}^{\log_{P_i} j} F(P_{i}^{e}) g_{i - 1, \lfloor j / P_i^e\rfloor}$

對於時間複雜度，考慮對於每個$j$的轉移的枚舉量，能夠獲得大概是$O(\frac{\sqrt{j}}{\log j})$。<del>(不會積分，告辭)</del> 

和Part 1同樣，獲得暴力計算這一部分的複雜度是$O(\frac{n}{\log n})$。

優化

由於質因子是從大到小枚舉的，因此當$j < P_{i}$的時候，$g_{i, j} = 1$。

因此當知足$P_i \leqslant j < P_{i}^2$，知足$g_{i, j} = g_{i - 1, j} + F(P_i)$。

和Part 1同樣，分紅三段，一段暴力轉移，一段須要時前綴和，一段不須要維護。

複雜度也是一樣的計算方法

  1 /**
  2  * SPOJ
  3  * Problem#DIVCNT3
  4  * Accepted
  5  * Time: 35190ms
  6  * Memory: 33792k
  7  */
  8 #include <bits/stdc++.h>
  9 using namespace std;
 10 typedef bool boolean;
 11 
 12 template <typename T>
 13 void pfill(T* pst, const T* ped, T val) {
 14     for ( ; pst != ped; *(pst++) = val);
 15 }
 16 
 17 #define ll long long
 18 
 19 typedef class Euler {
 20     public:
 21         int num;
 22         boolean *vis;
 23         int *pri, *d3, *mp;
 24         
 25         Euler(int n) : num(0) {
 26             d3 = new int[(n + 1)];
 27             mp = new int[(n + 1)];
 28             pri = new int[(n + 1)];
 29             vis = new boolean[(n + 1)];
 30             pfill(vis, vis + n + 1, false);
 31             d3[1] = 1;
 32             for (int i = 2; i <= n; i++) {
 33                 if (!vis[i]) {
 34                     d3[i] = 4, pri[num++] = i, mp[i] = 1;
 35                 }
 36                 for (int j = 0, _ = n / i, x; j < num && pri[j] <= _; j++) {
 37                     vis[x = pri[j] * i] = true;
 38                     if (!(i % pri[j])) {
 39                         d3[x] = (d3[i / mp[i]] + 3) * d3[mp[i]];
 40                         mp[x] = mp[i];
 41                         break;
 42                     } else {
 43                         mp[x] = i, d3[x] = d3[i] << 2;
 44                     }
 45                 } 
 46             }
 47         }
 48         ~Euler() {
 49             delete[] mp;
 50             delete[] vis;
 51             delete[] d3;
 52             delete[] pri;
 53         }
 54 } Euler;
 55 
 56 const int N = 330000;
 57 
 58 Euler *_euler;
 59 
 60 int T;
 61 vector<ll> ns;
 62 
 63 ll maxn, D;
 64 
 65 int cnt_prime;
 66 int *pri, *d3;
 67 
 68 inline void init() {
 69     scanf("%d", &T);
 70     for (ll x; T--; ) {
 71         scanf("%lld", &x);
 72         ns.push_back(x);
 73         maxn = max(maxn, x);
 74     }
 75     D = sqrt(maxn + 0.5) + 1;
 76     while (D * 1ll * D <= maxn)
 77         D++;
 78     _euler = new Euler(D + 23);
 79     cnt_prime = _euler->num;
 80     pri = _euler->pri, d3 = _euler->d3;
 81 }
 82 
 83 ll n;
 84 int cnt_P[N];
 85 int range0[N], range1[N];
 86 ll f0[N], f1[N], g0[N], g1[N];
 87 
 88 void precalc() {
 89     for (int i = 1; i < D; i++) {
 90         range0[i] = range0[i - 1];
 91         while (pri[range0[i]] * 1ll * pri[range0[i]] <= i) {
 92             range0[i]++;
 93         }
 94     }
 95     for (int i = D - 1; i; i--) {
 96         ll y = n / i;
 97         range1[i] = range1[i + 1];
 98         while (pri[range1[i]] * 1ll * pri[range1[i]] <= y) {
 99             range1[i]++;
100         }
101     }
102     pfill(cnt_P, cnt_P + D, 0);
103     for (int i = 0; pri[i] < D; i++) {
104         cnt_P[pri[i]]++;
105     }
106     for (int i = 1; i < D; i++) {
107         cnt_P[i] += cnt_P[i - 1];
108     }
109 }
110 
111 int nump;
112 ll get_value_f(int i, ll j) {
113     if (j >= D) {
114         int rj = n / j;
115         return f1[rj] + (max(0, nump - max(range1[rj], i)) << 2);
116     }
117     return f0[j] + (max(0, cnt_P[j] - max(range0[j], i)) << 2);
118 }
119 
120 void calculate_f() {
121     for (nump = 0; pri[nump] < D; nump++);
122     for (int i = 1; i < D; i++) {
123         f0[i] = 1, f1[i] = 1;
124     }
125     for (int i = nump - 1; ~i; i--) {
126         for (int j = 1; j < D && i < range1[j]; j++) {
127             ll m = n / j, coef = 1;
128             while (m /= pri[i]) {
129                 f1[j] += (coef += 3) * get_value_f(i + 1, m);
130             }
131         }
132         for (int j = D - 1; j && i < range0[j]; j--) {
133             ll m = j, coef = 1;
134             while (m /= pri[i]) {
135                 f0[j] += (coef += 3) * get_value_f(i + 1, m);
136             }
137         }
138     }
139     for (int i = 1; i < D; i++) {
140         f1[i] = get_value_f(0, n / i);        
141     }
142     for (int i = 1; i < D; i++) {
143         f0[i] = get_value_f(0, i);
144     }
145 }
146 
147 ll get_value_g(int i, ll j) {
148     if (j >= D) {
149         int rj = n / j;
150         return g1[rj] - max(0, i - range1[rj]);
151     }
152     return g0[j] - max(0, min(i, cnt_P[j]) - range0[j]);
153 }
154 
155 void calculate_g() {
156     for (int i = 1; i < D; i++) {
157         g0[i] = i, g1[i] = n / i;
158     }
159     int cp = 0;
160     for (int i = 0; pri[i] < D; i++, cp++) {
161         for (int j = 1; j < D && i < range1[j]; j++) {
162             g1[j] -= get_value_g(i, n / j / pri[i]);
163         }
164         for (int j = D - 1; j && i < range0[j]; j--) {
165             g0[j] -= get_value_g(i, j / pri[i]);
166         }
167     }
168     for (int i = 1; i < D; i++) {
169         (g1[i] = get_value_g(cp, n / i) - 1) <<= 2;
170     }
171     for (int i = 1; i < D; i++) {
172         (g0[i] = get_value_g(cp, i) - 1) <<= 2;
173     }
174 }
175 
176 void Solve(ll _n) {
177     n = _n;
178     D = sqrt(n + 0.5) + 1;
179     while (D * 1ll * D <= n)
180         D++;
181     precalc();
182     calculate_g();
183     calculate_f();
184     ll ans = f1[1];
185     for (int i = 1; i < D; i++) {
186         ans += d3[i] * g1[i];
187 //        cerr << d3[i] << " " << g1[i] << '\n'; 
188     }
189     printf("%lld\n", ans);
190 }
191 
192 inline void solve() {
193     for (int i = 0; i < (signed) ns.size(); i++) {
194         Solve(ns[i]);
195     }
196 }
197 
198 int main() {
199     init();
200     solve();
201     return 0;
202 }

洲閣篩

 

Min25篩

Part 1

假設小於等於$\sqrt{n}$的素數從小到大分別是$P_1, P_2, \cdots, P_{m}$。

設$g_{i, j} = \sum_{k = 2}^{j} [k的最小質因子大於P_i或者k是質數] F(k)$

轉移式子：

$$g_{i, j} =
\begin{cases}g_{i - 1, j} - F(P_i)(g_{i - 1, \lfloor j / P_i\rfloor} - g_{i - 1, P_i - 1}) & (j \geqslant P_i^2) \\ 
g_{i - 1, j} & (1 \leqslant j < P_i^2)
\end{cases}
$$

和洲閣篩相似，不過第一種狀況還須要摳掉素數的貢獻。

這裏也要求$F$是一個徹底積性函數。

咱們考慮求出$1, 2, \cdots, \sqrt{n}, \lfloor n / \sqrt{n}\rfloor, \cdots, \lfloor n / 2\rfloor, n$處的值。

Part 2

這裏考慮計算$S(n, i) = \sum_{i = 2}^{n}[n的最小質因子大於等於P_i]F(n)$。

那麼顯然答案等於$S(n, 1) + F(1)$。

這裏簡單粗暴一點，不可貴到這樣一個式子：

$$
S(n, i) = g_{m, n} - \sum_{j = 1}^{m - 1}F(P_j) + \sum_{j = i \wedge P_j^2 \leqslant n}\sum_{e = 1\wedge P_{j}^{e + 1}\leqslant n} F(P_{j}^{e}) S(\lfloor n/P_{j}^e\rfloor) + F(P_{j}^{e + 1})
$$

前一半是計算素數的貢獻。

考慮計算合數的貢獻，後一半是先枚舉最小質因子，由於它是一個合數最小質因子，因此$p^2 \leqslant n$

$p^{e + 1} \leqslant n$是相似的緣由。

時間複雜度被證實爲是$O(\frac{n}{Poly(\log n)})$。

（彷佛目前min_25在$10^{13}$範圍內吊打其餘篩）

  1 /**
  2  * loj
  3  * Problem#6053
  4  * Accepted
  5  * Time: 1442ms
  6  * Memory: 3124k
  7  */ 
  8 #include <bits/stdc++.h>
  9 using namespace std;
 10 typedef bool boolean;
 11 
 12 #define ll long long
 13 
 14 const int Mod = 1e9 + 7;
 15 
 16 template <const int Mod = :: Mod>
 17 class Z {
 18     public:
 19         int v;
 20 
 21         Z() : v(0) {    }
 22         Z(int v) : v(v) {    }
 23         Z(ll x) : v(x % Mod) {    }
 24 
 25         Z operator + (Z b) {
 26             int x = v + b.v;
 27             return Z((x >= Mod) ? (x - Mod) : (x));
 28         }
 29         Z operator - (Z b) {
 30             int x = v - b.v;
 31             return Z((x < 0) ? (x + Mod) : (x));
 32         }
 33         Z operator * (Z b) {
 34             return Z(v * 1ll * b.v);
 35         }
 36 };
 37 
 38 typedef Z<> Zi;
 39 
 40 template <typename T>
 41 void pfill(T* pst, const T* ped, T val) {
 42     for ( ; pst != ped; *(pst++) = val);
 43 }
 44 
 45 typedef class Euler {
 46     public:
 47         int *pri;
 48         boolean* vis;
 49 
 50         Euler(int n) {
 51             int num = 0;
 52             pri = new int[(n + 1)];
 53             vis = new boolean[(n + 1)];
 54             pfill(vis + 1, vis + n + 1, false);
 55             for (int i = 2; i <= n; i++) {
 56                 if (!vis[i]) {
 57                     pri[num++] = i;
 58                 }
 59                 for (int j = 0; j < num && pri[j] * 1ll * i <= n; j++) {
 60                     vis[i * pri[j]] = true;
 61                     if (!(i % pri[j])) {
 62                         break;
 63                     }
 64                 }
 65             }
 66         }
 67         ~Euler() {
 68             delete[] vis;
 69         }
 70 } Euler;
 71 
 72 const int N = 1e5 + 30;
 73 
 74 ll n;
 75 int D;
 76 int *pri;
 77 Zi fp[N];
 78 Zi g0[N], g1[N];
 79 Zi h0[N], h1[N];
 80 int range0[N], range1[N];
 81 
 82 inline void init() {
 83     scanf("%lld", &n);
 84     D = sqrt(n + 0.5);
 85     while (D * 1ll * D <= n) {
 86         D++;
 87     }
 88     Euler euler(D + 23);
 89     pri = euler.pri;
 90 }
 91 
 92 int nump = 0;
 93 
 94 Zi get_value_g(ll j) {
 95     return (j >= D) ? (g1[n / j]) : (g0[j]);
 96 }
 97 
 98 Zi get_value_h(ll j) {
 99     return (j >= D) ? (h1[n / j]) : (h0[j]);
100 }
101 
102 Zi S(ll n, int m) {
103     if (pri[m] > n) {
104         return 0;
105     }
106     Zi rt = get_value_g(n) - fp[m];
107     for (int j = m; pri[j] * 1ll * pri[j] <= n; j++) {
108         int p = pri[j];
109         ll pe = p;
110         for (int e = 1; pe * p <= n; pe = pe * p, e++) {
111             rt = rt + Zi(p ^ e) * S(n / pe, j + 1) + Zi(p ^ (e + 1));
112         }
113     }
114 //    assert(m <= nump);
115 //    cerr << n << " " << m << " " << rt.v << '\n';
116     return rt;
117 }
118 
119 Zi comb2(ll n) {
120     return (n & 1) ? (Zi((n + 1) >> 1) * n) : (Zi(n >> 1) * (n + 1));    
121 }
122 
123 inline void solve() {
124     for (int i = 1; i < D; i++) {
125         range0[i] = range0[i - 1];
126         while (pri[range0[i]] * 1ll * pri[range0[i]] <= i) {
127             range0[i]++;
128         }    
129     }
130     range1[D] = range0[D - 1];
131     for (int i = D - 1; i; i--) {
132         range1[i] = range1[i + 1];
133         while (pri[range1[i]] * 1ll * pri[range1[i]] * i <= n) {
134             range1[i]++;
135         }
136     }
137     for (int i = 1; i < D; i++) {
138         g0[i] = i - 1, g1[i] = Zi(n / i - 1);
139         h0[i] = (((i + 1) * 1ll * i) >> 1) - 1;
140         h1[i] = comb2(n / i) - 1;    
141     }
142     for (int i = 0; pri[i] < D; i++, nump++) {
143         for (int j = 1; j < D && i < range1[j]; j++) {
144             g1[j] = g1[j] - (get_value_g(n / j / pri[i]) - g0[pri[i] - 1]);
145             h1[j] = h1[j] - (get_value_h(n / j / pri[i]) - h0[pri[i] - 1]) * pri[i];
146         }
147         for (int j = D - 1; j && i < range0[j]; j--) {
148             g0[j] = g0[j] - (g0[j / pri[i]] - g0[pri[i] - 1]);
149             h0[j] = h0[j] - (h0[j / pri[i]] - h0[pri[i] - 1]) * pri[i];
150         }
151     }
152     for (int i = 1; i < D; i++) {
153         g0[i] = h0[i] - g0[i] + (i >= 2) * 2;
154         g1[i] = h1[i] - g1[i] + (n >= (i << 1)) * 2;
155     }
156     for (int i = 1; i <= nump; i++) {
157         fp[i] = fp[i - 1] + (pri[i - 1] ^ 1);
158     }
159     Zi ans = S(n, 0) + 1;
160     printf("%d\n", ans.v);
161 }
162 
163 int main() {
164     init();
165     solve();
166     return 0;
167 }

min_25篩