ZigZagK的博客
Min_25筛
2019年4月9日 18:18
查看标签

一类问题

已知积性函数 $f(n)$ ,其中 $f(p)$ 是简单多项式,且 $f(p^k)$ 可以快速计算( $p$ 是素数),求其前缀和。

上杜教筛?如果 $f(n)​$ 很奇怪就没法卷另外一个函数上去了,所以要用到Min_25筛。

预处理

定义 $h(n)​$ 表示把 $n​$ 强行当成素数代入 $f(n)​$ 中的值,比如 $f(p)=p​$ 那么 $h(n)=n​$ ,且 $h(n)​$ 需要是完全积性函数

定义 $g_{i,j}$ 表示 $[1,i]$ 中是素数或者最小质因数 $>p_j$ 的 $h(n)$ 的和,初值 $g_{i,0}=\sum_{i=2}^{i}h(i)$ ,考虑如何递推:

  • $p_j^2>i​$ :$g_{i,j}=g_{i,j-1}​$ ,因为最小质因数 $>p_j​$ ,所以 $g_{i,j}=g_{i,j-1}​$ 。
  • $p_j^2\le i$ :$g_{i,j}=g_{i,j-1}-h(p_j)[g_{\lfloor{i\over p_j}\rfloor,j-1}-\sum_{k=1}^{j-1}h(p_k)]$ ,$g_{i,j-1}$ 比 $g_{i,j}$ 多出来一部分即最小质因数为 $p_j$ 的部分,又因为 $h(n)$ 是完全积性函数,所以可以先把 $h(p_j)$ 提出来,乘上最小质因数 $\ge p_j$ 的部分。

这样我们就可以得到 $\sum_{i=1}^{n}[i\in P]f(i)​$ 也就是 $g_{n,|P|}​$ ,实际处理中只需要开一维数组,并且只需要 $O(\sqrt n)​$ 个。可以证明时间复杂度 $O({n^{3\over 4}\over log_2n})​$ 。证明?我当然不会。

所以可以用这个方法求出素数个数(即 $h(n)=1​$ ),HDU5901

#include<cmath>
#include<cstdio>
using namespace std;
typedef long long LL;
const int maxs=632456;

LL L,R,n;int p[maxs+5];bool pri[maxs+5];
int S,m,ID[2][maxs+5];LL a[maxs+5],g[maxs+5];

inline void Make(int n){
    for (int i=2;i<=n;i++){
        if (!pri[i]) p[++p[0]]=i;
        for (int j=1,t;j<=p[0]&&(t=i*p[j])<=n;j++)
            {pri[t]=true;if (!(i%p[j])) break;}
    }
}
inline LL Min25(LL x){
    n=x;S=sqrt(n);m=0;
    for (LL l=1,r;l<=n;l=r+1){
        r=n/(n/l);a[++m]=n/l;g[m]=a[m]-1;
        a[m]<=S?ID[0][a[m]]=m:ID[1][n/a[m]]=m;
    }
    for (int j=1;j<=p[0]&&p[j]<=S;j++)
        for (int i=1;i<=m&&(LL)p[j]*p[j]<=a[i];i++)
            g[i]+=j-1-g[a[i]/p[j]<=S?ID[0][a[i]/p[j]]:ID[1][n/(a[i]/p[j])]];
    return g[1];
}
int main(){
    freopen("program.in","r",stdin);freopen("program.out","w",stdout);
    for (Make(maxs);~scanf("%lld",&n);printf("%lld\n",Min25(n)));return 0;
}

求前缀和

定义 $S(n,j)$ 表示 $[1,n]$ 中最小质因数 $\ge p_j$ 的 $f(n)$ 的和,那么 $\sum_{i=1}^{n}f(i)$ 就是 $S(n,1)+f(1)$ 。

由于我们已经知道了素数的答案,只需要加上合数,枚举合数的最小质因子以及次数,则:

$$ S(n,j)=g_n-\sum_{i=1}^{j-1}f(p_i)+\sum_{i=j}\sum_{k=1}S(\lfloor{n\over p_i^k}\rfloor,i+1)f(p_i^k)+f(p_i^{k+1}) $$

不用记忆化,可以证明时间复杂度 $O({n^{3\over 4}\over log_2n})$ 。证明?我当然不会。

比如算莫比乌斯函数前缀和,由于莫比乌斯函数 $\mu(p)=-1$ 对应的 $h(n)=-1$ 不是完全积性函数,所以不能直接预处理,但是我们可以通过预处理 $h(n)=1$ 来求出 $S(n,j)$ ,51Nod1244

#include<cmath>
#include<cstdio>
using namespace std;
typedef long long LL;
const int maxs=200000;

LL L,R,n;int p[maxs+5];bool pri[maxs+5];
int S,m,ID[2][maxs+5];LL a[maxs+5],g[maxs+5];

inline void Make(int n){
    for (int i=2;i<=n;i++){
        if (!pri[i]) p[++p[0]]=i;
        for (int j=1,t;j<=p[0]&&(t=i*p[j])<=n;j++)
            {pri[t]=true;if (!(i%p[j])) break;}
    }
}
LL Sum(LL a,int b){
    if (a<1||p[b]>a) return 0;LL ans=-g[a<=S?ID[0][a]:ID[1][n/a]]+b-1;
    for (int i=b;i<=p[0]&&(LL)p[i]*p[i]<=a;i++) ans-=Sum(a/p[i],i+1);return ans;
}
inline LL Min25(LL x){
    n=x;S=sqrt(n);m=0;
    for (LL l=1,r;l<=n;l=r+1){
        r=n/(n/l);a[++m]=n/l;g[m]=a[m]-1;
        a[m]<=S?ID[0][a[m]]=m:ID[1][n/a[m]]=m;
    }
    for (int j=1;j<=p[0]&&p[j]<=S;j++)
        for (int i=1;i<=m&&(LL)p[j]*p[j]<=a[i];i++)
            g[i]+=j-1-g[a[i]/p[j]<=S?ID[0][a[i]/p[j]]:ID[1][n/(a[i]/p[j])]];
    return Sum(n,1)+1;
}
int main(){
    freopen("program.in","r",stdin);freopen("program.out","w",stdout);
    scanf("%lld%lld",&L,&R);Make(maxs);printf("%lld\n",Min25(R)-Min25(L-1));return 0;
}

跑得比杜教筛快多了,可能是因为这个常数小,而杜教筛要记忆化常数比较大的原因……

版权声明:本博客所有文章除特别声明外,均采用 CC BY 4.0 CN协议 许可协议。转载请注明出处!