Min25 筛法

之前学习的的确是太浅薄了，于是重新学习一下。

• 可以做什么？

对于满足条件的积性函数 $f(n)$，求其前 $n$ 项和 $\sum_{i=1}^nf(i)$。

• 需要准备些什么？

设 $N=\{x|x=\lfloor\frac{n}{i}\rfloor\}$，即为所有整除的不同点值。

设 $B=\sqrt n,p_{k}$ 代表 $\le B$ 的所有质数。

设 $mn(n)$ 为 $n$ 的最小质因子。

第一部分#

我们需要对于 $i\in N$，求出 $\sum_{p_k\leq i}f(p_k)$。

这里要求：$f(p_k)=\sum_{j=0}c_j(p_k)^j$，即可以表示成一个关于 $p_k$ 的低次多项式，那么我们把多项式拆成若干个单项式分别求和即可，记为 $g(n)$。

设 $G_k(n)$ 为用 $p_{1}\to p_{k}$ 做完线性筛后剩下的数的 $g(i)$ 之和，也就是 $G_k(n)=\sum_{j=1}^k g(p_j)+\sum_{i=p_k+1}^n[mn(i)>p_k]g(i)$。

我们让 $k$ 从 $0$ 依次增大，这样枚举完 $\leq B$ 的质数后就是要求的了。

$G_0(n)$ 就是一个自然数幂和，不多赘述，现在我们看怎么从 $G_{k-1}$ 推到 $G_k$。

相比于 $G_{k-1}(n)$，$G_k(n)$ 筛掉了所有 $mn(i)=p_k$ 的 $i$，我们从 $i$ 中提出一个 $p_k$，然后一定有 $mn(\frac{i}{p})> p_{k-1}$，这个东西就是 $G_{k-1}(n)$ 去掉前边质数的贡献，因此我们有：

$$G_k(n)=G_{k-1}(n)-g(n)(G_{k-1}(\lfloor \frac{n}{p_k}\rfloor )-sum_{k-1}))$$

其中 $sum_k$ 是前 $k$ 个质数的贡献。

这样就完成了递推。

第二部分#

设 $S_k(n)$ 为 $\sum_{i=1}^n[mn(i)>p_k]f(i)$，我们要求出 $S_0(n)$。

$S_k(n)$ 的贡献分成质数和合数两部分，质数就是 $G_0(n)-sum_{k}$。

对于合数，我们枚举最小质因子 $p_{j}(j>k)$ 及其次数 $c$，显然这里 $p_j\leq B$。

那么贡献就是

$$S_k(n)=G_0(n)-sum_k+\sum_{j>k}\sum_cf(p_j^c)(S_{j,\lfloor\frac{n}{p_j^c}\rfloor}+[c\neq 1])$$

直接递归即可。

两者的复杂度都可以证明是 $O(\frac{n^{\frac{3}{4}}}{\ln\ln n})$，实际上比杜教筛还要快。

拓展#

我们能不能对于 $i\in N$，求出 $S_k(i)$ 呢？

我们知道杜教筛可以记忆化，但是这里即使记忆化复杂度依旧爆炸，我们考虑把递归改成递推的。

首先改一下 $S_k(n)$ 的定义：我们只考虑合数的贡献不考虑质数，这样做的好处是我们定义域满足 $p_k^2\leq n$ 。

那么式子就是：

$$S_k(n)=\sum_{j>k}\sum_cf(p_j^c)(S_{j,\lfloor\frac{n}{p_j^c}\rfloor}+G_0(\lfloor\frac{n}{p_j^c}\rfloor)-sum_j+[c\neq 1])$$

$$S_k(n)=S_{k+1}(n)+\sum_cf(p_{k+1}^c)(S_{k+1,\lfloor\frac{n}{p_{k+1}^c}\rfloor}+G_0(\lfloor\frac{n}{p_{k+1}^c}\rfloor)-sum_{k+1}+[c\neq 1])$$

