min_25筛

问题引入

求${\sum }_{i=1}^{n}F(i)$，且需要满足：

• $F(i)$是一个积性函数且$F(p)$为关于$p$低阶多项式，例如$F(p)=p(p-1)$
• $F(p^k)$可以快速算出

min_25筛 = 数论+DP+离散化

min_25筛

因为多项式可以拆分成多个单项式，所以我们只需考虑求出$f(p)=p^k$的前缀和，注意到$f(i)$是一个完全积性函数，在下面的运算中有着重要作用。还有一点就是如下的讨论不包含$f(1)=1$，只需在最后加上$1$即为答案。

一些记号

$P_i$，第$i$个素数

$lpf(n)$：$n$的最小质因子

分类

首先对每一个$i$按质数和合数分类：

${\sum }_{i=1}^{n}f(i)={\sum }_{p\in primes,p\le n}f(p)+{\sum }_{j\notin primes,j\le n}f(j)$

然后枚举合数的小于$n$的最小质因子的幂次，不难知道范围内合数的最小质因子都一定小于等于$\sqrt{n}$：

${\sum }_{i=1}^{n}f(i)={\sum }_{p\in primes}f(p)+{\sum }_{p\le \sqrt{n}}{\sum }_{e=1,p^e\le n}f(p^e){\sum }_{j\le \frac{n}{p^e}\&\&lpf(j)>p}f(j)$

整个式子被分为两个部分，第一个部分是所有质数的和，第二部分为枚举范围内最小质因子幂次后，所有$j$的最小质因子大于这个质因子的$f(j)$之和。

第一部分

考虑一个递推公式$g(j,n)$表示$n$以内筛过前$j$个质数的倍数的合数的所有数之和，即$g(j,n)={\sum }_{i\in priems||lpf(i)>P_j}f(i)={\sum }_{i\in priems||lpf(i)>P_j}{i}^k$

实际上就是埃氏筛的筛合数部分并减去这些合数，如下代码的$sum$的值即为$f(p)=p$下的$g(j,n)$：

//设f(p)=p
int sum = n * (n + 1) / 2;
for (int i = 1; i <= j; i++) {
	for (int k = prime[j] * prime[j]; k <= n; k += prime[j]) {
		sum -= k;
	}
}

考虑如何状态转移，显然$g(j,n)$比$g(j-1,n)$筛的合数更多，那么其和会变小，因此应该考虑减去哪些才能正确转移，这些数应该是最小质因子恰好为$P_j$的合数，用$\lfloor \frac{n}{P_j}\rfloor$得到的是$n$以内有多少$P_j$的倍数，那就减去$g(j,\lfloor \frac{n}{P_j}\rfloor )$。但是有一部分的最小质因子并不是$P_j$而是它前面的更小的，即$f(P_j){\sum }_{i=1}^{j-1}\ast f(P_i)$这部分多减了，那么就要再加上。因此递推式就得到了：

$g(j,n)=g(j-1,n)-f(P_j)[g(j-1,\lfloor \frac{n}{P_j}\rfloor )-{\sum }_{i=1}^{j-1}\ast f(P_i)]=g(j-1,n)-P_j^k[g(j-1,\lfloor \frac{n}{P_j}\rfloor )-{\sum }_{i=1}^{j-1}\ast P_i^k]$

根据上述埃氏筛的代码不难发现当$P_j^2>n$时，$g(j,n)=g(j-1,n)$，因此我们只需计算出$\sqrt{n}$以内的素数。

但是这个二维的$DP$由于范围太大还是无法写，因为$g(j-1,\lfloor \frac{n}{P_j}\rfloor )$较难处理。这时引入一个数论中常见且重要的性质：$\lfloor \frac{n}{\frac{a}{b}}\rfloor =\lfloor \frac{n}{ab}\rfloor$

联系到整除分块，也就是说对于给定的$n$，$\lfloor \frac{n}{x}\rfloor$这样形式的结果最多只有$\sqrt{n}$个，但是又来了一个问题，这$\sqrt{n}$个结果如何访问，这时需要离散化，但是我们还需要知道离散化后的下标访问，这时可以使用$unordered\_map$，但还是会比较慢。我们考虑使用两个数组$id1[maxn],id2[maxn]$建立两个映射代替$map$，前者代表$x$这个数的下标，后者代表$\frac{n}{x}$这个数的下标，这样两个数组最多都只会到$\sqrt{n}$。

对于$g(j-1,n),g(j-1,\lfloor \frac{n}{P_j}\rfloor )$，我们只需在状态转移时采用滚动数组得到，设$g(i)$为当前考虑的素数下从大到小第$i$个$\lfloor \frac{n}{x}\rfloor$对应的值（因为$w[tot]$是非递增的），那么从小到大枚举每个质数，然后根据如上推导的状态转移。

