洲阁筛学习总结

3 篇文章 0 订阅

洲阁筛学习总结

By Semiwaker

问题描述

给出一个积性函数 F(x) ,满足积性函数的基本性质:
F(1)=1
如果 a 和 b 互质,则 F(ab)=F(a)F(b)
x=pkii ,则有

F(x)=F(pkii)

特别的
p为质数,F(p) 是一个低阶多项式, F(pc) 可以快速算出。
现在要求
i=1nF(i)

n的范围大概在 1011 左右。


小范围问题的解决方法

接下来的部分是线性筛,可以跳过
假如n在 107 范围内,考虑怎么线性地求出每一个 F(x)。
分析线性筛的过程,大概分为几个要点:
1、质数的时候需要直接计算。
2、对于某个非质数 x,会分解为两部分:最小的质因数 p 和剩余的部分 y 。有 py=x

质数的部分就不必多说了,只能直接算。
非质数的部分分为3类。
1、p和y互质,那么显然 F(x)=F(p)F(y)
2、y是p的幂 y=pc ,也就是 y 中只有 p 一个质因子,此时需要知道次数 c 然后直接计算。(特殊情况下可以从 F(y) 推过来)
3、y是p的倍数 p|y ,但是y中还有其他质因子,设 y=pcz ,则 F(x)=F(z)F(pc+1) 。如果我们知道 pc ,那么我们就可以这么计算: F(x)=F(ypc)F(pc×p)
我们分析需要知道什么:第二类需要知道每个数最小的质因数的次数 Cnt[x] ,第三类需要知道每个数中最小的质因子构成的部分 Fir[x]。
幸运的是,这两样都可以随着线性筛的过程得出。

具体来说,线性筛的过程如下:
1、由之前筛的结果判断 i 是否质数,质数的情况下,直接算出F(i),Fir[i]=i,Cnt[i]=1。非质数的数据已经在之前算出了。
2、枚举已经求出的质数 P[j] 。
i×P[j] 设为非质数。
3、判断 P[j] 是否是 i 的因数。
4、如果不是, F(i×P[j])=F(i)F(P[j])
Fir[i×P[j]]=P[j] Cnt[i×P[j]]=1
5、如果是,分判断 i 是否等于 Fir[i] ,当 i=Fir[i] 时,i中只有P[j]一个质因数,此时 F(i×P[j])=F(P[j]Cnt[i]+1) ,按照定义直接计算。
否则 F(i×P[j])=F(iFir[i])F(Fir[i]×P[j]) ,显然两项都小于i,已经计算过了。
然后无论如何 Fir[i×P[j]]=Fir[i]×P[j] Cnt[i×P[j]]=Cnt[i]+1


前置技能

引理1

如果 x>a,b
xab=xab
证明:
y=xa
k=xab ,则有 x=kab+c ,显然 c<ab
那么 y=kb+ca ,因为 c<ab ,所以 ca<b ,所以 yb=k

引理2

nx 最多有 2n 种取值。
证明
xn 时,只有 n 种取值。
x>n 时, nxn ,也只有 n 种取值。

这个上限是比较精准的,大多数情况下是 2n1 或者刚好 2n 个。

定理:从n开始整除任意数任意次得到的不同结果只有 2n 种。

根据引理2,n整除1次的结果只有 2n 种。
根据引理1,n整除1次的结果包含整除2次的结果,以此类推,整除任意次的结果都属于n整除1次的结果。

洲阁筛核心思想

因为 n=1011 ,所以我们需要在低于线性时间的复杂度内计算。
考虑线性筛的瓶颈在哪里?为什么线性筛必须是O(n)的?
因为线性筛一定要把1~n的每一个质数都单独处理

n范围内的质数个数大概是 O(nlnn) 范围的,由于ln n只有几十,所以枚举所有质数需要的时间至少O(n)级别。

现在的关键点是:怎样才能不用枚举所有的质数?
核心思想:把质数分类
质数可以简单地分为两类:第一类是小于等于 n 的质数,另外一类是大于 n 的。

