神仙 Min_25 发明的神奇筛子，用来筛积性函数前缀和。

对于要筛前缀和的积性函数 \(f(x)\)，其具体要求是：\(f(p)\) 是一个简单多项式，\(f(p^e)\) 可以快速计算。

复杂度是亚线性的，但我不会证。

Description

对于一个积性函数 \(f(x)\) 的前缀和，我们可以做如下变换：

\[\begin{aligned} \sum_{i=1}^nf(x)=&\sum_{p\le n}f(p)+\sum_{\substack{i\le n\\i\not\in\mathbb P}}f(i)\\ =&\sum_{p\le n}f(p)+\sum_{\substack{p\le \sqrt n\\ p^e\le n}}\sum_{\substack{i\le\lfloor\frac n{p^e}\rfloor\\ \operatorname{lpf}(i)>p}}f(i) \end{aligned} \]

第一行就是把 \(1\sim n\) 的数的贡献分为素数和合数分别进行计算。

第二行就是枚举每个合数的最小质因子，再计算剩下的，其中 \(\operatorname{lpf}(i)\) 代表 \(i\) 的最小质因子。

然后就可以分为两个部分来计算了！

Part 1 前半部分

已知 \(f(p)\) 是关于 \(p\) 的简单多项式，那么我们可以把这个多项式拆成若干个单项式来进行计算。也就是说，我们只需求出 \(\sum_{p\le n}p^k\) 即可，其中 \(p^k\) 是原多项式的一部分。

考虑进行一个神奇的 dp，定义 dp 数组 \(g(n,j)\) 为：

\[g(n,j)=\sum_{i\le n}[i\in \mathbb P\operatorname{or}\operatorname{lpf}(i)>p_j]i^k \]

如果把算法过程想象成埃氏筛的话，\(g(n,j)\) 计算的其实就是筛去 \(p_j\) 及倍数后剩余的所有“素数”对答案的贡献。

考虑从 \(g(n,j-1)\) 转移过来，也就是说我们要计算筛掉 \(p_j\) 对答案的贡献。把 \(p_j\) 分成两类讨论：

• \(p_j>\sqrt n,p_j^2>n\)：可知最小质因数为 \(p_j\) 的最小合数为 \(p_j^2\)，因此筛去 \(p_j\) 对答案没有贡献，可得 \(g(n,j)=g(n,j-1)\)。
• \(p_j\le\sqrt n,p_j^2\le n\)：枚举除去最小质因子 \(p_j\) 后被筛掉的合数的值，可知筛去 \(p_j\) 对答案的贡献为 \(p_j^kg(\lfloor\frac n{p_j}\rfloor,j-1)\)。但是这样我们会多筛一次 \(p_1\sim p_{j-1}\)，因此答案要再加上 \(p_j^kg(p_j-1,j-1)\)。

于是我们得到了 \(g(n,j)\) 的递推式：

\[g(n,j)=g(n,j-1)-p_j^k \begin{cases} 0&p_j>\sqrt n,p_j^2>n\\ g(\lfloor\frac n{p_j}\rfloor,j-1)-g(p_j-1,j-1)&p_j\le\sqrt n,p_j^2\le n \end{cases} \]

其中，\(g(p_j-1,j-1)\) 其实就是 \(\sum_{i=1}^{j-1}p_j^k\)，可以线性筛出，把它记为 \(\textit{Sp}_j\)。

注意到当 \(p_j>\sqrt n\) 时，\(g(n,j)\) 的值都是相同的，于是这样的 \(g(n,j)\) 的值可以不用筛。同时，对于每个 \(x\)，我们需要的只有所有 \(g(\lfloor\frac nd\rfloor,x)\) 的值，而这样的值只有 \(O(\sqrt n)\) 个，存储时映射一下下标即可。

例题：P5493 质数前缀统计

设 \(S(n)\) 为 \(n\) 以内质数的 \(k\) 次方和。
给定 \(N\)，求：

