Prime number The Ulam spiral

Question

I want a simple method of visualizing the prime numbers that reveals the apparent tendency of certain quadratic polynomials to generate unusually large numbers of primes.

I was able to display the prime numbers using this

Table[Prime[n], {n, 41538}] 

I am new to Mathematica and I would like a guide on how I can do this in mathematica. I am working on something similar and i wanted to see some of the different done ulam spiral.. is there is any. Thanks in advance.

Answer 1

The fastest way is probably a compiled sieve of Eratosthenes.

PrimesUpTo = Compile[{{n, _Integer}},
  Block[{S = Range[2, n]},
    Do[
      If[S[[i]] != 0,
        Do[
          S[[k]] = 0,
          {k, 2i+1, n-1, i+1}
        ]
      ],
      {i, Sqrt[n]}
    ];
    Select[S, Positive]
  ],
  CompilationTarget -> "C",
  RuntimeOptions -> "Speed"
];

It's fast for your input:

PrimesUpTo[500000] // Length // AbsoluteTiming

{0.009409, 41538}

It's fast(ish) for larger inputs too, but if you go any higher than the following examples, you'll probably want to use a segmented sieve.

PrimesUpTo[10^7] // Length // AbsoluteTiming

{0.337327, 664579}

PrimesUpTo[10^8] // Length // AbsoluteTiming

{3.546054, 5761455}

To find all primes less than $n$ in a segmented fashion, you'll need to first find all primes $\leq \sqrt n$.

SegmentedPrimeSieve = Compile[{{l, _Integer}, {hi, _Integer}, {primes, _Integer, 1}},
  Module[{lo = Max[l, 2], S, sqrt = Round[Sqrt[hi]], p, st},
    S = Range[lo, hi];
      Do[
        p = primes[[i]];

        If[p > sqrt, Break[]];
        st = lo + Mod[-lo, p];

        If[st > hi, Continue[]];
        If[st == p, st += p];
        Do[
          S[[k]] = 0,
          {k, st - lo + 1, hi - lo + 1, p}
        ],
        {i, 1, Length[primes]}
      ];

      Select[S, Positive]
    ],
    CompilationTarget -> "C",
    RuntimeOptions -> "Speed"
];

Now let's count all primes less than 1 billion in batches of length 10^4. This takes ~17.5 sec on my computer.

n = 10^9;
ps = PrimesUpTo[Ceiling[Sqrt[n]]];
Δ = 10^4;

Total @ Table[
  Length[SegmentedPrimeSieve[1 + i*Δ, (i + 1)Δ, ps]], 
  {i, 0, n/Δ - 1}] // AbsoluteTiming

{17.582762, 50847534}

If you have a lot of cores you can speed things up with ParallelTable. With 12 cores, I get about a 7.25 times speedup.

LaunchKernels[];
DistributeDefinitions[SegmentedPrimeSieve];

Total @ ParallelTable[
  Length[SegmentedPrimeSieve[1 + i*Δ, (i + 1)Δ, ps]], 
  {i, 0, n/Δ - 1}] // AbsoluteTiming

{2.413624, 50847534}

Let's try 10 billion in batches of length 10^5.

n = 10^10;
ps = PrimesUpTo[Ceiling[Sqrt[n]]];
Δ = 10^5;

Total @ ParallelTable[
  Length[SegmentedPrimeSieve[1 + i*Δ, (i + 1)Δ, ps]], 
  {i, 0, n/Δ - 1}] // AbsoluteTiming

{30.046231, 455052511}

The only two ways I can think of to get more substantial speedups are

• Have factors of 2, 3, and 5 pre-sieved (aka wheel sieve). This will reduce the amount of work by 70%!
• Low level type optimizations you can't really do in Mathematica would provide substantial speed ups :(.

The fastest prime generator I know of is primesieve, and it employs the above ideas. On my machine it's ~60 times faster to generate all primes < 10^9 and ~71.5 times faster to generate all primes < 10^10... Impressive!