Loading

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;
}


posted @ 2024-06-12 19:54  Larunatrecy  阅读(32)  评论(0编辑  收藏  举报