前言
杜教筛学了，顺便把min25筛也学了吧= =刚好多校也有一道题需要补。
下面推荐几篇博客，我之后写一点自己的理解就是了。
传送门1
传送门2
传送门3
这几篇写得都还是挺好的，接下来我就写下自己对min25筛的理解吧 。
正文
简介：
min25筛同杜教筛类似，是用来解决一类积性函数的前缀和，即(sum_{i=1}^nF(i))，并且这里的(n)可以达到(10^{10})的规模。
但所求积性函数要求满足以下条件：
- (F(p))可以表示为简单多项式的形式，比如(p_1^{k_1}+p_2^{k_2}+p_3^{k_3})，(p)为质数；
- (F(p^e))的值要较容易求得。
为什么是这样？往下看就知道啦。
求解：
我们可以将上面说的多项式拆成的每项单独来算前缀和，最后加起来即可。
设(f(p)=p^k)，我们现在就来求(f(i))的前缀和。注意这里(f)函数是完全积性函数，后面我们需要用到这一性质。
Part 1:
首先设(g(n,j))表示：(1,2...n)中，满足条件的(f)之和，满足条件是指要么为质数，要么其最小质因子大于第(j)个质数。
为什么要这么设？后面就知道啦。
之后我们考虑递推求解：
- 若(p_j^2>n)，根据定义，(g(n,j)=g(n,j-1))；
- 否则，我们要减去最小质因子为(p_j)的数，那么就有递推式：(g(n,j)=g(n,j-1)-f(p_j)*(g(frac{n}{p_j},j-1)-sum_{j-1}))。
- 解释一下上式，就相当于首先提取一个(p_j)出来，那么只要减去剩下的质因子大于等于(p_j)的就行啦；但根据定义，(g)中还包含了(sum_{i=1}^{j-1}f(p_i))，我们减去就行了。另外上式还利用了完全积性函数的性质。
这里是min25筛的关键。
我们还可以用筛法的思想去考虑，我们从(p_{j-1})递推到(p_j)的过程，其实就是在埃氏筛中，把(p_j)的倍数筛去。每次得到的(g(n,j))，其实就是利用前(j)个素数来进行埃氏筛，剩下的数以及所有质数的(f)之和。
这也是为啥min25后面有个“筛”的原因~
从筛法的角度来考虑，那么初始化时我们将所有的数都看作素数，然后一个一个来筛求解即可（1除外，我们先不考虑1，最后考虑即可）。
Q:那为什么要求出这个呢？
A:求出(g)之后，我们可以方便地求出(sum_pf(p))的值，即(g(n,|P|))，其中(|P|)为素数集大小。
那合数怎么办？
Part 2:
上面我们把质数的情况考虑了，利用了类似埃氏筛的思想递推得到了(g)函数。
接下来我们考虑合数的情况，因为(F)为积性函数，那么我们直接枚举最小质因子以及其次数来求解即可。
类似地，定义(S(n,j))表示(1,2,cdots,n)中，满足条件的(F)之和，满足条件是指最小质因子大于等于(这里有个等于，其实也可以没有~)第(j)个质数。
那么可以直接暴力求解(S):
稍微解释一下：前面部分就相当于我们求出了满足条件的(F)之和，这里的条件指的是大于等于(p_j)的质数。
接下来我们就相当于暴力枚举最小的质因子以及其次方，并且根据其积性函数的性质，将枚举值提出来，然后递归求解子问题。
但我们对于(p_k^e)这种数没有统计到，所以后面加上就行了。
那么最终答案就为(S(n,1)+F(1))。
所以...min25筛就这么完了。似乎也不是很难。
但我还想加个
Part 3:
虽然min25筛的思想说得差不多了，但我还想说一下其实现的一些细节。毕竟代码也是很重要的。
下面以模板题为例：
Code
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e6 + 5, MOD = 1e9 + 7, inv3 = 333333336;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { // sqrt
chk[1] = 1;
for(int i = 1; i <= n; i++) {
if(!chk[i]) {
prime[++tot] = i;
sum1[tot] = (sum1[tot - 1] + i) % MOD;
sum2[tot] = (sum2[tot - 1] + 1ll * i * i % MOD) % MOD;
}
for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
chk[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
}
void calc_g() {
int z = sqrt(n);
for(ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
w[++cnt] = n / i;
g1[cnt] = w[cnt] % MOD;
g2[cnt] = g1[cnt] * (g1[cnt] + 1) / 2 % MOD * (2 * g1[cnt] + 1) % MOD * inv3 % MOD - 1;
g1[cnt] = g1[cnt] * (g1[cnt] + 1) / 2 % MOD - 1;
if(n / i <= z) ind1[n / i] = cnt;
else ind2[n / (n / i)] = cnt;
}
for(int i = 1; i <= tot; i++) {
for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
ll tmp = w[j] / prime[i], k;
if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
(g1[j] -= prime[i] * (g1[k] - sum1[i - 1] + MOD) % MOD) %= MOD;
(g2[j] -= prime[i] * prime[i] % MOD * (g2[k] - sum2[i - 1] + MOD) % MOD) %= MOD;
if(g1[j] < 0) g1[j] += MOD;
if(g2[j] < 0) g2[j] += MOD;
}
}
}
ll S(ll x, int y) { // 2~x >= P_y
if(x <= 1 || prime[y] > x) return 0;
ll z = sqrt(n);
ll k = x <= z ? ind1[x] : ind2[n / x];
ll ans = (g2[k] - g1[k] + MOD - (sum2[y - 1] - sum1[y - 1]) + MOD) % MOD;
for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
ll pe = prime[i], pe2 = prime[i] * prime[i];
for(int e = 1; pe2 <= x; ++e, pe = pe2, pe2 *= prime[i]) {
ans = (ans + pe % MOD * ((pe - 1) % MOD) % MOD * S(x / pe, i + 1) + pe2 % MOD * ((pe2 - 1) % MOD) % MOD) % MOD;
}
}
return ans % MOD;
}
int main() {
cin >> n;
int tmp = sqrt(n);
pre(tmp);
calc_g();
cout << (S(n, 1) + 1) % MOD ;
return 0;
}
接下来我会说说代码里面的一些要点：
- 预处理部分
没什么好说的，不同的问题预处理也不同。
- 求(g)部分
我先把式子抄下来...懒得翻上去了。
观察这个式子可以发现，第二维每次都是从(j-1)优化过来，所以我们可以直接滚动掉一维。并且，每次递推时，都是从(frac{n}{p_j})转移过来，我们联想到了(lfloorfrac{n}{i}
floor)这个形式。
显然，因为(n)可能会很大，我们不能直接将(n)求出来，质数也是，筛那么大的数，线性筛也会超时。
下面有两个观察：
- 观察1：递推中，有用的素数只有不超过(sqrt{n})的部分；
- 观察2：因为(lfloorfrac{n}{i} floor)只有(sqrt{n})个不同的值，所以我们只需要预处理这(sqrt{n})个值就行了。
因为有了观察1，我们递推到(sqrt{n})就相当于算出了当(j=|P|)的情况。
但可能某些值还是很大，怎么办？
首先我们肯定需要一个数组(w)来储存所有的(lfloorfrac{n}{i}
floor)，之后用了两个(ind)数组来存储下标，如果(lfloorfrac{n}{i}
floor<sqrt{n})，那么直接存储；否则就在另外一个数组存储(lfloorfrac{n}{lfloorfrac{n}{i}
floor}
floor)的下标。
那么之后我们就可以通过(lfloorfrac{n}{i}
floor)的形式来访问相关值了。
这里类似于新加了一个映射关系。(w)数组映射到(lfloorfrac{n}{i}
floor)，而(ind)数组用了点trick把(lfloorfrac{n}{i}
floor)映射回(w)数组。
- 求(S)部分
有一个问题是，这里怎么得到(|P|)。
因为我们把(g)滚动了的，最后虽然求出的是为(sqrt{n})的情况，但稍加思考就会发现，其值等于(|P|)下的值。
怎么得到(i)的下标？注意我们是从(n)开始递推，每次会除以一个数，并且有这样一个形式：(lfloorfrac{lfloorfrac{n}{a}
floor}{b}
floor=lfloorfrac{n}{ab}
floor)，那么每次的(i)都一定可以表示为(lfloorfrac{n}{k}
floor)的形式。所以直接根据刚才的映射关系来找就行啦。这也是为什么需要映射关系的原因
感觉代码部分也说得差不多了，可能有些地方我理解的有问题或者有点复杂，各位看官也不要太过于纠结字眼...自己理解也是不错的办法。
反正我感觉这里利用(lfloorfrac{n}{i}
floor)的性质来求解很关键，优化了代码时间复杂度以及编写难度（求(g)和(S)都用到了它相关性质）。十分巧妙。
时间复杂度：(我也不会)O(感觉能过)(O(frac{n^frac{3}{4}}{log_n}))。
例题：
思路：
因为有(f(p^c)=p^c)，注意到除开(2)以外的素数都为奇数，那么其余(f(p)=p-1)。我们处理的时候将(f(2))的值也当作(1)来算，最后加回来就行。
基本上也是一道模板题吧~
Code
#include <bits/stdc++.h>
#define heyuhhh ok
using namespace std;
typedef long long ll;
const int N = 1e6 + 5, MOD = 1e9 + 7, inv2 = (MOD + 1) / 2;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { // sqrt
chk[1] = 1;
for(int i = 1; i <= n; i++) {
if(!chk[i]) {
prime[++tot] = i;
sum1[tot] = (sum1[tot - 1] + i) % MOD;
}
for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
chk[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
}
void calc_g() {
int z = sqrt(n);
for(ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
w[++cnt] = n / i;
g2[cnt] = (w[cnt] - 1) % MOD;
g1[cnt] = (w[cnt] % MOD) * ((w[cnt] + 1) % MOD) % MOD * inv2 % MOD; --g1[cnt];
if(n / i <= z) ind1[n / i] = cnt;
else ind2[n / (n / i)] = cnt;
}
for(int i = 1; i <= tot; i++) {
for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
ll tmp = w[j] / prime[i], k;
if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
(g1[j] -= prime[i] * (g1[k] - sum1[i - 1] + MOD) % MOD) %= MOD;
(g2[j] -= (g2[k] - i + 1 + MOD) % MOD) %= MOD;
if(g1[j] < 0) g1[j] += MOD;
if(g2[j] < 0) g2[j] += MOD;
}
}
}
ll S(ll x, int y) { // 2~x >= P_y
if(x <= 1 || prime[y] > x) return 0;
ll z = sqrt(n);
ll k = x <= z ? ind1[x] : ind2[n / x];
ll ans = (g1[k] - sum1[y - 1] + MOD - g2[k] + y - 1 + MOD) % MOD;
if(y == 1) ans += 2;
for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
ll pe = prime[i];
for(int e = 1; pe * prime[i] <= x; ++e, pe = pe * prime[i]) {
ans = (ans + (prime[i] ^ e) * S(x / pe, i + 1) % MOD + (prime[i] ^ (e + 1))) % MOD;
}
}
return ans;
}
int main() {
#ifdef heyuhhh
freopen("input.in", "r", stdin);
#else
ios::sync_with_stdio(false); cin.tie(0);
#endif
cin >> n;
int tmp = sqrt(n);
pre(tmp);
calc_g();
cout << (S(n, 1) + 1) % MOD ;
return 0;
}
牛客多校第七场K.Function
题意：
求(sum_{i=1}^nf(i))，其中，(f(i))为：
思路：
考虑(min25)筛求解。
我们不管1，2的存在，最后加上其贡献即可。
首先求出(g_1,g_2)，分别质数表示(\% 4)为1和3的情况总数，然后一个一个来筛。筛的时候要同时考虑当前质数模4的情况以及对(g_1,g_2)的影响，它们都会减去某个值。
然后就直接求和就行了，基本都是套用板子。
最主要的还是求(g)的思路。
Code
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 2e6 + 5;
int T;
ll n;
ll sum1[N], sum2[N], prime[N];
ll w[N], ind1[N], ind2[N];
ll g1[N], g2[N];
bool chk[N];
int tot, cnt;
void pre(int n) { // sqrt
chk[1] = 1;
for(int i = 1; i <= n; i++) {
if(!chk[i]) {
prime[++tot] = i;
sum1[tot] = sum1[tot - 1];
sum2[tot] = sum2[tot - 1];
if(i % 4 == 1) ++sum1[tot];
if(i % 4 == 3) ++sum2[tot];
}
for(int j = 1; j <= tot && prime[j] * i <= n; j++) {
chk[i * prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
}
void calc_g() {
int z = sqrt(n);
for(ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
w[++cnt] = n / i;
g1[cnt] = w[cnt] / 4 + (w[cnt] % 4 >= 1) - 1;
g2[cnt] = w[cnt] / 4 + (w[cnt] % 4 >= 3);
if(n / i <= z) ind1[n / i] = cnt;
else ind2[n / (n / i)] = cnt;
}
for(int i = 1; i <= tot; i++) {
for(int j = 1; j <= cnt && prime[i] * prime[i] <= w[j]; j++) {
ll tmp = w[j] / prime[i], k;
if(tmp <= z) k = ind1[tmp]; else k = ind2[n / tmp];
if(prime[i] % 4 == 1) {
g1[j] -= (g1[k] - sum1[i - 1]);
g2[j] -= (g2[k] - sum2[i - 1]);
} else if(prime[i] % 4 == 3) {
g1[j] -= (g2[k] - sum2[i - 1]);
g2[j] -= (g1[k] - sum1[i - 1]);
}
}
}
}
int f(int p, int q) {
if(p % 4 == 1) return 3 * q + 1;
return 1;
}
ll S(ll x, int y) { // 2~x >= P_y
if(x <= 1 || prime[y] > x) return 0;
ll z = sqrt(n);
ll k = x <= z ? ind1[x] : ind2[n / x];
ll ans = 4 * g1[k] - 4 * sum1[y - 1] + g2[k] - sum2[y - 1];
if(y == 1) ans++;
for(int i = y; i <= tot && prime[i] * prime[i] <= x ; i++) {
ll pe = prime[i], pe2 = prime[i] * prime[i];
for(int e = 1; pe2 <= x; ++e, pe = pe2, pe2 *= prime[i]) {
ans += f(prime[i], e) * S(x / pe, i + 1) + f(prime[i], e + 1);
}
}
return ans;
}
int main() {
freopen("input.in", "r", stdin);
cin >> T;
while(T--) {
memset(chk, 0, sizeof(chk)); tot = cnt = 0;
cin >> n;
int tmp = sqrt(n);
pre(tmp);
calc_g();
cout << S(n, 1) + 1 << '
';
}
return 0;
}
后记
写得还是比较匆忙，感觉也没有很好地把心静下来，可能有些地方写得冗余、累赘或者有错误，请谅解。
总得来说，min25筛的思想以及代码实现部分都是很巧妙的~为什么会想到搞个(g)来啊？为什么这样复杂度就比较低啊？这些都是问题...
另外加一些我学习时的草稿：
因为转移的原因，所以有用的素数只有小于等于(sqrt{n})的。最后算出来的值就是(g(n,|P|))啦。
处理(g)时，离散化(lfloorfrac{n}{i} floor)，因为第一维为(n)不会变，并且有个这样性质：(lfloorfrac{lfloorfrac{n}{a} floor}{b} floor=lfloorfrac{n}{ab} floor)，所以有用的值也就只有(lfloorfrac{n}{i} floor)这些。
离散化要点：利用(w)递减记录离散化的值，并且利用(ind1,ind2)来记录相关数的下标，不超过(sqrt{n})。涉及整除的性质。巧妙！