洲阁筛学习总结

洲阁筛学习总结

By Semiwaker

问题描述

给出一个积性函数$F(x)$，满足积性函数的基本性质：
F(1)=1
如果 a 和 b 互质，则$F(ab)=F(a)F(b)$
设$x=\prod p_i^{k_i}$，则有

$$F(x)=\prod F(p_i^{k_i})$$
特别的
p为质数，F(p) 是一个低阶多项式，$F(p^c)$可以快速算出。
现在要求 $$\sum_{i=1}^nF(i)$$
n的范围大概在$10^{11}$左右。

---

小范围问题的解决方法

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

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

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

---

前置技能

引理1

如果$x>a,b$
$\lfloor\frac{\lfloor\frac xa\rfloor}b\rfloor=\lfloor\frac x{ab}\rfloor$
证明：
设$y=\lfloor\frac xa\rfloor$。
设$k=\lfloor\frac x{ab}\rfloor$，则有$x=kab+c$，显然$c<ab$。
那么$y=kb+\lfloor\frac ca\rfloor$，因为$c<ab$，所以$\lfloor\frac ca\rfloor<b$，所以$\lfloor\frac yb\rfloor=k$。

引理2

$\lfloor\frac nx\rfloor$最多有$2\sqrt n$种取值。
证明
当$x\le \sqrt n$时，只有$\sqrt n$种取值。
当$x>\sqrt n$时，$\lfloor\frac nx\rfloor\le\sqrt n$，也只有$\sqrt n$种取值。

这个上限是比较精准的，大多数情况下是$2\sqrt n-1$或者刚好$2\sqrt n$个。

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

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

洲阁筛核心思想

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

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

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

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

---

算法

总述

我们把１～n的每一个数分类，有大于$\sqrt n$的质数的，和没有的。
那么可以写成这样：

$$\sum_{i=1}^nF(i)=\sum_{i=1}^{\sqrt n}F(i)\Bigl(\sum_{p>\sqrt n}^{\lfloor\frac ni\rfloor}F(p)[p为质数]\Bigr)+\sum_{i=1}^nF(i)[i无大于\sqrt n的质因子]$$
也就是说，有大于$\sqrt n$的质因子的，我们枚举$\le \sqrt n$的质因子的部分，再枚举$>\sqrt n$的质因子。
没有的就直接枚举。

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

$$\sum_{p>\sqrt n}^{\lfloor\frac ni\rfloor}F(p)[p为质数]$$
第三个部分是求 $$\sum_{i=1}^nF(i)[i无大于\sqrt n的质因子]$$

第二部分

求

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

