# 洲阁筛 & min_25筛学习笔记

### 洲阁筛

$\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)$。

#### 优化

• 维护一个$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$，所以这一部分的值不会改变。

#### 优化

  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

$$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}$$

#### Part 2

$$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^{e + 1} \leqslant 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筛
posted @ 2019-02-23 23:24 阿波罗2003 阅读(...) 评论(...) 编辑 收藏