\[\sum_{i=1}^{\lfloor\sqrt N\rfloor}i^2S(\lfloor\frac Ni\rfloor)\bmod p \]

\(0 \le k \le 10\)，\(1 \le N \le 4\times {10}^{10}\)，\({10}^9 < p < 1.01 \times {10}^9\) 。

板子题，按照上面的过程实现即可。计算自然数 \(k\) 次方和可参考 CF622F。

code（用了 FastMod 优化取模）：

#include <cstdio>
#include <cmath>

const int maxn = 1e6 + 5;
typedef long long ll;

namespace FastModnmsp {
    typedef unsigned int uint;
    typedef __uint128_t uL;
    typedef unsigned long long ull;

    struct FastMod {
        ull b,m;
        FastMod() {}
        FastMod(ull _b): b(_b),m(ull((uL(1) << 64) / _b)) {}
        inline void init(ull _b) { b = _b,m = ull((uL(1) << 64) / _b); }
        friend inline ull operator % (const ull& a,const FastMod& mod) {
            ull r = a - (uL(mod.m) * a >> 64) * mod.b;
            return r >= mod.b ? r - mod.b : r;
        }
    } Mod;
} using FastModnmsp::Mod;

ll n,sq,k,ans,P,finv[15],prsk[15];
ll cntP,pri[maxn],prik[maxn],sumk[maxn];
ll cntS,w[maxn],g[maxn];
int id1[maxn],id2[maxn];
bool isp[maxn];

inline ll read() {
#define gc c = getchar()
    ll d = 0; int f = 0,gc;
    for(;c < 48 || c > 57;gc) f |= (c == '-');
    for(;c > 47 && c < 58;gc) d = (d << 1) + (d << 3) + (c ^ 48);
#undef gc
    return f ? -d : d;
}

inline ll fpow(ll a,ll b) {
    ll res = 1;
    for(;b;a = a * a % Mod,b >>= 1) if(b & 1) res = res * a % Mod;
    return res;
}

inline ll getPowSum(ll n) {
    ll ans = 0,p = 1,pre = 0;
    for(int i = 1;i <= k + 2;i ++) {
        p = p * (n - i) % Mod,pre = (pre + prsk[i]) % Mod;
        if(n == i) return pre;
        ll inv = finv[i - 1] * finv[k + 2 - i] % Mod * fpow(n - i,P - 2) % Mod;
        ll res = pre * inv % Mod;
        if((k - i) & 1) res = P - res;
        ans = (ans + res) % Mod;
    }
    return p * ans % Mod;
}

inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void init() {
    finv[0] = isp[1] = 1; sq = sqrt(n);
    for(ll i = 1;i <= k + 2;i ++) finv[i] = finv[i - 1] * fpow(i,P - 2) % Mod,prsk[i] = fpow(i,k);
    for(ll i = 2;i <= sq;i ++) {
        if(!isp[i]) pri[++ cntP] = i,
                    prik[cntP] = fpow(i,k),
                    sumk[cntP] = (sumk[cntP - 1] + prik[cntP]) % Mod;
        for(int j = 1;j <= cntP && i * pri[j] <= sq;j ++) {
            isp[i * pri[j]] = true;
            if(!(i % pri[j])) break;
        }
    }
    for(ll r,l = 1;l <= n;l = r + 1) {
        r = n / (n / l); w[++ cntS] = n / l;
        g[cntS] = getPowSum(w[cntS] % Mod) - 1;
        if(n / l <= sq) id1[n / l] = cntS;
        else id2[r] = cntS;
    }
    for(int i = 1;i <= cntP;i ++)
        for(int j = 1;j <= cntS && pri[i] * pri[i] <= w[j];j ++) {
            ll h = getId(w[j] / pri[i]);
            g[j] -= prik[i] * (g[h] - sumk[i - 1] + P) % Mod;
            if(g[j] < 0) g[j] += P;
        }
}