PS：可以对整除分块过程中的$n/l,n/(n/l)$分别输出，可以发现前者是从$n$开始递减的，后者是递增的，且极值会在$\sqrt{n}$处取到。那么在$\sqrt{n}$前我们将$n/(n/l)$作为下标保存到$id2$，之后的将$n/l$作为下标保存到$id1$。

第二部分

我们已经求出了所有质数的和，接下来的步骤仍然是$DP$的思想，设$S(n,j)$表示$n$以内所有的最小质因子大于等于$P_j$的$f(x)$的和，即$S(n,j)={\sum }_{i=1}^{n}F(i),lpf(i)\ge P_j$。按照上述分类之后，我们能得到关于$S$的如下转移方程：

S ( n , j ) = ∑ { g t ( i d [ n ] ) − s u m t ( j − 1 ) } + ∑ k ≥ y , P k 2 ≤ n ∑ e = 1 , P k e + 1 ≤ n F ( P k e ) ∗ S ( n P k e , k + 1 ) + F ( P k e + 1 ) S(n,j)=\sum\{g_t(id[n])-sum_t(j-1)\}+\sum_{k \geq y,P_k^2 \leq n}\sum_{e=1,P_k^{e+1}\leq n}F(P_k^e)*S(\frac{n}{P_k^e},k+1)+F(P_k^{e+1}) S(n,j)=∑{gt​(id[n])−sumt​(j−1)}+∑k≥y,Pk2​≤n​∑e=1,Pke+1​≤n​F(Pke​)∗S(Pke​n​,k+1)+F(Pke+1​)

其中$su{m}_t(j-1)={\sum }_{i=1}^{j-1}{P}_i^t$为前$j-1$个质数的$t$次幂之和， ∑ { g t ( i d [ n ] ) − s u m t ( j − 1 ) } \sum\{g_t(id[n])-sum_t(j-1)\} ∑{gt​(id[n])−sumt​(j−1)}为多个单项式加减运算凑成多项式的质数和${\sum }_{p\in primes}F(x),p<P_j$；后面的是直接计算的多项式函数的$F(P_k^e),F(P_k^{e+1})$。需要注意这个递归代码并不需要记忆化，否则会$TLE$。

最终的答案即$S(n,1)+1$

模板代码

设$F(p^k)=p^k(p^k-1)$，那么求${\sum }_{i=1}^{n}F(i)$的代码为：

const ll Mod = 1e9 + 7, inv2 = 5e8 + 4, inv3 = 333333336;
const int maxn = 1e6 + 10;

ll prime[maxn], sum1[maxn], sum2[maxn];  //p的前缀和以及p^2的前缀和
ll w[maxn], id1[maxn], id2[maxn], g1[maxn], g2[maxn];
ll n, cnt;
int m, tot;  //m=sqrt(n)
bitset<maxn> vis;

void init() {
    vis.reset();
    m = sqrt(n + 0.5), tot = cnt = 0;
    //线性筛
    for (int i = 2; i <= m; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            sum1[cnt] = (sum1[cnt - 1] + i) % Mod;
            sum2[cnt] = (sum2[cnt - 1] + 1LL * i * i) % Mod;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= m; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }

    //预处理每种n/l的前缀和
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;  //在w中保存的数不能取模
        g1[tot] = w[tot] % Mod;
        g2[tot] = g1[tot] * (g1[tot] + 1) / 2 % Mod * (g1[tot] * 2 + 1) % Mod * inv3 % Mod - 1;  //x=n/l,g2(tot)=x(x+1)(2x+1)/6-1
        g1[tot] = g1[tot] * (g1[tot] + 1) / 2 % Mod - 1;  //x=n/l,g1(tot)=x(x+1)/2-1
        if (n / l <= m) id1[n / l] = tot;
        else id2[r] = tot;
    }

    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= tot && prime[i] * prime[i] <= w[j]; j++) {
            ll k = (w[j] / prime[i]) <= m ? id1[w[j] / prime[i]] : id2[n / (w[j] / prime[i])];
            g1[j] -= prime[i] * (g1[k] - sum1[i - 1] + Mod) % Mod;
            g2[j] -= prime[i] * prime[i] % Mod * (g2[k] - sum2[i - 1] + Mod) % Mod;
            g1[j] %= Mod, g2[j] %= Mod;
            if (g1[j] < 0) g1[j] += Mod;
            if (g2[j] < 0) g2[j] += Mod;
        }
    }
}

