一种计算 π(n) 的组合方法
这里
π(n)
指:不大于
n
的素数个数
注:本文方法来自于维基百科对素数计数的一个组合方法:https://en.wikipedia.org/wiki/Prime-counting_function
由于最后的一些优化技巧还未掌握。文中的方法还有待优化。
与原文的方法略微不同。与洲阁筛的方法近似。。。
常识:一般来说,x 以内的素数大约有:
π(x)=O(xln x)
那么我们可以猜测埃式筛法的复杂度:
T(n)=O(n∑p≤n1p)
可以认为是线性的。
在
n
不是很大的时候。埃式筛法是个不错的选择。(不知道素数筛法的可以自行百度。)
π(x)=O(xln x)
T(n)=O(n∑p≤n1p)
埃式筛法在计算π(n) 时。还会得到素数表。
像埃式筛法这样:通过得到所有不大于
n
的素数,来计算π(n) 。算法复杂度的下界是
O(nlogn)
。不会很快。
<script type="math/tex; mode=display" id="MathJax-Element-12"></script>
组合方法:
对于一个整数
q
,q∈(n‾√,n]
q
如果不能被[2,n‾√] 的任何整数整除。
则:
q
是一个素。
很容易得到一个基于容斥原理的式子:
[1,n] 中,不能被前
m
个素数整除的数字个数为:
C(m,n)=⌊n⌋−∑1≤i≤m⌊nPi⌋+∑1≤i<j≤m⌊nPiPj⌋−∑1≤i<j<k≤m⌊nPiPjPk⌋...
Pi
指第
i
个素数
则:
C(π(n‾‾√),n)=π(n)−π(n‾‾√)+1
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)
对于第二个式子。
[1,n]
中不能被前
m−1
个素数整除的数字包含:
可以被
Pm
整除的数字。但不能被前
m−1
个素数整除的数字
对于
Pm
的所有倍数:
kPm, 其中 k∈[1,⌊nPm⌋]
不能被前
m−1
个素数整除的
k
的数量为:
C(m−1,⌊nPm⌋)
所以:
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)
特别的:
C(0,n)=⌊n⌋
方便起见,下文所有数均为整数。
对于
C(m,n)
定义:
D(m,d)=C(m,⌊nd⌋)
因为:
⌊⌊na⌋b⌋=⌊nab⌋
所以:
D(m,d)=D(m−1,d)−D(m−1,dPm)
那么,通过D数组与C数组交替递推。并且选择一个恰当的分界
B
我们计算D[][1...B] 和
C[][1...⌊nB+1⌋]
有:
dPm≤B:D[m][d]=D[m−1][d]−D[m−1][dPm]dPm>B:D[m][d]=D[m−1][d]−C[m−1][⌊ndPm⌋]
<script type="math/tex; mode=display" id="MathJax-Element-44"></script>
C[m][j]=C[m−1][j]−C[m][⌊jPm⌋]
在
B=⌊n‾‾√⌋
时,这样计算的复杂度:
O(nlogn)
并不比前面说到的方法快多少。
优化1:
通过上面的容斥原理很容易有:
<script type="math/tex; mode=display" id="MathJax-Element-48"></script>
当
Pm+1>n
时:
C(m,n)=1
则:
P2m>n
时:
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)=C(m−1,n)−1
优化2:
当
P2m≥⌊(nB+1)14⌋
时:
C(m,j)=π(max(Pm,j))−m+1
综上 。我们分两段递推:
第一段,前:
Pm≤⌊(nB+1)14⌋
。正常计算
第二段。此时
P2m>⌊(nB+1)14⌋
对于第二个优化。预处理出前
⌊n‾√⌋
个
π[]
.
此时不在更新
C[][]
,
如果需要C数组的信息。利用:优化2的式子得到。
对于D数组,应用优化1:
P2m>⌊nd⌋
时,即:
d>⌊nP2m⌋
D(m,d)=D(m−1,d)−1
可以肯定的是。不对
d∈[⌊nP2m⌋+1,B]
更新
仅仅维护
d∈[1,⌊nP2m⌋]
时。
我们记录最早开始不更新的哪个素数标号。并预处理前缀和。必要时刻查前缀和的表。即可。对于没有记录的区间有需要更新的区间。我们暴力更新。时间复杂度不会增加。(可以自行证明。很简单)
还有一种笨的方法。也不会很慢
笨的方法就是:使用区间数据结构来维护(总复杂度多了个 log?)。
建议使用数状数组。(常数小。也就慢了700ms。。。。
所以第二阶段。我们仅仅计算了
d∈[1,⌊nP2m⌋]
那么总多时间复杂度:
T(n)=∑P≤(nB+1)1/4(O(B)+O(nB))+∑(nB+1)1/4<P≤n1/2O(np2)
在
B=n‾‾√
时:
T(n)=O(n34logn)
但是中间维护咱们不是
O(1)
。
使用数据结构的话是:
T(n)=O(n34)
注意:修改对C的定义。让其变为[2,n]上的数并不会优化计算。这是因为中间的递推依然会多出来一个余项。
下面是代码。。。(可以快速计算100亿以内的答案。)
#include <algorithm>
#include <stdio.h>
#include <string.h>
#include <cmath>
#define MAXN 1111111
using namespace std;
typedef long long LL;
struct arry
{
int A[MAXN];
int n;
arry()
{
memset(A,0,sizeof A);
}
void clear(int m)
{
memset(A,0,(m+1)*sizeof(int));
n=m;
}
int lowbit(int x)
{
return x&(-x);
}
void add(int x,int key)
{
while(x<=n)
{
A[x]+=key;
x+=lowbit(x);
}
}
int sum(int x)
{
int ans=0;
while(x)
{
ans+=A[x];
x-=lowbit(x);
}
return ans;
}
}S;
int prim[MAXN],deep=1;
int pi[MAXN];
LL G[MAXN];
LL C[MAXN];
void init()
{
for(int i=2;i<MAXN;i++)pi[i]=1;
for(int i=2;i<MAXN;i++)
{
if(!pi[i])continue;
prim[deep++]=i;
for(int j=i<<1;j<MAXN;j+=i) pi[j]=0;
}
for(int i=2;i<MAXN;i++) pi[i]+=pi[i-1];
}
void clat_1(int m,int k,LL n)
{
LL p=prim[k];
for(int i=1;i<m;i++)
{
LL d=i*p;
LL u=n/d;
if(u<m)
G[i]-=C[u];
else
G[i]-=G[d];
}
for(int i=m;i;i--)C[i]-=C[i/p];
}
LL slove(LL n)
{
if(n<MAXN)return pi[n];
int m=sqrt(n)+1.1;
int n_4=pow(n,1.0/4.0)+1.1;
for(int i=1;i<=m;i++)
{
C[i]=i;
G[i]=n/i;
}
int k;
for(k=1;prim[k]<n_4;k++)clat_1(m,k,n);
S.clear(m);
while(prim[k]<m)
{
LL p=(LL)prim[k]*prim[k];
LL lim=n/p+1;
for(int d=1;d<lim;d++)
{
LL u=(LL)d*prim[k];
LL b=n/u;
if(b<m)
{
if(b<=prim[k-1])
G[d]-=1;
else
G[d]-=pi[b]-k+2;
}
else
G[d]-=G[u]-S.sum((int)u);
}
S.add((int)lim,1);
k++;
}
return k+G[1]-2;
}
int main ()
{
init();
LL n;
while(scanf("%lld",&n)==1) printf("%lld\n",slove(n));
return 0;
}
增加内容(与维基百科上的一致):
上面的递推有一个简化和推广。(其实是针对上面递推变形出来的一个)
记:
Gk(i,j)
表示:
[1,j]
上,由
k
个大于Pi 的素数组成的数的数量。
那么有
C(i,j)=∑k=0∞Gk(i,j)
那么:
C(π(j13),j)=∑k=0∞Gk(π(j13),j)=∑k=02Gk(π(j13),j)+∑k=3∞Gk(π(j13),j)
明显1:
∑k=3∞Gk(π(j13),j)=0
明显2:
G0(i,j)=1G1(i,j)=π(max(j,Pi))−i
<script type="math/tex; mode=display" id="MathJax-Element-79"></script>
对于
G2(π(j13),j)
.它的可选素数是有
P>j13
那么对于另一个素数
P′>j13
,且
PP′≤j , P≤P′
.
这样的素数
P′
个数为:
π(jP)−π(P)+1
所以:
G2(π(j13),j)=∑j1/3<P≤j(π(jP)−π(P)+1)
综上:
π(j)=C(π(j13),j)−G0(π(j13),j)−G2(π(j13),j)+π(j13)
这种方法 。感觉前一部还是要打表出
n1/4
比较好。我比较弱。还是不知道优化的优雅方法。
这种方法空间复杂度比较高。有些人貌似是部分记忆话。效率还很高呢。
#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
typedef long long ll;
#define INF 30000000000
#define chk(pos, i) (((pos[i/64])&(1<<((i>>1)&31))))
#define set(pos, i) (((pos[i/64])|=(1<<((i>>1)&31))))
#define check(x) ((x&&(x&1)&&(!chk(pos, x)))||(x==2))
const int N = 101;
const int M = 49500;
const int P = 700000;
const int UP = 5000000;
ll tmp[N][M];
unsigned int pos[157000] = {0};
int len = 0, pm[P], cnt[UP];
void init(){
set(pos, 0), set(pos, 1);
for(ll i = 3; i*i < UP; i += 2){
if(!chk(pos, i)){
ll k = i<<1;
for (ll j = (i * i); j < UP; j += k) set(pos, j);
}
}
for(int i = 1; i < UP; ++i){
cnt[i] = cnt[i-1];
if(check(i)) pm[len++] = i, cnt[i]++;
}
for(ll n = 0; n < N; ++n){
for(int m = 0; m < M; ++m){
if(!n){ tmp[n][m] = m; continue; }
tmp[n][m] = tmp[n-1][m]-tmp[n-1][m/pm[n - 1]];
}
}
}
ll euler(ll m, int n){
if(n == 0) return m;
if(pm[n - 1] >= m) return 1;
if(m < M && n < N) return tmp[n][m];
return euler(m, n - 1) - euler(m / pm[n - 1], n - 1);
}
ll solve(ll m){
if(m < UP) return cnt[m];
int s = sqrt(m+0.9);
int y = pow(m+0.9,1.0/3.0);
ll res = euler(m, cnt[y])+cnt[y]-1;
for (int i = cnt[y]; pm[i] <= s; i++){
res += - solve(m/pm[i]) + solve(pm[i]) - 1;
}
return res;
}
int main() {
init();
ll n;
while(scanf("%lld",&n)==1)printf("%lld\n",solve(n));
return 0;
}
则:
C(π(n‾‾√),n)=π(n)−π(n‾‾√)+1
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)
kPm, 其中 k∈[1,⌊nPm⌋]
C(m−1,⌊nPm⌋)
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)
C(0,n)=⌊n⌋
D(m,d)=C(m,⌊nd⌋)
⌊⌊na⌋b⌋=⌊nab⌋
D(m,d)=D(m−1,d)−D(m−1,dPm)
dPm≤B:D[m][d]=D[m−1][d]−D[m−1][dPm]dPm>B:D[m][d]=D[m−1][d]−C[m−1][⌊ndPm⌋]
C[m][j]=C[m−1][j]−C[m][⌊jPm⌋]
C(m,n)=1
C(m,n)=C(m−1,n)−C(m−1,⌊nPm⌋)=C(m−1,n)−1
C(m,j)=π(max(Pm,j))−m+1
D(m,d)=D(m−1,d)−1
T(n)=∑P≤(nB+1)1/4(O(B)+O(nB))+∑(nB+1)1/4<P≤n1/2O(np2)
#include <algorithm>
#include <stdio.h>
#include <string.h>
#include <cmath>
#define MAXN 1111111
using namespace std;
typedef long long LL;
struct arry
{
int A[MAXN];
int n;
arry()
{
memset(A,0,sizeof A);
}
void clear(int m)
{
memset(A,0,(m+1)*sizeof(int));
n=m;
}
int lowbit(int x)
{
return x&(-x);
}
void add(int x,int key)
{
while(x<=n)
{
A[x]+=key;
x+=lowbit(x);
}
}
int sum(int x)
{
int ans=0;
while(x)
{
ans+=A[x];
x-=lowbit(x);
}
return ans;
}
}S;
int prim[MAXN],deep=1;
int pi[MAXN];
LL G[MAXN];
LL C[MAXN];
void init()
{
for(int i=2;i<MAXN;i++)pi[i]=1;
for(int i=2;i<MAXN;i++)
{
if(!pi[i])continue;
prim[deep++]=i;
for(int j=i<<1;j<MAXN;j+=i) pi[j]=0;
}
for(int i=2;i<MAXN;i++) pi[i]+=pi[i-1];
}
void clat_1(int m,int k,LL n)
{
LL p=prim[k];
for(int i=1;i<m;i++)
{
LL d=i*p;
LL u=n/d;
if(u<m)
G[i]-=C[u];
else
G[i]-=G[d];
}
for(int i=m;i;i--)C[i]-=C[i/p];
}
LL slove(LL n)
{
if(n<MAXN)return pi[n];
int m=sqrt(n)+1.1;
int n_4=pow(n,1.0/4.0)+1.1;
for(int i=1;i<=m;i++)
{
C[i]=i;
G[i]=n/i;
}
int k;
for(k=1;prim[k]<n_4;k++)clat_1(m,k,n);
S.clear(m);
while(prim[k]<m)
{
LL p=(LL)prim[k]*prim[k];
LL lim=n/p+1;
for(int d=1;d<lim;d++)
{
LL u=(LL)d*prim[k];
LL b=n/u;
if(b<m)
{
if(b<=prim[k-1])
G[d]-=1;
else
G[d]-=pi[b]-k+2;
}
else
G[d]-=G[u]-S.sum((int)u);
}
S.add((int)lim,1);
k++;
}
return k+G[1]-2;
}
int main ()
{
init();
LL n;
while(scanf("%lld",&n)==1) printf("%lld\n",slove(n));
return 0;
}
C(i,j)=∑k=0∞Gk(i,j)
C(π(j13),j)=∑k=0∞Gk(π(j13),j)=∑k=02Gk(π(j13),j)+∑k=3∞Gk(π(j13),j)
∑k=3∞Gk(π(j13),j)=0
G0(i,j)=1G1(i,j)=π(max(j,Pi))−i
π(jP)−π(P)+1
G2(π(j13),j)=∑j1/3<P≤j(π(jP)−π(P)+1)
π(j)=C(π(j13),j)−G0(π(j13),j)−G2(π(j13),j)+π(j13)
#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
using namespace std;
typedef long long ll;
#define INF 30000000000
#define chk(pos, i) (((pos[i/64])&(1<<((i>>1)&31))))
#define set(pos, i) (((pos[i/64])|=(1<<((i>>1)&31))))
#define check(x) ((x&&(x&1)&&(!chk(pos, x)))||(x==2))
const int N = 101;
const int M = 49500;
const int P = 700000;
const int UP = 5000000;
ll tmp[N][M];
unsigned int pos[157000] = {0};
int len = 0, pm[P], cnt[UP];
void init(){
set(pos, 0), set(pos, 1);
for(ll i = 3; i*i < UP; i += 2){
if(!chk(pos, i)){
ll k = i<<1;
for (ll j = (i * i); j < UP; j += k) set(pos, j);
}
}
for(int i = 1; i < UP; ++i){
cnt[i] = cnt[i-1];
if(check(i)) pm[len++] = i, cnt[i]++;
}
for(ll n = 0; n < N; ++n){
for(int m = 0; m < M; ++m){
if(!n){ tmp[n][m] = m; continue; }
tmp[n][m] = tmp[n-1][m]-tmp[n-1][m/pm[n - 1]];
}
}
}
ll euler(ll m, int n){
if(n == 0) return m;
if(pm[n - 1] >= m) return 1;
if(m < M && n < N) return tmp[n][m];
return euler(m, n - 1) - euler(m / pm[n - 1], n - 1);
}
ll solve(ll m){
if(m < UP) return cnt[m];
int s = sqrt(m+0.9);
int y = pow(m+0.9,1.0/3.0);
ll res = euler(m, cnt[y])+cnt[y]-1;
for (int i = cnt[y]; pm[i] <= s; i++){
res += - solve(m/pm[i]) + solve(pm[i]) - 1;
}
return res;
}
int main() {
init();
ll n;
while(scanf("%lld",&n)==1)printf("%lld\n",solve(n));
return 0;
}