与素数玩耍

设 $sum(x)$ 表示小于等于 x 的素数个数。

假设我很蠢 ~~（这件事根本不用假设好吗)~~，连 10 以内的素数有哪些都不知道，只知道 1 不是素数。那么我就会令 $sum(x)=x-1$

然后我就会做一些事情来扔掉合数，于是我从小到大处理素数。假设我在处理素数 $p$ ，并且想要扔掉 $sum(i)$ 中以 p 为最小素因子的合数，我就会看向 $sum(\lfloor \frac{i}{p} \rfloor)$ ，由于前期扔掉了一些数，这个里面应该只有素数和素因子大于等于 p 的数，排除掉小于 p 的素数，剩下的数乘以 p 应该就是我们这一轮操作要扔掉的合数，也就是说， $sum(i)-=sum(\lfloor \frac{i}{p} \rfloor)-sum(p-1)$

设 $f1(x)=sum(x)$ , $f2(x)=sum(\lfloor \frac{n}{x} \rfloor)$ 。f1 和 f2 都只用开根号级别的空间。

那么我们就可以灵活运用 f1 和 f2 来完成筛法了。复杂度可以近似拟合为 $O(n^{0.7})$

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1000005;
LL f1[N],f2[N],n,lim;
int main()
{
    scanf("%lld",&n);
    lim=sqrt(n);
    for(LL i=1;i<=lim;++i) f1[i]=i-1,f2[i]=n/i-1;
    for(LL p=2;p<=lim;++p) {
        if(f1[p]==f1[p-1]) continue;//不是素数，则跳过
        //以下都是按照公式 sum[i]-=sum[i/p]-sum[p-1] 进行操作
        for(LL i=1;i<=lim/p;++i) f2[i]-=f2[i*p]-f1[p-1];
        for(LL i=lim/p+1;i<=n/p/p&&i<=lim;++i) f2[i]-=f1[n/i/p]-f1[p-1];
        for(LL i=lim;i>=p*p;--i) f1[i]-=f1[i/p]-f1[p-1];
    }
    printf("%lld\n",f2[1]);
    return 0;
}

你已经学会了与素数玩耍，可以来水一道题了：loj6202 叶氏筛法

其实和上一题差不多，只是公式变成了 $sum(x)-=(sum(\lfloor \frac{i}{p} \rfloor)-sum(p-1)) \times p$

此外，如果你不想写高精的话，最后一个样例是过不去的……

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double db;
const int N=1000005;
const db eps=1e-12;
db f1[N],f2[N];LL l,r;
db sum(LL x) {return (db)(x+1)*(db)x/2.0;}
db query(LL x) {
    if(x==0) return 0;
    LL lim=sqrt(x);
    for(LL i=1;i<=lim;++i) f1[i]=sum(i),f2[i]=sum(x/i);
    for(LL p=2;p<=lim;++p) {
        if(fabs(f1[p]-f1[p-1])<eps) continue;
        db tmp=f1[p-1];
        for(LL i=1;i<=lim/p;++i) f2[i]-=(f2[i*p]-tmp)*(db)p;
        for(LL i=lim/p+1;i<=x/p/p&&i<=lim;++i)
            f2[i]-=(db)(f1[x/i/p]-tmp)*(db)p;
        for(LL i=lim;i>=p*p;--i) f1[i]-=(f1[i/p]-tmp)*(db)p;
    }
    return f2[1];
}
int main()
{
    scanf("%lld%lld",&l,&r);
    printf("%.0lf\n",query(r)-query(l-1));
    return 0;
}

与积性函数玩耍

例题:spoj-DIVCNTK

…… 看代码吧。关键思想是每个 f(x) 的贡献是在处理 x 的最大质因子时算上的，x 的最大质因子次数是 1 的时候可以用前缀和快速处理，否则可以枚举处理。

