[hdu5901]Count primes

最简单的是利用Min25筛求$h(n)$的过程，即

#include<bits/stdc++.h>
using namespace std;
#define N 1000005
#define ll long long
int K,vis[N],P[N],h[N<<1];
ll n;
int id(ll k){
    if (k<=K)return k;
    return 2*K+1-n/k;
}
ll get(int k){
    if (k<=K)return k;
    return n/(2*K+1-k);
}
void init(int K){
    for(int i=2;i<=K;i++){
        if (!vis[i])P[++P[0]]=i;
        for(int j=1;(j<=P[0])&&(i*P[j]<=K);j++){
            vis[i*P[j]]=1;
            if (i%P[j]==0)break;
        }
    }
}
int main(){
    init(N-5);
    while (scanf("%lld",&n)!=EOF){
        K=(int)sqrt(n);
        for(int i=1;i<=2*K;i++)h[i]=get(i);
        for(int i=1;P[i]<=K;i++)
            for(int j=2*K;(j)&&(1LL*P[i]*P[i]<=get(j));j--)h[j]-=h[id(get(j)/P[i])]-h[id(P[i]-1)];
        printf("%lld\n",h[id(n)]-1);
    }
} 

事实上，还可以进一步优化

定义$Min(i)$为$i$的最小素因子（特别的，$Min(1)=\infty$），$Max(i)$为$i$的最大素因子，$p_{i}$为从小到大第$i$个素数，$\pi_{k}(n)=\sum_{p_{i}\le n}p_{i}^{k}$（本题中$k=0$，但由于$k$并不影响下述算法，因此以下设$k$为较小的常数）

记$cnt(x)$为$x$的可重素因子个数，具体来说，满足：

1.$cnt(1)=0$，$\forall p_{i},cnt(p_{i})=1$

2.$\forall x,y>0,cnt(xy)=cnt(x)+cnt(y)$

选择一个素数$p_{B}$，满足$p_{B}^{3}\ge N$且$p_{B}^{2}\le N$，关于$P_{B}$最为合适的值以后再分析

定义$\phi_{s}(n,m)=\sum_{1\le i\le n,p_{m}<Min(i),cnt(i)=s}i^{k}$，$\phi(n,m)=\sum_{s}\phi_{s}(n,m)$

当$p_{m+1}^{s}>n$时有$\phi_{s}(n,m)=0$，因此，不难得到有$\phi(N,B)=\sum_{s=0}^{2}\phi_{s}(N,B)$，分别考虑这三项，又有
$$
\phi(N,B)=1+(\pi_{k}(N)-\pi_{k}(P_{B}))+\phi_{2}(N,B)
$$
简单移项，又有
$$
\pi_{k}(N)=\pi_{k}(P_{B})-\phi_{2}(N,B)+\phi(N,B)-1
$$
关于$\pi_{k}(P_{B})$可以$o(P_{B})$的求出，$\phi_{2}(N,B)$不妨暴力枚举具体的两个素数，即
$$
\phi_{2}(N,B)=\sum_{B<i\le j,p_{i}p_{j}\le N}(p_{i}p_{j})^{k}=\sum_{B<i}p_{i}^{k}\sum_{i\le j,p_{i}\cdot p_{j}\le N}p_{j}^{k}
$$
由于$i>B$，后者即$\pi_{k}$的前缀和差分，用$o(\frac{N}{P_{B}})$的复杂度预处理，同时$p_{i}^{2}\le N$，关于$i$的枚举即$o(\sqrt{N})$

$\phi(n,m)$的计算-递推

类似前面关于$\phi(n,m)$的计算，不难得到$\phi(n,m)=\phi(n,m-1)-p_{m}^{k}\phi(\lfloor\frac{n}{p_{m}}\rfloor,m-1)$

考虑递归计算，那么就构成了一棵满二叉树，更具体的，有：

1.用二元组表示一个节点，通常记作$(x,h)$（$h\ge 0$），根为$(1,m)$，且边有边权

2.对于$(x,h)$的状态，分类讨论——

(1)若$h=0$，则$(x,h)$的为叶子；

(2)若$h>0$，则$(x,h)$的左儿子为$(x,h-1)$，边权为1，右儿子为$(x\cdot p_{h},h-1)$，边权为$-p_{h}^{k}$