筛的过程可用递推实现。
设小于等于$\sqrt n$的质数有Pn个，从小到大排列为$P_{1..Pn}$。
设在 [1,j] 中，与$P_{1..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个质数不互质。
$\sum (P[i]x)^k=P[i]^k\sum x$
又$P[i]x\le j$，那么有$x\le\lfloor\frac j{P[i]}\rfloor$
所以得到：

$$g[i][j]=g[i-1][j]-P_i^kg[i-1][\lfloor\frac j{P[i]}\rfloor]$$

根据之前提到的定理，j的取值有$2\sqrt n$种，i的取值有$O(\frac {\sqrt n}{ln \sqrt n})$种，乘起来还是O(n)级别。需要优化。

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

怎么利用这个性质呢？
如果$P_{i-1}>\lfloor\frac j{P[i]}\rfloor$，那么$g[i-1][\lfloor\frac j{P[i]}\rfloor]=1$.
放宽一点$P_i>\lfloor\frac j{P[i]}\rfloor$，此时的判断条件转化为$P_i^2>j$。
换句话说，
当$P_i^2>j$时，有$g[i][j]=g[i-1][j]-P_i^k$
那么设$P_{L[j]}^2>j\ge P_{L[j]-1}^2$，即L[j]是第一个满足条件的i。
对于$i<L[j]$，我们需要老老实实地去算。对于$i\ge L[j]$，我们只需要一个前缀和就可以了。

还有一点要注意的：
如果$P_i>j$，此时$\lfloor\frac j{P[i]}\rfloor=0$，那么有

$$g[i][j]=g[i-1][j]$$
设小于等于j的质数有CntP[j]个，那么对于$i\ge CntP[j]$，是不需要去计算的。

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

观察到递推式只和i-1有关，是可以滚动数组的，所以具体的实现方法如下：
当$i<L[j]$时，直接$g[j]=g[j]-P_i^kg[\lfloor\frac j{P[i]}\rfloor]$
当$i\ge L[j]$时，停止对$g[j]$的递推。
如果在某个i调用到了g[j]，那么我们再算上缺少的L[j]~min(i,CntP[j])的部分的前缀和即可。
设前缀和为 $SumP[i]=\sum_{t=0}^ip_t^k$ ，那么我们可以这样表示 g[i][j]：
$g[j]+max(0,SumP[min(i,CntP[j])]-SumP[L[j]])$

时间复杂度

由于每一个j对应的需要枚举的i是$\frac {\sqrt j}{ln j}$范围的，所以我们可以列式：

$$\sum_{i=1}^{\sqrt n}\frac {\sqrt \frac ni}{ln \frac ni}\approx O(\frac {n^{\frac 34}}{ln n})$$

细节

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

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

第三部分

求

$$\sum_{i=1}^nF(i)[i无大于\sqrt n的质因子]$$

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

没有大于$\sqrt n$的质因子，那么就意味着只有$\sqrt n$以内的质因子。
所以我们设：
x为在[1,j]的范围内，只有$P_{i..Pn}$这些质因子的数。（即i以后的质因子。）
$f[i][j]=\sum F(x)$

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

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

那么有：

$$f[i][j]=f[i+1][j]+\sum F(P_i^c)f[i+1][\lfloor\frac j{P_i^c}\rfloor]$$
简单来说，直接枚举当前质因子的幂，然后考虑乘起来不会超过j的，只包含$P_{i+1..Pn}$的数的和为多少。因为是积性函数，所以可直接乘起来。

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

同样考虑$P_i^2>j$，则$\lfloor\frac j{P_i}\rfloor<P_i$，此时$f[i+1][\lfloor\frac j{P_i}\rfloor]=1$（因为只包含大于Pi的质数）。
而当$c\ge2$，$\lfloor\frac j{P_i^c}\rfloor=0$。
所以有

$$f[i][j]=f[i+1][j]+F(P_i)$$

同样当$P_i>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]=\sum_{t=1}^iF(t)$
那么，在某个 i 时，f[j]的实际值为 f[j]+Max(0,SumP[L[i]-1]-SumP[Max(i,CntP[i])-1])

时间复杂度

同样为

$$O(\frac {n^{\frac 34}}{ln n})$$

细节

同样，当$j\le\sqrt n$时，记为 f0[j]，否则设$j=\lfloor\frac nx\rfloor$，记为f1[j]

计算答案

现在回到原本的式子：

$$\sum_{i=1}^nF(i)=\sum_{i=1}^{\sqrt n}F(i)\Bigl(\sum_{p>\sqrt n}^{\lfloor\frac ni\rfloor}F(p)[p为质数]\Bigr)+\sum_{i=1}^nF(i)[i无大于\sqrt n的质因子]$$
实际需要计算的是：
$Ans=\sum_{i=1}^{\sqrt n}F(i)g(\lfloor\frac ni\rfloor)+f(n)$
即
$Ans=\sum_{i=1}^{\sqrt n}F(i)g1(i)+f1(1)$

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

---

例题 SPOJ DIVCNT3

$F(x)=\delta_0(x^3)$
那么有：
$F(p)=4$
$F(p^c)=3c+1$
因为F(p)是常数，所以只用算０次方和。

代码

#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