为什么要这样分类呢?
因为n以内的数,最多只有一个大于 n 的质因数。
那么,如果我们把小于等于 n 的质数筛掉,剩下的就全部都是质数。
所以,我们可以通过分类来减少对质数的枚举


算法

总述

我们把1~n的每一个数分类,有大于 n 的质数的,和没有的。
那么可以写成这样:

i=1nF(i)=i=1nF(i)(p>nniF(p)[p])+i=1nF(i)[in]

也就是说,有大于 n 的质因子的,我们枚举 n 的质因子的部分,再枚举 >n 的质因子。
没有的就直接枚举。

那么我们的算法现在可以分为三个部分。
第一个部分是枚举 n 以内的每一个i,求F(i),这部分我们可以线性筛。
第二个部分是对于每一个 ni ,求

p>nniF(p)[p]

第三个部分是求
i=1nF(i)[in]

第二部分

p>nniF(p)[p]

我们只需要筛掉 n 内的所有质数,那么剩下的就是我们要求的。
考虑怎样简化问题,因为 F(p) 是一个低阶的多项式,所以我们只需要对于每一项单独求和,第k项求出 pk ,然后合并即可。
即,现在的问题是,对于某个范围[1,m]内的所有大于 n 的质数,求 pk

筛的过程可用递推实现。
设小于等于 n 的质数有Pn个,从小到大排列为 P1..Pn
设在 [1,j] 中,与 P1..i 互质的所有数的k次方和为:g[i][j]。

显然边界为g[0][j]=j。

我们要求的为g[pn][j]。

考虑怎样从前 i-1 个质数推到前 i 个质数。
我们只需要算出与前i-1个质数互质的数,再减去与前i-1个质数互质,但是不与第i个质数互质的数即可。
设某个与前 i-1 个质数互质的数为x,那么 P[i]x 必然与第i个质数不互质。
(P[i]x)k=P[i]kx
P[i]xj ,那么有 xjP[i]
所以得到:

g[i][j]=g[i1][j]Pkig[i1][jP[i]]

根据之前提到的定理,j的取值有 2n 种,i的取值有 O(nlnn) 种,乘起来还是O(n)级别。需要优化。

注意到,如果 Pi>j ,那么除了1以外的所有数都被筛去了。
显然只要j不为0,那么g[i][j]=1。

怎么利用这个性质呢?
如果 Pi1>jP[i] ,那么 g[i1][jP[i]]=1 .
放宽一点 Pi>jP[i] ,此时的判断条件转化为 P2i>j
换句话说,
P2i>j 时,有 g[i][j]=g[i1][j]Pki
那么设 P2L[j]>jP2L[j]1 ,即L[j]是第一个满足条件的i。
对于 i<L[j] ,我们需要老老实实地去算。对于 iL[j] ,我们只需要一个前缀和就可以了。

还有一点要注意的:
如果 Pi>j ,此时 jP[i]=0 ,那么有

g[i][j]=g[i1][j]

设小于等于j的质数有CntP[j]个,那么对于 iCntP[j] ,是不需要去计算的。

综合起来,根据i的大小,可以分为3段。
第一段0到L[j]-1,直接算。
第二段L[j]到CntP[j],预处理 pk 的前缀和。
第三段CntP[j]到Pn,不需要算。

观察到递推式只和i-1有关,是可以滚动数组的,所以具体的实现方法如下:
i<L[j] 时,直接 g[j]=g[j]Pkig[jP[i]]
iL[j] 时,停止对 g[j] 的递推。
如果在某个i调用到了g[j],那么我们再算上缺少的L[j]~min(i,CntP[j])的部分的前缀和即可。
设前缀和为 SumP[i]=it=0pkt ,那么我们可以这样表示 g[i][j]:
g[j]+max(0,SumP[min(i,CntP[j])]SumP[L[j]])

时间复杂度

由于每一个j对应的需要枚举的i是 jlnj 范围的,所以我们可以列式:

i=1nnilnniO(n34lnn)

细节

1、CntP[j]只需要保存 jn 的。因为当 j>n 时,所有i都满足 Pij 。此时只要当CntP[j]=Pn即可。
2、关于j如何储存的问题。肯定不是暴力哈希。
由于j的范围是[1,n],所有直接开数组是不可行的。
我们可以分类存放。
jn ,我们直接开数组,存放在g0[j]。
j>n ,设 nx=j ,存放在g1[x]。