定义$(x,h)$的权值为其到根路径上所有边边权乘积$\times \phi(\lfloor\frac{n}{x}\rfloor,h)$，根据$\phi$的递推式，一个点的权值即为其两个儿子的权值和（都乘了一个$x$到根的系数）

根据这个，我们只需要选择若干个节点，使得每一个叶子到根路径上恰有一个节点被选择，那么不难证明根节点的权值即为所有选择的节点权值和

下面，我们构造一种选择方式，并简单证明一下这样选择的正确性——

选择$(x,h)$当且仅当$h=0$且$x\le p_{B}$，或$p_{B}<x\le p_{B}\cdot Min(x)$且$p_{h+1}=Min(x)$

根据前面，也就是要求叶子到根的路径恰好存在一个节点被选择

从构造树的过程，可以发现：对于一个非根节点$(x,h)$，若$p_{h+1}=Min(x)$则其父亲为$(\frac{x}{Min(x)},h+1)$，否则为$(x,h+1)$

1.对于$(x,0)$（其中$x\le p_{B}$），首先其必然直接被选择，同时对于其到祖先的链上节点权值必然因此不大于$x$即不大于$p_{B}$，不可能有其余点被选

2.对于$(x,0)$（其中$x>p_{B}$），不断向上找到其中第一个节点$y$使得$y$的父亲不超过$p_{B}$，由于$y$的父亲要比$y$小（否则$y$的儿子也可以），那么$y$的父亲即$\frac{y}{Min(y)}$，那么$p_{h+1}=Min(y)$

同时，$y$自己是可以的，因此就得到了$p_{B}<y\le p_{B}\cdot Min(y)$

对于这个右儿子的父亲以上的节点与第一点相同，都不超过$p_{B}$，右儿子以下的节点中分类讨论：

(1)若等于$y$，那么$h$减小，不满足$p_{h+1}=Min(y)$

(2)若不等于$y$（大于），记作$z$，来考虑第一个比其小的祖先，其权值不小于$y$（否则可以选$y$），不难证明这个节点就是$\frac{z}{Min(z)}\ge y>p_{B}$，不符合条件

接下来，考虑一个$x$所对应的系数：首先，每一次乘上素数$p$都有一个$p^{k}$，总共即为$x^{k}$，符号取决于素因子个数，恰好是$\mu(x)$

另外还有$x$的条件，在上面的基础上，我们还需要保证$(x,h)$这个节点会出现，即$Max(x)\le p_{m}$以及$x$中不含有素数平方的因子，后者由于$\mu(x)$恰好为0，可以不判定

更具体的，有下式
$$
\phi(n,m)=\sum_{1\le x\le p_{B},Max(x)\le p_{m}}\mu(x)x^{k}\phi(\lfloor\frac{n}{x}\rfloor,0)+\sum_{p_{B}<x\le p_{B}\cdot Min(x),Max(x)\le p_{m}}\mu(x)x^{k}\phi(\lfloor\frac{n}{x}\rfloor,\pi(Min(x))-1)
$$

$\phi(n,m)$的计算-分类

关于$\phi(n,0)=\sum_{i=1}^{n}i^{k}$不难推导出通项公式，因此前者可以$o(p_{B})$计算

