J Primes Enumeration

Question

J will answer the n-th prime via p:n.

If I ask for the 100 millionth prime I get an almost instant answer. I cannot imagine J is sieving for that prime that quickly, but neither looking it up in a table as that table would be around 1GB in size.

There are equations giving approximations to the number of primes to a bound, but they are only approximations.

How is J finding the answer so quickly ?

Corody

I'm no expert, but looking at v2.c, it looks like J has a lot of big prime tables built into it, and maybe sieves an additional bunch of them on start-up. Also, I bet on a modern computer, with a clever algorithm that makes a lot of "space-for-time sacrifices", finding the 100 millionth prime takes a lot less time than you think! – Octavla

Some insight can be gained from this 2001 forum post by Roger Hui. I would guess that more advanced techniques have been put into play in the last 14 years. jsoftware.com/pipermail/general/2001-August/007334.html – Schutzstaffel

Answer 1

J uses a table to start, then calculates

NOTE! This is speculation, based on benchmarks (shown below).

If you want to quickly try for yourself, try the following:

p:1e8   NB. near-instant
p:1e8-1 NB. noticeable pause

The low points on the graph are where J looks up the prime in a table. After that, J is calculating the value from a particular starting point so it doesn't have to calculate the entire thing. So some lucky primes will be constant time (simple table lookup) but generally there's first a table lookup, and then a calculation. But happily, it calculates starting from the previous table lookup instead of calculating the entire value.

Benchmarks

I did some benchmarking to see how p: performs on my machine (iMac i5, 16G RAM). I'm using J803. The results are interesting. I'm guessing the sawtooth pattern in the time plots (visible on the 'up to 2e5' plot) is lookup table related, while the overall log-ish shape (visible on the 'up to 1e7' plot) is CPU related.

NB. my test script
ts=:3 : 0
  a=.y
  while. a do.
  c=.timespacex 'p:(1e4*a)'  NB. 1000 times a
  a=.<:a
  b=.c;b
  end.
  }:b
)

a =: ts 200

require'plot'
plot >0{each a  NB. time
plot >1{each a  NB. space

(p: up to 2e5)

time

space

(p: up to 1e7)

time

space

During these runs one core was hovering around 100%:

Also, the voc page states:

Currently, arguments larger than 2^31 are tested to be prime according to a probabilistic algorithm (Miller-Rabin).

And in addition to a prime lookup table as @Mauris points out, v2.c contains this function:

static F1(jtdetmr){A z;B*zv;I d,h,i,n,wn,*wv;
 RZ(w=vi(w));
 wn=AN(w); wv=AV(w);
 GA(z,B01,wn,AR(w),AS(w)); zv=BAV(z);
 for(i=0;i<wn;++i){
  n=*wv++;
  if(1>=n||!(1&n)||0==n%3||0==n%5){*zv++=0; continue;}
  h=0; d=n-1; while(!(1&d)){++h; d>>=1;}
  if     (n< 9080191)*zv++=spspd(31,n,d,h)&&spspd(73,n,d,h);
  else if(n<94906266)*zv++=spspd(2 ,n,d,h)&&spspd( 7,n,d,h)&&spspd(61,n,d,h);
  else               *zv++=spspx(2 ,n,d,h)&&spspx( 7,n,d,h)&&spspx(61,n,d,h);
 }
 RE(0); R z;
}    /* deterministic Miller-Rabin */ 

Whirlybird