0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ How would you like us to help? You don't tell us what the pattern is, nor what you want to achieve... Table[Prime[n], {n, 41538}] will generate all the primes up to 500,000. Maybe you can start from there and explain more about what you would like to do? $\endgroup$ – MarcoB Jun 11 '15 at 16:52
  • $\begingroup$ Thanks @MarcoB am actually new .. i didnt know where to start now i have an idea. Thanks for this . I believe once i try this i can come up with what i want to see. $\endgroup$ – Madona Syombua Jun 11 '15 at 16:55
  • $\begingroup$ Alternatively, you can use Prime /@ Range[41538] or since Prime is Listable: Prime@Range[41538] $\endgroup$ – Bob Hanlon Jun 11 '15 at 18:33
  • $\begingroup$ I managed to edit my question, i don't know if i can find help on the edited one? $\endgroup$ – Madona Syombua Jun 14 '15 at 18:15
5
$\begingroup$

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!

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ Chip, how does your function compare to the built in? $\endgroup$ – MarcoB Jun 12 '15 at 1:15
  • 1
    $\begingroup$ @MarcoB To find all primes less than $n$, this sieve runs in $O(n \log\log n)$ time. Prime[Range[PrimePi[n]]] (the built in version) runs in essentially $O(n^{5/3})$ time for larger $n$ because Prime[x] uses a $O(n^{2/3})$ time algorithm (see reference.wolfram.com/language/tutorial/…). For n == 10^8, this sieve is about 4 times faster, and the larger $n$ grows, the faster the sieve will be compared to the built in way. $\endgroup$ – Chip Hurst Jun 12 '15 at 2:12
  • $\begingroup$ Chip, thank you for the details, and for the link: I had not actually come across that page yet, and it made for interesting reading. $\endgroup$ – MarcoB Jun 12 '15 at 4:31

Not the answer you're looking for? Browse other questions tagged or ask your own question.