动态规划公式来自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);
}