Min_25 筛

用来求积性函数 \(f(x)\) 的前缀和。要求其在质数 \(p\) 处的取值为多项式，并且 \(f(p^k)\) 可以快速计算。

因为多项式可以拆成单项式，所以只用求出形如 \(f(p)=p^k\) 的前缀和，然后加起来就行。

设 \(g(n,j)\) 为小于等于 \(n\) 的所有质数和最小质因子大于第 \(j\) 个质数的合数的 \(k\) 次方和，即：

\[\large g(n,j)=\sum_{i=1}^n [i \in Prime \or \text{minp}(i) > p_j] i^k \]

考虑递推，得：

\[\large\begin{aligned} g(n,j)&=g(n,j-1)-\sum_{i=1}^n[i \not\in Prime \and \text{minp}(i)=p_j]i^k \\ &=g(n,j-1)-p_j^k\sum_{i=1}^{\frac{n}{p_j}}[\text{minp}(i) \geqslant p_j]i^k \\ &=g(n,j-1)-p_j^k\left( g(\frac{n}{p_j},j-1)-g(p_{j-1},j-1) \right)[p_j^2\leqslant n] \end{aligned} \]

\(g(n,j-1)\) 比 \(g(n,j)\) 多的部分为小于等于 \(n\) 的最小质因子为 \(p_j\) 的合数的 \(k\) 次方和，提出 \(p_j\) 后，剩下的部分要求最小质因子大于等于 \(p_j\)，即 \(g(\frac{n}{p_j},j-1)\)，但这里多算了质数，所以还要减去 \(g(p_{j-1},j-1)\)。考虑到只有 \(p_j \leqslant \sqrt n\) 时，后一项才有值，因此设 \(sp(n)=\sum_{i=1}^n p_i^k\)，得 \(g(p_{j-1},j-1)=sp(j-1)\)，这一部分可以线性筛预处理。

一开始处理出 \(g(n,0)\)，然后递推计算到 \(g(n,x)\)，其中 \(x=\text{maxp}(n)\)，也就是 \(g(n,x)\) 为小于等于 \(n\) 的所有质数的 \(k\) 次方和。

因为 \(n\) 过大，无法求出所有 \(g(n,x)\)，但是发现形如 \(\left\lfloor \frac{n}{a} \right\rfloor\) 的取值只有 \(O(\sqrt n)\) 个，因此只用求这 \(O(\sqrt n)\) 个位置的 \(g(n,x)\) 即可。

处理出 \(g(n,x)\) 的复杂度为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)。

求前缀和时，类似的设 \(s(n,j)\) 为小于等于 \(n\) 的所有最小质因子大于 \(p_j\) 的数处的函数值的和，即：

\[\large s(n,j)=\sum_{i=1}^n [\text{minp}(i)>p_j]f(i) \]

将其分为两部分来考虑，分为大于 \(p_j\) 的质数和最小质因子大于 \(p_j\) 的合数，质数部分即为 \(g(n,x)-sp(j)\)，合数部分枚举最小质因子即可，得：

\[\large s(n,j)=g(n,x)-sp(j)+\sum_{k>j \and p_k^e\leqslant n}f(p_k^e)\left( s(\frac{n}{p_k^e},k)+[e \ne 1] \right) \]

因为当 \(e \ne 1\) 时没有计算 \(f(p_k^e)\) 的值，所以还要加上 \([e \ne 1]\)。实现时就直接递归计算，不用记忆化。

递归的复杂度为 \(O(n^{1-\epsilon })\)。

求积性函数 \(f(x)\) 的前缀和，其满足 \(f(p^k)=p^k(p^k-1)\)：

#include<bits/stdc++.h>
#define maxn 200010
#define p 1000000007
#define inv2 500000004
#define inv6 166666668
using namespace std;
typedef long long ll;
template<typename T> inline void read(T &x)
{
    x=0;char c=getchar();bool flag=false;
    while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
    while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
    if(flag)x=-x;
}
ll n,t,tot,cnt;
ll pri[maxn],g1[maxn],g2[maxn],sp1[maxn],sp2[maxn],id1[maxn],id2[maxn],q[maxn];
bool tag[maxn];
int id(ll x)
{
    return x<=t?id1[x]:id2[n/x];
}
void init()
{
    for(int i=2;i<=t;++i)
    {
        if(!tag[i]) pri[++tot]=i;
        for(int j=1;j<=tot;++j)
        {
            int k=i*pri[j];
            if(k>t) break;
            tag[k]=true;
            if(i%pri[j]==0) break;
        }
    }
}
ll s(ll n,int j)
{
    if(pri[j]>=n) return 0;
    int x=id(n);
    ll v=((g2[x]-sp2[j]-(g1[x]-sp1[j]))+p)%p;
    for(int k=j+1;k<=tot&&pri[k]*pri[k]<=n;++k)
    {
        ll now=pri[k];
        for(int e=1;now<=n;++e,now*=pri[k])
        {
            ll val=now%p;
            v=(v+val*(val-1+p)%p*(s(n/now,k)+(e!=1)))%p;
        }
    }
    return v;
}
int main()
{
    read(n),t=sqrt(n),init();
    for(int i=1;i<=tot;++i)
        sp1[i]=(sp1[i-1]+pri[i])%p,sp2[i]=(sp2[i-1]+pri[i]*pri[i])%p;
    for(ll l=1,r;l<=n;l=r+1)
    {
        ll v=n/l;
        r=n/v,q[++cnt]=v,v%=p;
        g1[cnt]=((v+1)*v%p*inv2-1+p)%p;
        g2[cnt]=(v*(v+1)%p*(2*v+1)%p*inv6-1+p)%p;
        if(q[cnt]<=t) id1[q[cnt]]=cnt;
        else id2[n/q[cnt]]=cnt;
    }
    for(int j=1;j<=tot;++j)
    {
        for(int i=1;i<=cnt&&pri[j]*pri[j]<=q[i];++i)
        {
            int x=id(q[i]/pri[j]);
            g1[i]=(g1[i]-pri[j]*(g1[x]-sp1[j-1]+p)+p)%p;
            g2[i]=(g2[i]-pri[j]*pri[j]%p*(g2[x]-sp2[j-1]+p)+p)%p;
        }
    }
    printf("%lld",(s(n,0)+1)%p);
    return 0;
}