动态规划求小于等于n的质数个数

时间:2021-12-02 13:11:59

动态规划公式来自Codeforces 665F 这题的官方题解,没有看懂题解中是怎么计算的,于是自己写了个计算的程序。


首先,P[ t ] 表示第 t 个质数,P[1]=2 , P[2]=3, P[3]=5 等

状态: D[ n , m] 表示 1到n这n个数中,有多少个数的所有质因子都大于等于P[m]。

边界条件:D[n,1]=n    ,  D[1,m]=1。

(1没有质因子,但是每个D[n,m]中都要算上1,不然会导致递推式出错)。


递推公式 :D[n,m+1] = D[ n , m ] - D[ n/P[m] , m ]   其中的除法是整数除法(即对除法结果下取整)

翻译一下:【1到n中】质因数都大于等于P[m+1]的数的数量 = 【1到n中】质因数都大于等于P[m]的数的数量 

-     【1到n中】最小质因数为P[m]的数的数量

其中,【1到n中】最小质因数为P[m]的数的数量   = 【1到(n/P[m])中】质因数都大于等于P[m]的数的数量

于是就有了递推式。


有了上述定义之后, 

对于n,求最小的m使得P[m]^2 > n

那么,小于等于n的质数数量 = (小于P[m]的质数数量)+(大于等于P[m]又小于等于n的质数数量)

=(m-1)+( D[n,m] - 1)

= D[n,m] + m - 2


小于P[m]的质数数量为m-1因为P[m]是第m个质数。

大于等于P[m]的质数数量为D[n,m]-1,因为P[m]^2 > n,所以所有质因子都大于等于P[m]的数都只有一个因数,

也就是说这些都是质数。当然还要减去1.因为每个D中都计算了数字1.

所以小于等于n的质数数量 = D[n,m] + m - 2 。


剩下就是计算D[n,m]了。


首先当P[m]>n时,D[n,m]为1.

当P[m1]<=n 且P[m2]<=n 且 P[m1]^2>n 且 P[m2]^2>n 时,有公式D[n,m1]+m1=D[n,m2]+m2

所以m > k 时,可以直接计算D[n,k],然后D[n,m]=D[n,k]- (m - k)。 其中k是最小的使 P[k]^2 > n成立的数

动态规划求小于等于n的质数个数

求D值的代码:

LL GetD(LL n,int m){
assert(m <= Pn);
//到达边界直接返回
if(P[m] > n || n == 1) return 1;
if(m == 1) return n;

//不是边界就递归
int k = GetK(n);
if(m > k) return GetD(n,k) - (m - k);
else return GetD(n,m-1) - GetD(n/P[m-1],m-1);
}




尝试了很多种优化,但是效果都不明显。直接递归也不算太慢。

计算n之前,要保证 n < 质数表中最大质数的平方。


运行截图:

动态规划求小于等于n的质数个数


代码如下:

#include <iostream>
#include <map>
#include <queue>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#define assert(a) if(!(a)) {fprintf(stderr,#a);fprintf(stderr,"\n");exit(1);}
#define LL long long
using namespace std;
int MxPF[1000000];//最大质因子表
int P[100000],Pn;//P质数表(P[1]=2) Pn:质数个数
LL P2[100000];//P2[i]=P[i]*P[i] 加速计算

void Init();//初始化质数表
LL GetPI(LL n);//计算PI(n)
int GetK(LL n);//计算K(n)
LL GetD(LL n,int m);//递归计算D[n,m]

int main(void){
Init();//打质数表
LL n;
while(cin>>n){
//计算答案,计时
clock_t start = clock();
LL ANS = GetPI(n);
clock_t end = clock();

//输出
cout<<"PI("<<n <<")="<<ANS<<" 用时: "<<(double)(end-start)/CLOCKS_PER_SEC<<" 秒" <<endl;
}
return 0;
}

void Init(){//初始化质数表
memset(MxPF,-1,sizeof(MxPF));
Pn=0;
for(int i=2;i<=120000;++i){
if(~MxPF[i]) continue;
P[++Pn]=i;
P2[Pn]=(LL)P[Pn]*P[Pn];
for(int j=i;j<=120000;j+=i) MxPF[j]=i;
}
}
LL GetPI(LL n){//返回有多少个小于等于n的质数
assert(n>0);
assert(n < P2[Pn]);
int m = GetK(n);
return GetD(n,m) + m - 2;
}
int GetK(LL n){//求最小的k使得 P[k]^2 > n
int L = 1,R = Pn,M;//[L,R] first ^2 > n
while(L ^ R){
M = (L + R) >> 1;
if(P2[M] > n) R = M;
else L = M + 1;
}
return L;
}
LL GetD(LL n,int m){
assert(m <= Pn);
//到达边界直接返回
if(P[m] > n || n == 1) return 1;
if(m == 1) return n;

//不是边界就递归
int k = GetK(n);
if(m > k) return GetD(n,k) - (m - k);
else return GetD(n,m-1) - GetD(n/P[m-1],m-1);
}