int main() {
    n = read(),k = read(),P = read(); Mod.init(P);
    init();
    for(ll i = 1;i <= sq;i ++) ans = (ans + i * i % Mod * g[i] % Mod) % Mod;
    printf("%lld\n",ans);
    return 0;
}

Part 2 求解答案

还是 dp。考虑借鉴前半部分，定义 dp 数组 \(S(n,j)\)：

\[S(n,j)=\sum_{i\le n}[\operatorname{lpf}(i)>p_j]f(i) \]

答案显然就是 \(S(n,0)+1\)。这里的 \(1\) 是 \(f(1)\) 的贡献，把它提出来计算。

接下来分两部分讨论贡献：

• 素数：所有素数的贡献为 \(g(n,j)\)，但是要减掉 \(p_1\sim p_{j-1}\) 的贡献，因此答案是 \(g(n,j)-\textit{Sp}_j\)。
• 合数：枚举最小质因数及其次数，可得贡献为 \(\sum_{\substack{k>j\\p_k^e\le n}}f(p_k^e)(S(\lfloor\frac n{p_k^e},k)+[e\not=1])\)，其中 \([e\not=1]\) 计算的是除去 \(p_k^e\) 后 \(f(1)\) 的贡献。

所以最终的式子就是：

\[S(n,j)=g(n,x)-\textit{Sp}_j+\sum_{\substack{k>j\\p_k^e\le n}}f(p_k^e)(S(\lfloor\frac n{p_k^e},k)+[e\not=1]) \]

其中 \(x\) 为 \(n\) 以内使得 \(p_x^2\le n\) 的最大的 \(x\)，即 \(x=\pi(\lfloor\sqrt n\rfloor)\)。

很神奇的是，这玩意直接爆搜就行了，根本不需要记忆化，而且复杂度还是优秀的 \(O(\frac{n^{0.75}}{\log n})\)。

例题：

1. P5325 【模板】Min_25筛

令积性函数 \(f(p^k)=p^k(p^k-1)\)。给定 \(n\)，求：

\[\sum_{i=1}^nf(i)\bmod (10^9+7) \]

\(1\le n\le 10^{10}\)。

在质数处，显然有 \(f(p)=p(p-1)=p^2-p\)，可以 Min_25 筛。

code：

#include <cstdio>
#include <cmath>

typedef long long ll;
const int maxn = 1e6 + 5;
const ll P = 1e9 + 7,inv3 = (P + 1) / 3;

ll n,sq,tot,cnt;
ll pri[maxn]; bool isp[maxn];
ll Sp1[maxn],Sp2[maxn];
ll g1[maxn],g2[maxn],w[maxn];
int id1[maxn],id2[maxn];

inline ll getLiSum(ll n) { return n * (n + 1) / 2 % P; }
inline ll getSqSum(ll n) { return n * (n + 1) / 2 % P * (2 * n + 1) % P * inv3 % P; }

inline ll Add(ll a,ll b) { a += b; return a >= P ? a - P : a; }

inline int getId(ll k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void Init() {
    isp[1] = true,sq = sqrt(n);
    for(ll i = 2;i <= sq;i ++) {
        if(!isp[i]) pri[++ cnt] = i,
                    Sp1[cnt] = Add(Sp1[cnt - 1],i),
                    Sp2[cnt] = Add(Sp2[cnt - 1],i * i % P);
        for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
            isp[i * pri[j]] = true;
            if(!(i % pri[j])) break;
        }
    }
    for(ll r,l = 1;l <= n;l = r + 1) {
        r = n / (n / l);
        w[++ tot] = n / l;
        g1[tot] = getLiSum(w[tot] % P) - 1;
        g2[tot] = getSqSum(w[tot] % P) - 1;
        if(n / l <= sq) id1[n / l] = tot;
        else id2[r] = tot;
    }
    for(int i = 1;i <= cnt;i ++)
        for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
            int k = getId(w[j] / pri[i]);
            g1[j] -= pri[i] * (g1[k] - Sp1[i - 1] + P) % P;
            g2[j] -= pri[i] * pri[i] % P * (g2[k] - Sp2[i - 1] + P) % P;
            if(g1[j] <= 0) g1[j] += P; if(g2[j] <= 0) g2[j] += P;
        }
}

