Nまでの素数の個数をO(N)で求める方法を思いついたので投稿します。
以下プログラムです。
#include <iostream>
#include <vector>
#define ll long long int
using namespace std;
vector<ll> p;
ll D(ll n,ll k);
int main(void){
ll ans=0,N=0;
int flag=0;
p.push_back(2);
while(1){
while(N<=3){
cout<<"Nまでの素数の個数を求めます。(N>3)"<<endl;
cout<<"N:";
cin>>N;
}
for(int i=p[p.size()-1]+1; p[p.size()-1]*p[p.size()-1]<=N; i++){
for(int j=0; p[j]*p[j]<=i; j++){
if(i%p[j]==0)flag=1;
}
if(flag==0)p.push_back(i);
flag=0;
}
ans = N-1-D(N,p.size()-2);
cout<<N<<"までの素数の個数: "<<ans<<endl;
ans=0;
N=0;
}
return 0;
}
ll D(ll n,ll k){
ll tem,ret;
if(n<=3)return 0;
if(p[k]*p[k]>n) return D(n,k-1);
if(k==0)return (n/2)-1;
tem = n/p[k];
ret = D(n,k-1)+tem-1-k-D(tem,k-1);
return ret;
}
それぞれ説明していきます。
まずこのアルゴリズム自体は「√Nまでの素数が全てわかればNまでの素数の個数が求まる」というものです。どうやって求めているかを大雑把に説明すると、Nまでの合成数は全て√N以下の素数で割り切れるというのを利用して合成数の個数を求め、それをNから引くことで求めています。（厳密には1が残るので１を引かないといけません）
√N以下の素数を求めている場所はここです
for(int i=p[p.size()-1]+1; p[p.size()-1]*p[p.size()-1]<=N; i++){
for(int j=0; p[j]*p[j]<=i; j++){
if(i%p[j]==0)flag=1;
}
if(flag==0)p.push_back(i);
flag=0;
}
pは素数を表すvectorです。（p[0]=2, p[1]=3, p[2]=5, . . .）
初期の状態でp[0]だけ2になってるのでこれで求めることができます。
ここは本題ではないので細かい解説は省きます。
次は本題のD関数です。
D(n,k)はn以下の合成数のうちk番目までの素数で割り切れるものの個数です。
例えばD(10,1)なら10以下の合成数のうち1番目までの素数(2と3)で割り切れるものの個数となります。なので4,6,8,9,10の5つで5となります。
そしてこの関数からは次の漸化式が成り立ちます。
D(n,0) = n/2 - 1
D(n,k) = D(n,k-1) + n/p[k] - 1 - k - D(n/p[k], k-1)
割り算は全てintの除算ですので、切り捨てです。
一つ目の式は簡単に理解できると思います。
単純に2で割れるものの個数がn/2で素数である2を含んでいるので1を引いてます。
次の式は少し難しいです。
k-1番目までの素数で割り切れる合成数の個数にk番目の素数のみで割り切れる合成数の個数を足してます。k番目の素数のみで割り切れる合成数というのは、p[k]で割り切れてp[t]（t<k)で割り切れない合成数のことをさしています。
k-1番目までの素数で割り切れる合成数の個数はD(n,k-1)ですね。定義の通りです。
k番目の素数のみで割り切れる合成数は、まず一つ目の式のようにn/p[k] - 1でp[k]で割り切れる合成数の個数を求め、それから重複した数を引いています。
重複した数はk + D(n/p[k], k-1)で求められます。
重複した数というのは例えばnが30、kが2の場合、10,15,20,25,30の中から2と3で割り切れるものをなくせば良いのでまずこれらを5で割り2,3,4,5,6にします。
5で割った後に2,3で割り切れるなら最初から2,3で割り切れるので2,3,4,5,6の中から2,3で割り切れるものをなくせば良いわけです。そしてそれらの個数はD(6,1)+2となるのでk + D(n/p[k], k-1)と同じ形になっているのがわかると思います。
これによって再帰的に求めることができるというのがわかっていただけたと思います。
計算量については√n>p[k]である事を考えると最大でO(√n√√n√√√n√√√√n. . .)となります。また指数部に注目すると1/2+1/4+1/8+1/16+. . .=1となるのでO(n)になります~~(多分)~~
実際Dが実行される回数を数えるとNを2倍すると2倍くらいになってたのであってると思います。
いずれ時間があればメモ化再帰にしたりkが小さい時の処理の工夫をしたりすることで高速化ができるという事について書くかもしれないです。