[hdu5901]Count primes
最简单的是利用Min25筛求$h(n)$的过程,即
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 1000005 4 #define ll long long 5 int K,vis[N],P[N],h[N<<1]; 6 ll n; 7 int id(ll k){ 8 if (k<=K)return k; 9 return 2*K+1-n/k; 10 } 11 ll get(int k){ 12 if (k<=K)return k; 13 return n/(2*K+1-k); 14 } 15 void init(int K){ 16 for(int i=2;i<=K;i++){ 17 if (!vis[i])P[++P[0]]=i; 18 for(int j=1;(j<=P[0])&&(i*P[j]<=K);j++){ 19 vis[i*P[j]]=1; 20 if (i%P[j]==0)break; 21 } 22 } 23 } 24 int main(){ 25 init(N-5); 26 while (scanf("%lld",&n)!=EOF){ 27 K=(int)sqrt(n); 28 for(int i=1;i<=2*K;i++)h[i]=get(i); 29 for(int i=1;P[i]<=K;i++) 30 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)]; 31 printf("%lld\n",h[id(n)]-1); 32 } 33 }
事实上,还可以进一步优化
定义$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})$的算法有较大的改善,但会被卡空间
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 3000005 4 #define ll long long 5 vector<int>v[N]; 6 int K,B,P[N],vis[N],pi[N],mu[N],Min[N],Max[N],f[N]; 7 ll n,ans,h[N]; 8 ll sqr(int k){ 9 return 1LL*k*k; 10 } 11 ll cube(int k){ 12 return 1LL*k*k*k; 13 } 14 int id(ll k){ 15 if (k<=K)return k; 16 return 2*K+1-n/k; 17 } 18 ll get(int k){ 19 if (k<=K)return k; 20 return n/(2*K+1-k); 21 } 22 void init(int K){ 23 mu[1]=1; 24 for(int i=2;i<=K;i++){ 25 pi[i]=pi[i-1]; 26 if (!vis[i]){ 27 P[++P[0]]=i; 28 mu[i]=-1; 29 pi[i]++; 30 Min[i]=Max[i]=i; 31 } 32 for(int j=1;(j<=P[0])&&(i*P[j]<=K);j++){ 33 vis[i*P[j]]=1; 34 Min[i*P[j]]=P[j]; 35 Max[i*P[j]]=Max[i]; 36 if (i%P[j])mu[i*P[j]]=-mu[i]; 37 else{ 38 mu[i*P[j]]=0; 39 break; 40 } 41 } 42 } 43 } 44 int lowbit(int k){ 45 return (k&(-k)); 46 } 47 void update(int k){ 48 while (1LL*k*P[B]<=n){ 49 f[k]++; 50 k+=lowbit(k); 51 } 52 } 53 int query(int k){ 54 int ans=0; 55 while (k>0){ 56 ans+=f[k]; 57 k-=lowbit(k); 58 } 59 return ans; 60 } 61 int query(int l,int r){ 62 return query(r)-query(l-1); 63 } 64 int main(){ 65 init(N-5); 66 while (scanf("%lld",&n)!=EOF){ 67 K=3e6; 68 if (n<=K){ 69 printf("%d\n",pi[n]); 70 continue; 71 } 72 int pB=(int)pow(n*pow(log(n),2),1.0/3); 73 for(int i=1;i<=P[0];i++) 74 if (P[i]>=pB){ 75 B=i; 76 break; 77 } 78 ans=pi[P[B]]-1; 79 for(int i=B+1,j=P[0];sqr(P[i])<=n;i++){ 80 while (1LL*P[i]*P[j]>n)j--; 81 ans-=pi[P[j]]-pi[P[i]-1]; 82 } 83 for(int i=1;i<=P[B];i++)ans+=mu[i]*(n/i); 84 K=(int)sqrt(n); 85 for(int i=1;i<=2*K;i++)h[i]=get(i); 86 for(int i=1;sqr(P[i])<=P[B];i++){ 87 for(int j=1;j<=P[B];j++) 88 if ((P[B]<1LL*j*P[i])&&(Min[j]>P[i]))ans-=mu[j]*h[id(n/P[i]/j)]; 89 for(int j=2*K;j;j--)h[j]-=h[id(get(j)/P[i])]; 90 } 91 for(int i=B;cube(P[i])>n;i--)ans+=pi[P[B]]-pi[P[i]]; 92 int l=0,r=0; 93 for(int i=1;i<=B;i++){ 94 if ((!l)&&(sqr(P[i])>P[B]))l=i; 95 if (cube(P[i])<=n)r=i; 96 } 97 for(int i=l;i<=r;i++)ans+=B-i; 98 for(int i=l;i<=r;i++){ 99 ll nn=n/P[i]; 100 for(ll x=P[i],y;P[i]<nn/x;x=y+1){ 101 y=nn/(nn/x); 102 ans+=1LL*(pi[y]-pi[x-1])*(pi[min(nn/x,(ll)P[B])]-pi[P[i]]); 103 } 104 } 105 memset(f,0,sizeof(f)); 106 for(int i=1;i<=B;i++) 107 if (sqr(sqr(P[i]))<=n)r=i; 108 for(int i=1;i<=P[0];i++)v[P[i]].clear(); 109 for(int i=1;1LL*i*P[B]<=n;i++) 110 if (vis[i]){ 111 if (Min[i]>P[r])update(i); 112 else v[Min[i]].push_back(i); 113 } 114 for(int i=r;i>=l;i--){ 115 for(int j=0;j<v[P[i]].size();j++)update(v[P[i]][j]); 116 ll nn=n/P[i]; 117 for(ll x=1,y;P[i]<nn/x;x=y+1){ 118 y=nn/(nn/x); 119 ans+=1LL*query(x,y)*(pi[min(nn/x,(ll)P[B])]-pi[P[i]]); 120 } 121 } 122 printf("%lld\n",ans); 123 } 124 }