3、由于使用了滚动数组,所以 j 要从大到小枚举。而 L[j] 是随着 j 减小而减小的,而我们枚举的条件是 i<L[j] ,所以在 iL[j] 直接停止枚举即可。

第三部分

i=1nF(i)[in]

同样我们用递推去求。
这次为了方便,我们从大的质数往小的质数推。

没有大于 n 的质因子,那么就意味着只有 n 以内的质因子。
所以我们设:
x为在[1,j]的范围内,只有 Pi..Pn 这些质因子的数。(即i以后的质因子。)
f[i][j]=F(x)

边界为 f[Pn+1][j]=1

要求的为f[0][j]。

那么有:

f[i][j]=f[i+1][j]+F(Pci)f[i+1][jPci]

简单来说,直接枚举当前质因子的幂,然后考虑乘起来不会超过j的,只包含 Pi+1..Pn 的数的和为多少。因为是积性函数,所以可直接乘起来。

同样,直接算这个递推式也是O(n)级别的,需要优化。

同样考虑 P2i>j ,则 jPi<Pi ,此时 f[i+1][jPi]=1 (因为只包含大于Pi的质数)。
而当 c2 jPci=0
所以有

f[i][j]=f[i+1][j]+F(Pi)

同样当 Pi>j 时, f[i][j]=f[i+1][j]

所以,同样根据i的大小分为3段。(注意是从大到小的。)
第一段 Pn 到 CntP[i]+1,这一段不用计算。
第二段 CntP[i] 到 L[i],这一段我们要计算F(p)的前缀和。
第三段 L[i]-1 到 0,这一段老老实实递推。

同样地,可以用滚动数组。

设前缀和为 SumP[i]=it=1F(t)
那么,在某个 i 时,f[j]的实际值为 f[j]+Max(0,SumP[L[i]-1]-SumP[Max(i,CntP[i])-1])

时间复杂度

同样为

O(n34lnn)

细节

同样,当 jn 时,记为 f0[j],否则设 j=nx ,记为f1[j]

计算答案

现在回到原本的式子:

i=1nF(i)=i=1nF(i)(p>nniF(p)[p])+i=1nF(i)[in]

实际需要计算的是:
Ans=ni=1F(i)g(ni)+f(n)

Ans=ni=1F(i)g1(i)+f1(1)

注意到 g0 和 f0 都没有用了。
那么我们在递推完之后,只需要把 g1 和 f1 没有加上去的前缀和补充完整即可。


例题 SPOJ DIVCNT3

F(x)=δ0(x3)
那么有:
F(p)=4
F(pc)=3c+1
因为F(p)是常数,所以只用算0次方和。

代码
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>

using namespace std;

typedef long long LL;

const int MAXN=1000000,SQRTMAXN=1000100;
#define MMax(x,y) (((x)>(y))?(x):(y))
#define MMin(x,y) (((x)<(y))?(x):(y))

int Prime[MAXN/10],pn;
int D3[MAXN+10],Fir[MAXN+10];
bool vis[MAXN+10];

int CntP[SQRTMAXN],SumP1[SQRTMAXN];
int L0[SQRTMAXN],L1[SQRTMAXN];
LL g0[SQRTMAXN],g1[SQRTMAXN],f0[SQRTMAXN],f1[SQRTMAXN];

void GetPrime();
void GetG(int SqrtN,int PN,int n);
void GetF(int SqrtN,int PN,int n);
LL GetSum(LL n);

int main()
{
    freopen("divcnt3.in","r",stdin);
    freopen("divcnt3.out","w",stdout);

    GetPrime();
    int T;
    scanf("%d",&T);
    while (T--)
    {
        LL n;
        scanf("%lld",&n);
        printf("%lld\n",GetSum(n));
    }

    return 0;
}

