Sieve of Eratosthenes in q

Continuing my adventures in q, I tried writing the sieve of Eratosthenes in a functional style:

p1:{[n]
 s:{[n;a;d] if[a[d]; a[d*d _ til ceiling n%d]:0b]; a};
 2_ where (n#1b) s[n;;]/ 2+til "i"$sqrt n}
p2:{s:{if[y z;y[z*z _til ceiling x%z]:0b];y};2_where(x#1b)s[x]/2+til"i"$sqrt x}
p3:{2_where(x#1b){if[y z;y[z*z _til ceiling x%z]:0b];y}[x]/2+til"i"$sqrt x}
p4:{a:x#1b;i:1;do["i"$sqrt x;if[a i+:1;a[i*i _til ceiling x%i]:0b]];2_where a}

p1, p2 and p3 are essentially the same while p4 is a loopy code. Among them p4 ran the fastest:

q)\t p3 5000000
3509
q)\t p4 5000000
400

Well, it turned out that I was writing functional code in an inefficient way. On 25 Jun, I was informed of Arthur Whitney’s recursive sieve at http://nsl.com/k/sieves.k. Arthur’s code in q (translation courtesy of Attila) is:

s1:{@[x#1b;y*til each ceiling x%y;:;0b]}
p:{[s;n]$[n<4;enlist 2;r,1_where s[n]r:.z.s[s]ceiling sqrt n]}

His code is really elegant. It works by sieving over a set of recursively computed primes. (It is not quite correct for n = 0, 1 and 2 since it returns 2… but that is easily rectified.) To test my understanding, I derived two optimized sieves that pre-sieves multiples of 2 and 3.

s2:{@[x#01b;y*til each ceiling x%y:1_y;:;0b]}
s3:{@[x#010001b;y*til each ceiling x%y:2_y;:;0b]}
p:{[s;n]$[n<10;where 00110101b til n;r,1_where s[n]r:.z.s[s]ceiling sqrt n]}

These recursive sieves are all faster than p4:

q)\t p4 15000000
1314
q)\t p[s1] 15000000
1218
q)\t p[s2] 15000000
1096
q)\t p[s3] 15000000
999

Is it possible to implement the sieve without loops and recursion? After some thought, I came up with this:

s4:{x,where@[y#1b;0 1,x*til each ceiling y%x;:;0b]}
pp:{[s;n]$[n<2;0#0;()s/reverse{ceiling sqrt x}scan n]}

It uses “reduction” instead of recursion. The performance of pp[s4] and p[s1] are on par.

I also tested q against Python and C++. To generate the primes under 15 million, it took Python about 14845ms and C++ about 488ms. I find q’s performance impressive, considering it is interpreted and not compiled.

The codes above are available for download: soe_q. The Python and C++ code used are:

def sieve(n):
    p = [True]*n
    for i in xrange(2, int(math.sqrt(n))+1):
        if p[i]:
            for j in xrange(i*i, n, i):
                p[j] = False
    return [i for i in xrange(2,n) if p[i]]

vector<int> sieve_bool(int n) {
    vector<int> p(n, true);
    p[0] = p[1] = false;
    for(int i = 2; i < sqrt(n)+1; i++)
        if(p[i]) for(int j = i*i; j < n; j+=i) p[j] = false;
    return p;
}

vector<int> sieve(int n) {
    vector<int> p = sieve_bool(n);
    vector<int> pr;
    for(int i = 0; i < n; ++i)
        if(p[i]) pr.push_back(i);
    return pr;
}