We start with 1 and find primes. This is memory intensive. I can't do more than 4 steps.
Full specs at Mathematica.SE.
primes = {1};
Do[
limv = primes[[-1]]^2 + 2 primes[[-1]];
g = Graph[Flatten[Table[p -> p (m + 1), {p, primes}, {m, 1, limv/p - 1}]]];
primes = Join[primes, Flatten[Position[VertexDegree[g], 1]]];
Print[Length[primes] - 1 , " primes <= ", limv],
{4}]
PrimePi[limv]
2 primes <= 3
6 primes <= 15
44 primes <= 195
3986 primes <= 37635
3986
A slightly different version finds over 58000 primes in 7 seconds. However it requires Prime[] awareness.