显然和递归是一个复杂度的，那么我们就保持复杂度不变的同时完成了递推所有值。

下面是一个 $f(n)=2^{\omega(n)}$ 的代码，$\omega(n)$ 是 $n$ 的质因子个数。

#pragma GCC optimize(2)
#include<bits/stdc++.h>
using namespace std;
template <typename T>inline void read(T &x)
{
    x=0;char c=getchar();bool f=0;
    for(;c<'0'||c>'9';c=getchar())f|=(c=='-');
    for(;c>='0'&&c<='9';c=getchar())x=(x<<1)+(x<<3)+(c-'0');
    x=(f?-x:x);
}
const int N = 1e6+7;
#define int long long
typedef long long LL;
int sum[N];
namespace sub
{
    int v[N],g[N],prime[N],tot=0;
    int solve(int n)
    {
        g[1]=1;
        for(int i=2;i<=n;i++)
        {
            if(!v[i])
            {
                v[i]=i;
                prime[++tot]=i;
                g[i]=2;
            }
            for(int j=1;j<=tot;j++)
            {
                if(prime[j]>v[i]||1ll*i*prime[j]>n)break;
                int k=i*prime[j];
                v[k]=prime[j];
                if(i%prime[j]==0)
                {
                    g[i*prime[j]]=g[i];
                    break;
                }
                g[i*prime[j]]=g[i]*2;
            }
        }
        for(int i=1;i<=n;i++)sum[i]=sum[i-1]+g[i];
    }
}
int v[N],g[N];
int prime[N],tot=0;
int sg[N],R;
void init(int n)
{
    g[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!v[i])
        {
            v[i]=i;
            prime[++tot]=i;
            g[i]=2;
            sg[tot]=sg[tot-1]+g[i];
        }
        for(int j=1;j<=tot;j++)
        {
            if(prime[j]>v[i]||1ll*i*prime[j]>n)break;
            int k=i*prime[j];
            v[k]=prime[j];
            if(i%prime[j]==0)
            {
                g[i*prime[j]]=g[i];
                break;
            }
            g[i*prime[j]]=g[i]*2;

        }
    }
}
int num[N],cnt=0,G[N],id1[N],id2[N],B;
int F[N],S[N];
void init()
{
    for(int j=tot;j>=0;j--)
    {
        for(int u=1;u<=cnt&&prime[j]*prime[j]<=num[u];u++)
        {
            int x=num[u];
            F[u]=S[u];
        }
        if(j==0)break;
        for(int u=1;u<=cnt&&prime[j]*prime[j]<=num[u];u++)
        {
            int x=num[u];
            int c=1,cur=prime[j];
            for(;cur<=x;cur=cur*prime[j],c++)
            {
                int v=(x/cur<=B?id1[x/cur]:id2[R/(x/cur)]);
                int w=(x/cur>prime[j]?F[v]+G[v]-sg[j]:0);
                S[u]=S[u]+2*(w+(c!=1));
            }
        }
    }
}
signed main()
{
    freopen("a.in","r",stdin);
    //freopen("count.out","w",stdout);
    int n;
    read(n);R=n;
    B=sqrt(n);
    init(B);
    for(int l=1,r;l<=n;l=r+1)
    {
        r=(n/(n/l));
        num[++cnt]=(n/l);
        G[cnt]=(n/l)-1;
        if(n/l<=B)id1[n/l]=cnt;
        else id2[n/(n/l)]=cnt;
    }
    for(int j=1;j<=tot;j++)
    for(int i=1;i<=cnt&&1ll*prime[j]*prime[j]<=num[i];i++)
    {
        int x=num[i]/prime[j];
        if(x<=B)x=id1[x];
        else x=id2[n/x];
        G[i]-=(G[x]-(j-1));
    }
    for(int i=1;i<=cnt;i++)G[i]=G[i]*2;
    init();
    cout<<F[1]+G[1]+1;
    return 0;
}