一般化してみました。
√Ｎとして、１０１までの素数を準備しているので、
１から１００００までの範囲で素数の個数を求めることができます。
この計算式では、Ｎが大きくなると組合せの数が大きくなり、実用的ではありません。
PUBLIC NUMERIC N
LET N=1000
!１００までの素数
DATA 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101
PUBLIC NUMERIC P(26)
MAT READ P
PUBLIC NUMERIC W
LET W=0
LET R=0 !√Nまでの素数の個数
DO WHILE P(R+1)^2<=N
LET R=R+1
LOOP
PRINT "R=";R !debug
DIM A(R)
FOR d=1 TO R !(-1)^d*Σ{1≦i＜j＜k＜ … ≦r}[N/(Pi*Pj*Pk* … )]
CALL IntegerSolution(R,d,1, A)
NEXT d
PRINT "範囲[ 1,";N;"]で素数は"; N +W +R-1; "個"
END
EXTERNAL SUB IntegerSolution(m,d,s, A()) !1≦A1＜A2＜A3＜…＜Ad≦mを満たす整数解
IF s>1 THEN LET t=A(s-1)+1 ELSE LET t=1
FOR i=t TO m
LET A(s)=i
IF d=1 THEN !ひとつの組が求まったら
LET t=P(A(1)) !除数を求める
FOR k=2 TO s
LET t=t*P(A(k))
IF t>N THEN EXIT FOR !Ｎより大きいなら、INT(N/t)は０となる
NEXT k
IF t<=N THEN LET W=W+INT(N/t)*(-1)^s
ELSE
CALL IntegerSolution(m,d-1,s+1, A) !次へ
END IF
NEXT i
END SUB