This algorithm was published in Algorithmica in 2007. A glimpse is available here with other references to more recent work on this computation: link.springer.com/article/10.1007%2Fs00453-007-0162-8 . It is O(k*a1). It is similar to the Eratosthenes method of finding primes; I have also seen it referred to as a round-robin approach. I will annotate the code.