void GetPrime()
{
    D3[1]=1;
    for (int i=2;i<=MAXN;++i)
    {
        if (!vis[i])
            Prime[pn++]=i,D3[i]=4,Fir[i]=i;
        for (int j=0;j<pn;++j)
        {
            if (i*Prime[j]>MAXN) break;
            vis[i*Prime[j]]=1;
            if (i%Prime[j]==0) 
            {
                if (i/Fir[i]==1) D3[i*Prime[j]]=D3[i]+3;
                else D3[i*Prime[j]]=D3[i/Fir[i]]*D3[Prime[j]*Fir[i]];
                Fir[i*Prime[j]]=Fir[i]*Prime[j];break;
            }
                else D3[i*Prime[j]]=D3[i]*4,Fir[i*Prime[j]]=Prime[j];
        }
    }
}
void GetG(int SqrtN,int PN,LL n)
{
    for (int i=1;i<SqrtN;++i) g0[i]=i,g1[i]=n/i;
    for (int i=0;i<PN;++i)
    {
        for (int j=1;j<SqrtN && i<L1[j];++j)
        {
            LL y=n/j/Prime[i];
            g1[j]-=((y<SqrtN)?
                (g0[y]-MMax(0,MMin(i,CntP[y])-L0[y])):
                (g1[n/y]-MMax(0,i-L1[n/y])));
        }
        for (int j=SqrtN-1;j>=1 && i<L0[j];--j)
        {
            LL y=j/Prime[i];
            g0[j]-=(g0[y]-MMax(0,MMin(i,CntP[y])-L0[y]));
        }
    }
    for (int i=1;i<SqrtN;++i)
        g0[i]-=CntP[i]-L0[i]+1,
        g1[i]-=MMax(0,PN-L1[i])+1;
    for (int i=1;i<SqrtN;++i) g0[i]*=4,g1[i]*=4;
}
void GetF(int SqrtN,int PN,LL n)
{
    for (int i=1;i<SqrtN;++i) f0[i]=f1[i]=1;
    for (int i=PN-1;i>=0;--i)
    {
        for (int j=1;j<SqrtN && i<L1[j];++j)
        {
            LL y=n/j/Prime[i];
            for (int c=1;y;y /=Prime[i],++c)
                f1[j]+=(3*c+1)*((y<SqrtN)?
                    (f0[y]+4*MMax(0,CntP[y]-MMax(i+1,L0[y]))):
                    (f1[n/y]+4*MMax(0,PN-MMax(i+1,L1[n/y]))));
        }
        for (int j=SqrtN-1;j>=1 && i<L0[j];--j)
        {
            int y=j/Prime[i];
            for (int c=1;y;y/=Prime[i],++c)
                f0[j]+=
                    (3*c+1)*
                    (f0[y]+4*MMax(0,CntP[y]-MMax(i+1,L0[y])));
        }
    }
    for (int i=1;i<SqrtN;++i)
        f0[i]+=4*MMax(0,CntP[i]-L0[i]),
        f1[i]+=4*(PN-L1[i]);
}
LL GetSum(LL n)
{
    int SqrtN=1,PN=0;
    for (;LL(SqrtN)*SqrtN<=n;++SqrtN);
    for (;LL(Prime[PN])*Prime[PN]<=n;++PN);
    L0[0]=0;
    for (int i=1;i<SqrtN;++i) for (L0[i]=L0[i-1];LL(Prime[L0[i]])*Prime[L0[i]]<=i;++L0[i]);
    L1[SqrtN]=0;
    for (int i=SqrtN-1;i>=1;--i) 
    {
        LL x=n/i;
        for (L1[i]=L1[i+1];LL(Prime[L1[i]])*Prime[L1[i]]<=x;++L1[i]);
    }
    CntP[0]=0;
    for (int i=1;i<SqrtN;++i)
        for (CntP[i]=CntP[i-1];Prime[CntP[i]]<=i;++CntP[i]);

    GetG(SqrtN,PN,n);GetF(SqrtN,PN,n);
    LL Ans=f1[1];
    for (LL i=1;i<SqrtN;++i) Ans+=D3[i]*g1[i];
    return Ans;
}
By SemiWaker
  • 10
    点赞
  • 21
    收藏
    觉得还不错? 一键收藏
  • 9
    评论
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值