动态规划公式来自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成立的数

求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 < 质数表中最大质数的平方。

运行截图：

代码如下：

#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);
}