[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
",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_{1le ile 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<ile j,p_{i}p_{j}le N}(p_{i}p_{j})^{k}=sum_{B<i}p_{i}^{k}sum_{ile 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(lfloorfrac{n}{p_{m}} floor,m-1)$

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

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

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

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

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

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

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

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

选择$(x,h)$当且仅当$h=0$且$xle p_{B}$，或$p_{B}<xle 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)$（其中$xle 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}<yle 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_{1le xle p_{B},Max(x)le p_{m}}mu(x)x^{k}phi(lfloorfrac{n}{x} floor,0)+sum_{p_{B}<xle p_{B}cdot Min(x),Max(x)le p_{m}}mu(x)x^{k}phi(lfloorfrac{n}{x} floor,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_{xle p_{B}<xcdot p_{m'},Min(x)>p_{m'},Max(x)le p_{m}}mu(x)x^{k}phi(lfloorfrac{n}{xcdot p_{m'}} floor,m'-1)
$$
具体代入，所求的是$phi(N,B)$，即
$$
-sum_{m'=1}^{B}p_{m'}^{k}sum_{xle p_{B}<xcdot p_{m'},Min(x)>p_{m'}}mu(x)x^{k}phi(lfloorfrac{N}{xcdot p_{m'}} floor,m'-1)
$$
分为三类来计算——

1.当$p_{m'}^{2}le p_{B}$，此时对于后面的状态，由于$lfloorfrac{n}{xcdot p_{m'}} floorin S_{n}$，可以在$o(frac{sqrt{p_{B}N}}{log n})$的时间内预处理出来（根据$phi(n,m)=phi(n,m-1)-p_{m}^{k}phi(lfloorfrac{n}{p_{m}} floor,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}$，此时考虑$xle p_{B}$且$Min(x)>p_{m'}$，不难发现$cnt(x)<2$，即$x$为素数或是1，后者显然不满足$p_{B}<xcdot p_{m'}$，因此$x$为素数

接下来，原式即
$$
sum_{m'le B,p_{m'}^{2}> p_{B}}p_{m'}^{k}sum_{x=m'+1}^{B}p_{x}^{k}phi(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,m'-1)
$$
$phi(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,m'-1)=sum_{s}phi_{s}(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,m'-1)$，又有当$lfloorfrac{N}{p_{x}cdot p_{m'}} floor<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}(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,m'-1)=[Nge p_{x}cdot p_{m'}]$，又因为$p_{B}^{2}le N$，恒为1，即
$$
sum_{m'le B,p_{m'}^{3}>N}p_{m'}^{k}sum_{m'<xle 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'<xle B}p_{x}^{k}phi(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,m'-1)
$$
这一类比较难算，回顾$phi(lfloorfrac{N}{p_{x}cdot p_{m'}} floor,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_{qcdot 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_{qcdot p_{m'}^{2}le N,Min(q)ge p_{m'}}q^{k}sum_{m'<xle B,qcdot p_{x}cdot p_{m'}le N}p_{x}^{k}
$$
关于这件事情的计算方式：

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

由于$p_{m'}^{2}>p_{B}$，即$qcdot 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_{1le 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_{1le ile 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{Nlog 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'<xle B}p_{x}^{k}
$$
随便维护一下，复杂度为$o(frac{N^{frac{1}{3}}}{log N})$

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

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

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

首先，有这样一个性质$forall sge 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_{1le jle lfloorfrac{n}{p_{i}} floor}[cnt(j)=s-1]=frac{sum_{p_{i}le n}{lfloorfrac{n}{p_{i}}} floor}{log n}
$$
对于分子简单来看，就是调和级数*素数的概率，后者根据素数定理即$o(frac{1}{log n})$，总共即$o(n)$（更精确的来说是埃氏筛的复杂度，即为$o(nloglog 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)$（$icdot 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
",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
",ans);
    }
}