inline ll S(ll n,int j) {
    if(pri[j] >= n) return 0;
    ll x = getId(n),ans = Add(((g2[x] - Sp2[j]) - (g1[x] - Sp1[j])) % P,P);
    for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
        ll p = pri[k];
        for(int e = 1;p <= n;e ++,p *= pri[k]) {
            ll tmp = p % P;
            ans = Add(ans,tmp * (tmp - 1) % P * (S(n / p,k) + (e != 1)) % P);
        }
    }
    return ans;
}

int main() {
    scanf("%lld",&n); Init();
    printf("%lld\n",Add(S(n,0),1));
    return 0;
}

2. SP34096 DIVCNTK

令 \(\sigma_0(n)\) 为 \(n\) 的约数个数。给定 \(n,k\)，求：

\[\sum_{i=1}^n\sigma_0(i^k)\bmod 2^{64} \]

多测，\(1\le n,k\le10^{10}\)。

三倍经验题，写完这题就可以切了 DIVCNT2 和 DIVCNT3。

令 \(f(n)=\sigma_0(n^k)\)，则显然有：

\[f(n)=\begin{cases} 1&n=1\\ ke+1&n=p^k\\ f(a)f(b)&n=ab,a\perp b \end{cases} \]

可以 Min_25 筛出来。

code：

#include <cstring>
#include <cstdio>
#include <cmath>

typedef unsigned long long ull;
const int maxn = 1e6 + 5;

ull n,K,sq;
ull pri[maxn]; bool isp[maxn];
ull Sp[maxn],g[maxn],w[maxn];
int t,tot,cnt,id1[maxn],id2[maxn];

inline int getId(ull k) { return k <= sq ? id1[k] : id2[n / k]; }

inline void Init() {
    isp[1] = true,sq = sqrt(n);
    for(ull i = 2;i <= sq;i ++) {
        if(!isp[i]) pri[++ cnt] = i,
                    Sp[cnt] = Sp[cnt - 1] + K + 1;
        for(int j = 1;j <= cnt && i * pri[j] <= sq;j ++) {
            isp[i * pri[j]] = true;
            if(!(i % pri[j])) break;
        }
    }
    for(ull r,l = 1;l <= n;l = r + 1) {
        r = n / (n / l);
        w[++ tot] = n / l;
        g[tot] = w[tot] - 1;
        if(n / l <= sq) id1[n / l] = tot;
        else id2[r] = tot;
    }
    for(int i = 1;i <= cnt;i ++)
        for(int j = 1;j <= tot && pri[i] <= w[j] / pri[i];j ++) {
            int k = getId(w[j] / pri[i]);
            g[j] -= g[k] - i + 1;
        }
    for(int i = 1;i <= tot;i ++) g[i] *= (K + 1);
}

inline ull S(ull n,int j) {
    if(pri[j] >= n) return 0;
    ull x = getId(n),ans = g[x] - Sp[j];
    for(int k = j + 1;k <= cnt && pri[k] <= n / pri[k];k ++) {
        ull p = pri[k];
        for(ull e = 1;p <= n;e ++,p *= pri[k])
            ans += (K * e + 1) * (S(n / p,k) + (e != 1));
    }
    return ans;
}

int main() {
    scanf("%d",&t);
    for(;t;t --) {
        memset(isp,false,sizeof(isp)); tot = cnt = 0;
        scanf("%llu%llu",&n,&K); Init();
        printf("%llu\n",S(n,0) + 1);
    }
    return 0;
}