进一步的，将后者拆分，先枚举$p_{m'}=Min(x)$，即
$$
-\sum_{m'=1}^{m}p_{m'}^{k}\sum_{x\le p_{B}<x\cdot p_{m'},Min(x)>p_{m'},Max(x)\le p_{m}}\mu(x)x^{k}\phi(\lfloor\frac{n}{x\cdot p_{m'}}\rfloor,m'-1)
$$
具体代入，所求的是$\phi(N,B)$，即
$$
-\sum_{m'=1}^{B}p_{m'}^{k}\sum_{x\le p_{B}<x\cdot p_{m'},Min(x)>p_{m'}}\mu(x)x^{k}\phi(\lfloor\frac{N}{x\cdot p_{m'}}\rfloor,m'-1)
$$
分为三类来计算——

1.当$p_{m'}^{2}\le p_{B}$，此时对于后面的状态，由于$\lfloor\frac{n}{x\cdot p_{m'}}\rfloor\in S_{n}$，可以在$o(\frac{\sqrt{p_{B}N}}{\log n})$的时间内预处理出来（根据$\phi(n,m)=\phi(n,m-1)-p_{m}^{k}\phi(\lfloor\frac{n}{p_{m}}\rfloor,m-1)$来计算）

再以$o(\frac{p_{B}\sqrt{p_{B}}}{\log p_{B}})$的复杂度枚举$m'$和$x$，总复杂度即为$o(\frac{\sqrt{p_{B}N}}{\log N})$（枚举$m'$和$x$的复杂度较小）

2.当$p_{m'}^{2}>p_{B}$，此时考虑$x\le p_{B}$且$Min(x)>p_{m'}$，不难发现$cnt(x)<2$，即$x$为素数或是1，后者显然不满足$p_{B}<x\cdot p_{m'}$，因此$x$为素数

接下来，原式即
$$
\sum_{m'\le B,p_{m'}^{2}> p_{B}}p_{m'}^{k}\sum_{x=m'+1}^{B}p_{x}^{k}\phi(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)
$$
$\phi(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)=\sum_{s}\phi_{s}(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)$，又有当$\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor<p^{s}_{m'}$时为0，那么考虑令$s$只能为0来缩小$p_{m'}$的范围，即$p_{x}\cdot p_{m'}^{2}>N$，又因为$p_{x}>p_{m'}$，放缩为$p_{m'}^{3}>N$

此时该式即为$\phi_{0}(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)=[N\ge p_{x}\cdot p_{m'}]$，又因为$p_{B}^{2}\le N$，恒为1，即
$$
\sum_{m'\le B,p_{m'}^{3}>N}p_{m'}^{k}\sum_{m'<x\le B}p_{x}^{k}
$$
对$m'$用$o(B)$枚举，后面即一个$\pi_{k}$的前缀和差分，预处理即可，总复杂度为$o(B)$

3.$p_{m'}^{2}>p_{B}$且$p_{m'}^{3}\le N$，仍然有$x$为素数，与第二部分的第一个式子相同，即
$$
\sum_{p_{m'}^{2}>p_{B},p_{m'}^{3}\le N}p_{m'}^{k}\sum_{m'<x\le B}p_{x}^{k}\phi(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)
$$
这一类比较难算，回顾$\phi(\lfloor\frac{N}{p_{x}\cdot p_{m'}}\rfloor,m'-1)$的定义，由于$cnt(x)$是任意的，即在$\phi_{s}(n,m)$的基础上去掉$s$的限制即可，将其代入即
$$
\sum_{p_{m'}^{2}>p_{B},p_{m'}^{3}\le N}p_{m'}^{k}\sum_{x=m'+1}^{B}p_{x}^{k}\sum_{q\cdot p_{x}\cdot p_{m'}\le N,Min(q)\ge p_{m'}}q^{k}
$$
调换$x$和$q$的枚举顺序，即
$$
\sum_{p_{m'}^{2}>p_{B},p_{m'}^{3}\le N}p_{m'}^{k}\sum_{q\cdot p_{m'}^{2}\le N,Min(q)\ge p_{m'}}q^{k}\sum_{m'<x\le B,q\cdot p_{x}\cdot p_{m'}\le N}p_{x}^{k}
$$
关于这件事情的计算方式：

先枚举$m'$，之后对$\lfloor\frac{N}{p_{m'}}\rfloor$数论分块，通过$m'$和$\lfloor\frac{N}{q\cdot p_{m'}}\rfloor$显然可以确定$x$这一维度的和（$o(B)$预处理），接下来也就是求区间内有多少给$q$满足$Min(q)\ge p_{m'}$

由于$p_{m'}^{2}>p_{B}$，即$q\cdot p_{B}\le N$，可以预处理出每一个$q$的$Min(q)$，并从大到小枚举$m'$，加入对应的$q$，再用树状数组来维护即可

时间复杂度上，对$m'$的枚举是$o(\frac{N^{\frac{1}{3}}}{\log N})$，接下来数论分块是$\sqrt{\frac{N}{p_{m'}}}$，将其求和即
$$
\sum_{p_{m'}^{2}>p_{B},p_{m'}^{3}\le N}\sqrt{\frac{N}{p_{m'}}}=\sqrt{N}\sum_{1\le i^{3}\le N}(\pi(i)-\pi(i-1))i^{-\frac{1}{2}}
$$
根据素数定理，即$\pi(n)=o(\frac{n}{\log n})$，代入后即为
$$
\frac{\sqrt{N}}{\log N}\sum_{1\le i\le N^{\frac{1}{3}}}i^{-\frac{1}{2}}=\frac{\sqrt{N}}{\log N}\int_{1}^{N^{\frac{1}{3}}}i^{-\frac{1}{2}}=\frac{N^{\frac{2}{3}}}{\log N}
$$
之后还有一个树状数组，总复杂度为$o(N^{\frac{2}{3}})$，理论上优于$o(\frac{N^{\frac{3}{4}}}{\log N})$

同时，对于其修改总复杂度为$o(\frac{N\log N}{p_{B}})$，这取决于需要插入的$q$数量

$\phi(n,m)$的计算-优化

发现复杂度的瓶颈在于第三部分，考虑再对$q$分类讨论：

1.若$q=1$，即为
$$
\sum_{p_{m'}^{2}>p_{B},p_{m'}^{3}\le N}p_{m'}^{k}\sum_{m'<x\le B}p_{x}^{k}
$$
随便维护一下，复杂度为$o(\frac{N^{\frac{1}{3}}}{\log N})$

2.若$q$是素数，则即要求$q\ge p_{m'}$，这个就不需要用树状数组维护，即$o(\frac{N^{\frac{2}{3}}}{\log N})$

3.若$q$不是素数，即$q\ge Min(q)^{2}\ge p_{m'}^{2}$，即$p_{m'}^{4}\le N$，此时询问复杂度类似于前面是$o(N^{\frac{5}{8}})$

对于修改的复杂度，事实上是$o(\frac{N}{p_{B}})$的，以下来证明：

首先，有这样一个性质$\forall s\ge 1,\sum_{i=1}^{n}[cnt(i)=s]=o(\frac{n}{\log n})$

证明考虑归纳，$s=1$即素数定理，$s>1$时有
$$
\sum_{i=1}^{n}[cnt(i)=s]=\sum_{p_{i}\le n}\sum_{1\le j\le \lfloor\frac{n}{p_{i}}\rfloor}[cnt(j)=s-1]=\frac{\sum_{p_{i}\le n}{\lfloor\frac{n}{p_{i}}}\rfloor}{\log n}
$$
对于分子简单来看，就是调和级数*素数的概率，后者根据素数定理即$o(\frac{1}{\log n})$，总共即$o(n)$（更精确的来说是埃氏筛的复杂度，即为$o(n\log\log n)$，这里就简单看作$o(n)$了）

那么即证明了对于$s>1$，对应的个数也是$o(\frac{n}{\log n})$的

事实上，这件事情的正确理解方式应该是对于任意$s$，总存在足够大的$n$，使得此时$cnt(i)=s$的数量大约为$o(\frac{n}{\log n})$的的级别

考虑$p_{B}^{3}\ge N$，且$Min(q)^{2}\ge p_{m'}^{2}>p_{B}$，有$Min(q)^{6}\ge N$，因此$cnt(q)\le 6$，数量为$o(\frac{N}{p_{B}\log N})$（由于$N$并不足够大，这里的范围实际更接近$o(\frac{N}{p_{B}})$），算上树状数组后即$o(\frac{N}{p_{B}})$

时间复杂度

综合上述，复杂度包含——

0.线性筛，$o(\frac{N}{p_{B}})$预处理出素数（个数）、$\mu(i)$、$\pi_{k}(i)$、$Min(i)$和$Max(i)$（$i\cdot p_{B}\le N$）

1.$o(1)$求出$\pi_{k}(p_{B})$

2.$o(\sqrt{N})$求出$\phi_{2}(N,B)$

3.$o(\frac{\sqrt{p_{B}N}}{\log N})$求出$\phi(N,B)$的第一部分

4.$o(B)$求出$\phi(N,B)$的第二部分

5.$o(\frac{N^{\frac{2}{3}}}{\log N}+N^{\frac{5}{8}}+\frac{N}{p_{B}})$求出$\phi(N,B)$的第三部分

关于$p_{B}$的取值考虑$o(\frac{\sqrt{p_{B}N}}{\log N})$和$o(\frac{N}{p_{B}})$，相等即解出$p_{B}=N^{\frac{1}{3}}\log^{\frac{2}{3}}N$，总时间复杂度为$o((\frac{N}{\log N})^{\frac{2}{3}})$

空间复杂度与时间复杂度相同，也为$o((\frac{N}{\log N})^{\frac{2}{3}})$，劣于$o(\sqrt{N})$的空间复杂度

具体实现中，这个算法的时间相比$o(\frac{N^{\frac{3}{4}}}{\log N})$的算法有较大的改善，但会被卡空间

#include<bits/stdc++.h>
using namespace std;
#define N 3000005
#define ll long long
vector<int>v[N];
int K,B,P[N],vis[N],pi[N],mu[N],Min[N],Max[N],f[N];
ll n,ans,h[N];
ll sqr(int k){
    return 1LL*k*k;
}
ll cube(int k){
    return 1LL*k*k*k;
}
int id(ll k){
    if (k<=K)return k;
    return 2*K+1-n/k;
}
ll get(int k){
    if (k<=K)return k;
    return n/(2*K+1-k);
}
void init(int K){
    mu[1]=1;
    for(int i=2;i<=K;i++){
        pi[i]=pi[i-1];
        if (!vis[i]){
            P[++P[0]]=i;
            mu[i]=-1;
            pi[i]++;
            Min[i]=Max[i]=i;
        }
        for(int j=1;(j<=P[0])&&(i*P[j]<=K);j++){
            vis[i*P[j]]=1;
            Min[i*P[j]]=P[j];
            Max[i*P[j]]=Max[i];
            if (i%P[j])mu[i*P[j]]=-mu[i];
            else{
                mu[i*P[j]]=0;
                break;
            }
        }
    }
}
int lowbit(int k){
    return (k&(-k));
}
void update(int k){
    while (1LL*k*P[B]<=n){
        f[k]++;
        k+=lowbit(k);
    }
}
int query(int k){
    int ans=0;
    while (k>0){
        ans+=f[k];
        k-=lowbit(k);
    }
    return ans;
}
int query(int l,int r){
    return query(r)-query(l-1);
}
int main(){
    init(N-5);
    while (scanf("%lld",&n)!=EOF){
        K=3e6;
        if (n<=K){
            printf("%d\n",pi[n]);
            continue;
        } 
        int pB=(int)pow(n*pow(log(n),2),1.0/3);
        for(int i=1;i<=P[0];i++)
            if (P[i]>=pB){
                B=i;
                break;
            }
        ans=pi[P[B]]-1;
        for(int i=B+1,j=P[0];sqr(P[i])<=n;i++){
            while (1LL*P[i]*P[j]>n)j--;
            ans-=pi[P[j]]-pi[P[i]-1];
        }
        for(int i=1;i<=P[B];i++)ans+=mu[i]*(n/i);
        K=(int)sqrt(n);
        for(int i=1;i<=2*K;i++)h[i]=get(i);
        for(int i=1;sqr(P[i])<=P[B];i++){
            for(int j=1;j<=P[B];j++)
                if ((P[B]<1LL*j*P[i])&&(Min[j]>P[i]))ans-=mu[j]*h[id(n/P[i]/j)];
            for(int j=2*K;j;j--)h[j]-=h[id(get(j)/P[i])];
        }
        for(int i=B;cube(P[i])>n;i--)ans+=pi[P[B]]-pi[P[i]];
        int l=0,r=0;
        for(int i=1;i<=B;i++){
            if ((!l)&&(sqr(P[i])>P[B]))l=i;
            if (cube(P[i])<=n)r=i;
        }
        for(int i=l;i<=r;i++)ans+=B-i;
        for(int i=l;i<=r;i++){
            ll nn=n/P[i];
            for(ll x=P[i],y;P[i]<nn/x;x=y+1){
                y=nn/(nn/x);
                ans+=1LL*(pi[y]-pi[x-1])*(pi[min(nn/x,(ll)P[B])]-pi[P[i]]);
            }
        }
        memset(f,0,sizeof(f));
        for(int i=1;i<=B;i++)
            if (sqr(sqr(P[i]))<=n)r=i;
        for(int i=1;i<=P[0];i++)v[P[i]].clear();
        for(int i=1;1LL*i*P[B]<=n;i++)
            if (vis[i]){
                if (Min[i]>P[r])update(i);
                else v[Min[i]].push_back(i);
            }
        for(int i=r;i>=l;i--){
            for(int j=0;j<v[P[i]].size();j++)update(v[P[i]][j]);
            ll nn=n/P[i];
            for(ll x=1,y;P[i]<nn/x;x=y+1){
                y=nn/(nn/x);
                ans+=1LL*query(x,y)*(pi[min(nn/x,(ll)P[B])]-pi[P[i]]);
            }
        }
        printf("%lld\n",ans);
    }
}