min_25筛学习笔记

作用

对积性函数 (F(n)) 求前缀和。

需要满足的条件：

1. 对于 (p) 是素数，(F(p)) 是多项式。
2. 对于 (p) 是素数，(F(p^k)) 可以快速计算（由k从小到大递推计算也可以）。

复杂度据说是 (O(frac{n^{3/4}}{logn})) ，在 (n le 10^{12}) 基本上都能快速出解。空间复杂度是 (O(sqrt n))

原理

分为两个阶段：预处理、计算答案。思路是先求所有素数的函数和，再求合数的，最后考虑 (F(1)) 。

预处理

由于 (F(p)) 是多项式，所以可以拆分成若干个 (p^k) 的形式。记 (f(p)=p^k) （这里的 p 可以不是素数），下面考虑求 (sum_{i=1}^n f(i)[ ext{i is prime}])

(g(n,j)) 表示小于等于 n 中所有质数的 (f) 和最小素因子大于 (P_j) 的合数的 (f) 之和。易得

[g(n,j)= egin{cases} g(n,j-1)&P_j^2gt n\ g(n,j-1)-P_{j}^k[g(lfloorfrac{n}{P_j} floor,j-1)-g(P_j-1,j-1)]&P_j^2le n end{cases} ]

对于下面那个式子，减号后面利用了 (f) 是完全积性函数的性质。即把 (P_j^k) 提出来，剩下的数字的最小质因数也要大于 (P_j) ，于是中括号内是最小质因子大于 (P_j) 的所有数的贡献。实际上也可以写成

[g(n,j)=g(n,j-1)-P_{j}^k[g(lfloorfrac{n}{P_j} floor,j-1)-sum_{i=1}^{P_j-1}f(i)[ ext{i is prime}]] ]

由于 (lfloorfrac{n}{i} floor) 只会取根号级别个，所以只需要对这些数计算 (g) 即可。而第二维，由于 (P_j^2>n) 的所有 j 没有贡献，所以只会取到 (O(sqrt n / ln n)) 。于是以 j 为阶段就可以进行 dp。

最终要求的 (sum_{i=1}^n f(i)[ ext{i is prime}]) 就是 (g(n,|P|)) 。

初始化：(g(n,0) = sum_{i=2}^nf(i))

计算答案

与预处理中类似符号的定义不同，这里定义 (S(n,j)) 是小于等于 n 中最小素因子大于等于 (P_j) 的数的 (F) 之和。则

[S(n,j)=G(n,|P|)-sum_{i=1}^{j-1}F(P_i)+sum_{kge j}sum_{e=1}^{P_k^{e+1}le n}(F(P_k^e)S(frac{n}{P_k^e},k+1)+F(P_k^{e+1})) ]

其中 (G(n,|P|)) 为 (sum_{i=1}^n F(i)[ ext{i is prime}]) 。这其实是合并了预处理中若干个 (f(p)=p^k) 的结果。

这个式子即：大于等于 (P_j) 的素数的 (F) 之和，加上最小质因子大于等于 (P_j) 的合数的 (F) 之和。

于是最终答案为 (S(n,1) + F(1)) 。

应用

【模板】Min_25筛

#include<bits/stdc++.h>
using namespace std;

typedef long long LL;
typedef long double LD;
typedef pair<int,int> pii;
typedef pair<LL,int> pli;
const int SZ = 1e7 + 10;
const int INF = 1e9 + 10;
const int mod = 1e9 + 7;
const LD eps = 1e-8;

LL read() {
    LL n = 0;
    char a = getchar();
    bool flag = 0;
    while(a > '9' || a < '0') { if(a == '-') flag = 1; a = getchar(); }
    while(a <= '9' && a >= '0') { n = n * 10 + a - '0',a = getchar(); }
	if(flag) n = -n;
	return n;
}
// 此题：f(p^k)=p^k(p^k-1)
struct Min_25 {

    LL n;
    int gen; // sqrt(n)
    LL w[SZ]; // 所有的 n/i 下取整
    int id1[SZ],id2[SZ]; // w中对应数字的标号
    int sump[SZ],sump2[SZ]; // 前 i 个素数的和、平方和
    int g1[SZ],g2[SZ]; // 计算完毕后为小于等于 w[i] 的所有质数的对应函数和

    bool vis[SZ];
    int pri[SZ],tot;
    void shai(int n) {
        tot = 0;
        for(int i = 2;i <= n;i ++) {
            if(!vis[i]) pri[++ tot] = i;
            for(int j = 1,m;j <= tot && (m = i * pri[j]) <= n;j ++) {
                vis[m] = 1;
                if(i%pri[j] == 0) break;
            }
        }
        for(int i = 1;i <= tot;i ++) {
            sump[i] = (sump[i-1] + pri[i]) % mod; /// modify
            sump2[i] = (sump2[i-1] + 1ll*pri[i]*pri[i]%mod) % mod; /// modify
        }
    }

    int getid(LL x) {
        if(x <= gen) return id1[x];
        return id2[n/x];
    }

    int S(LL x,int y) {
        if(x <= 1 || pri[y] > x) return 0;
        int k = getid(x);
        int ans = (g2[k]-g1[k]-(sump2[y-1]-sump[y-1])) % mod; /// modify
        for(int i = y;i <= tot && 1ll * pri[i] * pri[i] <= x;i ++) {
            LL p = pri[i];
            LL t1 = p,t2 = p*p;
            for(int e = 1;t2 <= x;e ++) {
                (ans += (S(x/t1,i+1) * (t1 % mod)%mod * ((t1-1) % mod) % mod /// modify
                         + (t2%mod) * ((t2-1)%mod) % mod)%mod) %=mod;
                t1 *= p; t2 *= p;
            }
        }
        return ans;
    }

    int f1(LL x) {
        x %= mod;
        return x * (x+1) / 2 % mod;
    }

    int f2(LL x) {
        x %= mod;
        int ni6 = 166666668;
        return x * (x+1) % mod * (2*x+1) % mod * ni6 % mod;
    }

    int work(LL nn) {
        n = nn;
        gen = sqrt(n);
        shai(gen);
        int m = 0;
        for(LL i = 1,r;i <= n;i = r + 1) {
            r = n/(n/i); w[++m] = n/i;
            g1[m] = f1(w[m])-1; /// modify
            g2[m] = f2(w[m])-1; /// modify
            if(w[m] <= gen) id1[w[m]]=m;
            else id2[r]=m;
        }
        for(int j = 1;j <= tot;j ++) {
            LL p = pri[j];
            for(int i = 1;i<=m && p*p<=w[i];i ++) {
                int k = getid(w[i]/p);
                (g1[i] -= p * (g1[k] - sump[j-1]) % mod) %= mod; /// modify
                (g2[i] -= p * p % mod * (g2[k] - sump2[j-1]) % mod) %= mod; /// modify
            }
        }
        int ans = S(n,1) + 1;
        ans += mod; ans %= mod;
        return ans;
    }
}min25;

int main() {
    LL n = read();
    printf("%d
",min25.work(n));
}

参考资料

min_25筛