#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long uLL;
uLL read() {
    uLL q=0;char ch=' ';
    while(ch<'0'||ch>'9') ch=getchar();
    while(ch>='0'&&ch<='9') q=q*10+(uLL)(ch-'0'),ch=getchar();
    return q;
}
const int N=200100;
uLL n,K,lim,tot,T,ans;
uLL f1[N],f2[N],pri[N];
void init(uLL x) {//首先筛出素数个数前缀和
    lim=sqrt(x),tot=0;
    if(x<=1) return;//这句话挺重要的
    for(uLL i=1;i<=lim;++i) f1[i]=i-1,f2[i]=x/i-1;
    for(uLL p=1;p<=lim;++p) {
        if(f1[p]==f1[p-1]) continue;
        pri[++tot]=p;uLL t=f1[p-1];
        for(uLL i=1;i<=lim/p;++i) f2[i]-=f2[i*p]-t;
        for(uLL i=lim/p+1;i<=x/p/p&&i<=lim;++i) f2[i]-=f1[x/i/p]-t;
        for(uLL i=lim;i>=p*p;--i) f1[i]-=f1[i/p]-t;
    }
}
void work(uLL pos,uLL res,uLL rem) {//rem:remain
    uLL t=(rem<=lim?f1[rem]:f2[n/rem]);
    ans+=(K+1)*res*(t-f1[pri[pos-1]]);
    //一举计算多个 f(x) 的贡献，其中 x 的最大质因子次数为 1
    for(uLL i=pos;i<=tot;++i) {
        if(rem<pri[i]*pri[i]) return;
        //若满足这个 if 条件，说明 pri[i] 最多作为一个 x 的最大质因子，且次数只能为 1，没有重复计算的必要
        for(uLL now=pri[i],js=1;now<=rem;now*=pri[i],++js) {
            if(now*pri[i]<rem) work(i+1,res*(K*js+1),rem/now);
            //如果不加前面的 if 判断的话，后面计算的 (t-f1[pri[pos-1]]) 可能会变成负数
            if(js>=2) ans+=res*(K*js+1);//最大质因子的次数不为 1 的 x 对应的 f(x) 在此计算
        }
    }
}
int main()
{
    T=read();
    while(T--) {
        n=read(),K=read();
        ans=1,init(n),work(1,1,n);
        printf("%llu\n",ans);
    }
    return 0;
}

好吧，如果你看懂了上面那道题，那么再来一道？loj6053 简单的函数

这道题和上一道也差不多，就是我们要计算出质数异或 1 的前缀和。

于是我们可以筛出质数前缀和以及质数个数前缀和。由于质数中偶数只有 2，那么只有 2 异或 1 后会加 1，其他都会减 1，这样就能够把质数异或 1 的前缀和求出来了。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod=1000000007,div2=500000004;
const int N=100005;
LL n,lim,tot,ans,f1[N],f2[N],s1[N],s2[N],pri[N];
LL mul(LL a,LL b) {//快速乘，防止炸 long long
    LL re=0;
    for(;b;b>>=1,a=(a+a)%mod) if(b&1) re=(re+a)%mod;
    return re;
}
LL sum(LL x) {return mul(x+1,x)*div2%mod;}
void init() {
    if(n<=1) return;
    lim=sqrt(n);LL tf,ts;
    for(LL i=1;i<=lim;++i)
        f1[i]=sum(i)-1,f2[i]=sum(n/i)-1,s1[i]=i-1,s2[i]=n/i-1;
    for(LL p=2;p<=lim;++p) {
        if(s1[p]==s1[p-1]) continue;
        pri[++tot]=p,tf=f1[p-1],ts=s1[p-1];
        for(LL i=1;i<=lim/p;++i) {
            f2[i]=(f2[i]-(f2[i*p]-tf)*p)%mod;
            s2[i]-=(s2[i*p]-ts)%mod,s2[i]=(s2[i]+mod)%mod;
        }
        for(LL i=lim/p+1;i<=n/p/p&&i<=lim;++i) {
            f2[i]=(f2[i]-(f1[n/i/p]-tf)*p)%mod;
            s2[i]-=(s1[n/i/p]-ts)%mod,s2[i]=(s2[i]+mod)%mod;
        }
        for(LL i=lim;i>=p*p;--i) {
            f1[i]-=(f1[i/p]-tf)*p%mod,f1[i]=(f1[i]+mod)%mod;
            s1[i]-=(s1[i/p]-ts)%mod,s1[i]=(s1[i]+mod)%mod;
        }
    }
    for(LL i=1;i<=lim;++i) {
        f1[i]=(f1[i]-s1[i]+(i>=2LL)*2LL+mod)%mod;
        f2[i]=(f2[i]-s2[i]+(n/i>=2LL)*2LL+mod)%mod;
    }
}
void work(LL pos,LL res,LL rem) {
    LL t=(rem<=lim?f1[rem]:f2[n/rem]);
    ans=(ans+res*(t-f1[pri[pos-1]])%mod)%mod;
    for(LL i=pos;i<=tot;++i) {
        if(rem<pri[i]*pri[i]) return;
        for(LL t=pri[i],js=1;t<=rem;t*=pri[i],++js) {
            if(t*pri[i]<rem) work(i+1,res*(pri[i]^js)%mod,rem/t);
            if(js>=2) ans=(ans+res*(pri[i]^js)%mod)%mod;
        }
    }
}
int main()
{
    scanf("%lld",&n);
    ans=1,init(),work(1,1,n);
    printf("%lld\n",ans);
    return 0;
}

【算法】与扩展埃氏筛玩耍 ——by litble

于2018年3月15日2018年3月15日由litble发布

与素数玩耍

与积性函数玩耍

litble

1 条评论

konnyakuxzy · 2018年3月15日 2:45 下午

发表回复取消回复

【算法】与扩展埃氏筛玩耍 ——by litble

于2018年3月15日2018年3月15日由litble发布

与素数玩耍

与积性函数玩耍

litble

1 条评论

konnyakuxzy · 2018年3月15日 2:45 下午

发表回复 取消回复

发表回复取消回复