Create a 1xN vector containing the first N primes?

[Alex Stead]

Hi all,

This may be fairly simple, but I can't figure it out. Say I wanted to create a vector containing the first N primes, how could I write some code to do that? I am using version 13.1 if that's relevant. Thanks,

Alex

[Nick Cox]

If this were my problem I would take a program from another language I knew and translate into Stata. I don't know what other languages you know. Depending how big you want N to be -- million? billion? -- I wouldn't want to write one from scratch.

[Scott Merryman]

You might find Stanislav Kolenikov's -primes- as a useful start. Or use Mata (taken from http://www.econometricsbysimulation....a-as-well.html):

mata:
real matrix function primes(real scalar n) {
prime = 1
     set = 2::n
     while (length(set) > 0) {
          prime = prime \ set[1]
          set = select(set, floor(set / set[1]) :!= set /set[1])
     }
return(prime)
}
end

mata : primes(1000)

[Mike Lacy]

The concision of Scott's code really impressed me (not that I'm sure I understand it!), as I had recently translated to Mata some code I found somewhere for the Sieve of Eratosthenes. My sieve code was much longer (see below). For amusement, I timed the two approaches. For N = 10,000, the sieve code ran in about 15% of the time. For N = 1,000,000, the sieve code ran in less than 1% of the time! I doubt that time matters for the current application, but this nicely illustrates that a concise Mata code might not be fast. (That's contrary to my usual experience.)

// Sieve of Eratosthenes
mata
real matrix function sieve(real scalar N)  {
   //
    A= J(N,1,1)
    for ( i = 2; i <= floor(sqrt(N)); i++) {
      if (A[i,1] == 1) {
          j = i^2
          while (j <= N) {
              A[j,1] = 0
              j = j + i
          }
      }
    }
   // The preceding algorithm creates a matrix with A[i,1] == 1 ==> i is prime.
  // Use A to create a list of the primes. There's probably a better way to do this.
    pos = 1
    primelist = J(sum(A), 1, -1)
    for ( i = 1; i <= rows(A); i++) {
        if (A[i,1] ==1) {
            primelist[pos,1] = i
            ++pos
        }
    }
    return(primelist)
}
end
//
timer clear 1
timer on 1
forval i = 1/100 {
    mata : m= sieve(10000)
}
timer off 1
timer list 1