ll S(ll x, int y) {
    if (x <= 1 || prime[y] > x) return 0;
    ll k = x <= m ? id1[x] : id2[n / x];
    ll ans = (g2[k] - g1[k] + Mod - (sum2[y - 1] - sum1[y - 1]) + Mod) % Mod;
    for (int i = y; i <= cnt && prime[i] * prime[i] <= x; i++) {
        ll pe = prime[i];
        for (int e = 1; pe * prime[i] <= x; e++, pe *= prime[i]) {
            ll xx = pe % Mod, res = pe * prime[i] % Mod;
            ans = (ans + xx * (xx - 1) % Mod * S(x / pe, i + 1) % Mod + res * (res - 1) % Mod) % Mod;
        }
    }
    return ans % Mod;
}

ll solve() {
    init();
    return (S(n, 1) + 1) % Mod;
}

应用一：求$[1,n]$内素数的个数

LibreOJ-6235

此时$f(x)=[x\in primes]$已经不是积性函数了，但是我们仍可以利用$min\_25$筛的第一部分对素数个数进行筛选。设$g(j,n)$为要么为质数要么最小质因子大于$P_j$的数的个数。根据$min\_25$筛的思想，考虑$g(j-1,n)$应该减去多少才能得到$g(j,n)$，不难得知要减去$P_j$的倍数的合数$f(P_j)g(j-1,\lfloor \frac{n}{P_j}\rfloor )$，然后因为会多减了前面几个质数倍数的$P_i\ast P_j(i<j)$，那么就加上$f(P_j){\sum }_{i=1}^{j-1}f(P_i)$，得到如下的状态转移方程：

$g(j,n)=g(j-1,n)-f(P_j)[g(j-1,\lfloor \frac{n}{P_j}\rfloor )-{\sum }_{i=1}^{j-1}\ast f(P_i)]=g(j-1,n)-g(j-1,\lfloor \frac{n}{P_j}\rfloor )-(j-1)$

const ll Mod = 1e9 + 7, inv2 = 5e8 + 4, inv3 = 333333336;
const int maxn = 1e6 + 10;

ll prime[maxn];
ll w[maxn], id1[maxn], id2[maxn], g[maxn];
ll n, cnt;
int m, tot;  //m=sqrt(n)
bitset<maxn> vis;

void init() {
    vis.reset();
    m = sqrt(n + 0.5), tot = cnt = 0;
    for (int i = 2; i <= m; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= m; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }

    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;
        g[tot] = w[tot] - 1;
        if (n / l <= m) id1[n / l] = tot;
        else id2[r] = tot;
    }
}

ll solve() {
    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= tot && prime[i] * prime[i] <= w[j]; j++) {
            ll k = (w[j] / prime[i]) <= m ? id1[w[j] / prime[i]] : id2[n / (w[j] / prime[i])];
            g[j] -= g[k] - i + 1;
        }
    }
    return g[1];
}

应用二：求$[1,n]$内素数之和

2020CCPC网络预选赛

在$min\_25$筛的第一部分，我们设$g(j,n)$表示的是被前$j$个素数筛掉$n$以内其倍数的合数后的和，那么显然我们恰好需要的是$[1,n]$范围内所有合数被筛过的和，即$g(x,n)$，$x$是$\sqrt{n}$范围内的最大素数，考虑到滚动数组进行$DP$后$n$以内所有的合数都已经被筛去，$g[i]$表示第$i$大$\lfloor \frac{n}{x}\rfloor$范围内素数的和，显然$i=1$时$w[1]=n$，即$n$以内的素数之和（可以打表验证）

PS：这道题目卡的有点紧，因此在二重循环处将取模都去掉才能过

ll prime[maxn], sum[maxn];
ll w[maxn], id1[maxn], id2[maxn], g[maxn];
bitset<maxn> vis;
int m, cnt, tot;
ll n, Mod = 1e9 + 7;

void init() {
    vis.reset(), cnt = tot = 0;
    vis[1] = 1;
    m = sqrt(n + 0.5);
    for (int i = 2; i <= m; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            sum[cnt] = sum[cnt - 1] + i;
        }
        for (int j = 1; j <= cnt && i * prime[j] <= m; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }

    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;
        g[tot] = w[tot] % Mod;
        g[tot] = g[tot] * (g[tot] + 1) / 2 - 1;
        if (n / l <= m) id1[n / l] = tot;
        else id2[r] = tot;
    }

    for (int i = 1; i <= cnt; i++) {
        for (int j = 1; j <= tot && prime[i] * prime[i] <= w[j]; j++) {
            ll k = (w[j] / prime[i]) <= m ? id1[w[j] / prime[i]] : id2[n / (w[j] / prime[i])];
            g[j] -= (g[k] - sum[i - 1] + Mod) % Mod * prime[i] % Mod;
            g[j] %= Mod;
            if (g[j] < 0) g[j] += Mod;
        }
    }
}

ll solve() {
    init();
    return g[1];
}