Knowing the time limit helps determine what kind of effort is needed to make the code fast enough:
- simple trial division
- trial division with wheeled striding through candidates and divisors
- trial division with memoised divisor primes and wheeled striding through candidates
- Sieve of Eratosthenes
- odds-only Sieve of Eratosthenes
- segmented odds-only Sieve of Eratosthenes
- segmented odds-only Sieve of Eratosthenes with presieving
Each complication makes the code significantly faster, but also more complex. So how fast is fast enough?
Generally speaking, the primes up to any given k should only be computed once total instead of twice for each query. Also, given the fact that the query count limit is much higher than the piddling small numbers of primes involved (10000) it makes sense to compute all of them in one go, which is much more efficient and much less complicated than writing some logic that can extend the current list of prime-indexed primes up to some number k when it is too short.
The high number on queries also makes it imperative to benchmark the code for parsing inputs and producing outputs, since that can easily dominate the total program runtime. SPOJ's INOUTEST - and the same or similar tasks at other sites - can be a good guide there.
I've had several such tasks where my submissions clocked times of several seconds or even TLE, even though the raw problem-solving part took only a few milliseconds to run. Writing a dedicated input parser and (batched) output converter cut runtime by two orders of magnitude or more.
Based on my experience recently solving problems of similar structure and magnitude (SPOJ's TDPRIMES and TDKPRIME) I would recommend using a segmented odds-only Sieve of Eratosthenes for building the list of prime-indexed primes in one go. That gives you the best performance for the least amount of complications. It would be possible to speed up your trial division by several orders of magnitude but then it would be more complex than the sieving solution, and the sieve would still blow it out of the water despite using simpler code.
I tried quite a few approaches and strategies, and perhaps my experiences can help you avoid going down some blind alleys.
Recommendation: use an odds-only sieve with bool array that is the same size as your CPU's L1 cache, e.g. 32 KiByte.
Using packed bits in a byte array would allow sieving eight times as much per cache-sized window, and extraction of primes from the sieve would be much more efficient. However, the overhead of indexing (sieve[j>>3] |= 1 << (j&7) instead of a simple sieve[j] = true) means that the sieving using byte[] took slightly longer overall in my tests than sieving plus extraction using bool[], even though the bool code had to sieve eight times as many windows. I was able to make byte[] 30% faster than bool[] but the effort required was essentially the same as that for a modulo 30 wheel, tripling the total size of the source code. Using the stock bitsets of C# (BitArray) and C++ (std::vector<bool>) was an unmitigated performance disaster; the whole program took several times as long to run just because of their mind-boggling inefficiency.
Recommendation: use functions that focus on a specific part of the job and that do whole batches at a time with maximum efficiency.
For example, I found it much faster to extract all primes from a sieved window and then pick out every hundredth from that, instead of integrating the 'take only every 100th prime' bit into the code that extracts primes from the sieve. The reason is focus and simplified looping/branching structures. Modern CPU don't like unpredictable branches at all and levy heavy taxes on them in the form of wait cycles. Better factoring also makes it easier to try different approaches for different sub tasks, for example via function pointers or delegates that can be switched at runtime. And it makes it easier to compose new, completely different solutions from the existing bits and pieces.
For example, instead of complicating - and thus slowing down - the sieve code by the logic for memoising the small handful of primes needed for sieving, I used a separate function for that purpose. Here is the pseudocode (of a hashish persuasion):
internal static List<ushort> small_odd_primes_up_to (int n)
{
Trace.Assert(3 <= n && n <= ushort.MaxValue + 1);
// the '- 1' comes straight from the reversing the conversion formula n = (b << 1) + 1;
// it ensures that even n get treated the same as the immediately preceding odd number
// instead of tagging an unnecessary bit onto the end of the array
int max_bit = (n - 1) >> 1;
// the '- 1' comes from the number-to-bit conversion again and is not related to the sqrt
int sqrt_n_bit = (int)(Math.Sqrt(n) - 1) >> 1;
var sieve = new bool[max_bit + 1];
for (int i = 3 >> 1; i <= sqrt_n_bit; ++i)
if (!sieve[i])
for (int p = (i << 1) + 1, j = p * p >> 1; j <= max_bit; j += p)
sieve[j] = true;
// 3 needs to be handled separately because of the mod 3 stepping in the for loop
var result = new List<ushort>() { 3 };
for (int i = 5 >> 1, d = 1; i <= max_bit; i += d, d ^= 3)
if (!sieve[i])
result.Add((ushort)((i << 1) + 1));
return result;
}
In your case it would have to deliver only the 1168 odd sieve primes, i.e. those not greater than sqrt(1366661) where 1366661 is the 104729th prime and the 10000th prime-indexed prime. This means that the optimisations it contains - odds-only sieving, plus skipping multiples of 3 when reading out the sieved primes - are not really necessary for the problem at hand. However, I've used that function unchanged in almost all code challenges related to primes, and there's no reason to fiddle with tried and trusted code. Basically, it's about the simplest possible code that can do the job, with only a few sleight-of-hand optimisations. I've shown it to provide a contrast to the design of the main sieving logic, which adds quite a few complications.
For example, here's the hashish pseudocode of a function that gets handed a sieved window by the segmented sieve function, extracts all primes from it into a cache, and then adds every Kth prime from the cache to the list of results (which is a property of the class):
// phase and result[] are static class fields
internal static int[] cache = new int[MAX_PRIMES_PER_SIEVE];
internal static void sieve_to_array_post_filtered (bool[] sieve, int window_base, int window_bits)
{
// get things into step for the mod 3 loop stepping on top of the mod 2 wheel
// (the loop is anchored at residue 1)
int i = 0, k = 0;
switch (window_base % 3)
{
case 0:
{
i = 1; goto case 2;
}
case 2:
{
if (!sieve[i])
if (i < window_bits)
cache[k++] = window_base + (i << 1);
++i; break;
}
}
for (int e = window_bits - 5; i < e; i += 6)
{
if (!sieve[i])
cache[k++] = window_base + (i << 1);
if (!sieve[i + 2])
cache[k++] = window_base + (i << 1) + 4;
if (!sieve[i + 3])
cache[k++] = window_base + (i << 1) + 6;
if (!sieve[i + 5])
cache[k++] = window_base + (i << 1) + 10;
}
for ( ; i < window_bits; ++i)
if (!sieve[i])
cache[k++] = window_base + (i << 1);
for (i = K - phase - 1; i < k; i += K)
result[result_size++] = cache[i];
phase = (phase + k) % K;
}
Most of this does the exactly same job as the 2nd for loop in the earlier function, but much faster. The main loop uses mod 3 stepping on top of the odds-only sieve, twice unrolled. The loop after that cleans up leftovers, and the final loop extracts every Kth prime, with phase coordinating the correct counting from one sieve window to the next. This is from the TDPRIMES task where every 100th prime was to be printed; adapting the logic to the requirements of finding prime-indexed primes is trivial.
This example shows how the tiny 2nd loop of the first function got blown up quite lot by the tricks required for achieving performance goals. The code is significantly faster than all other variations I tried, and I tried an awful lot of them (because I'm currently learning the hashish language). The speed is made possible by batched operation, focussing on a single job, and reusing buffers; the loop for extracting the Kth prime was essentially tacked on to the end of the function without disturbing anything. That was required because the sieving function simply invokes this via a process_sieve delegate that is supposed to do whatever is necessary with the freshly sieved primes.
To complete the picture, here's the pseudo code of the segmented sieving logic. It's essentially a normal Sieve of Eratosthenes that sieves the target range in cache-sized increments, remembering the working offsets for each prime from one iteration to the next and pumping the sieved primes out via the process_sieve delegate to make place for a new batch.
internal static void process ()
{
const int N = 100000000;
var small_odd_primes = small_odd_primes_up_to((int)Math.Sqrt(N)).ToArray();
var offset = new int[small_odd_primes.Length];
// ... some preparatory stuff elided (e.g. setting presieve_pattern and presieve_modulus) ...
int window_base = 3, n = N - (~N & 1), last_prime_index = small_odd_primes.Length - 1;
for (int i = first_sieve_prime_index; i <= last_prime_index; ++i)
{
int prime = small_odd_primes[i], start = prime * prime;
offset[i] = (start - window_base) >> 1;
}
var sieve = new bool[1 << 15];
for (int bits_to_sieve = ((n - window_base) >> 1) + 1; bits_to_sieve > 0; )
{
int window_bits = Math.Min(bits_to_sieve, sieve.Length);
Buffer.BlockCopy(presieve_pattern, (window_base >> 1) % presieve_modulus, sieve, 0, window_bits);
for (int i = first_sieve_prime_index; i <= last_prime_index; ++i)
{
int p = small_odd_primes[i], j = offset[i];
for ( ; j < window_bits; j += p)
sieve[j] = true;
offset[i] = j - window_bits;
}
process_sieve(sieve, window_base, window_bits);
window_base += window_bits << 1;
bits_to_sieve -= window_bits;
}
}
Obviously you'd use N = 1366661 instead of 10^8, and the logic for copying the presieved pattern would be replaced by a simple call to Array.Clear. However, I left the presieve logic in place to show where it fits into the overall scheme, because this simple trick can add quite a performance boost. The reason is that the smallest few primes have the shortest strides and thus have to make many more hops in any window than the rest taken together, and blasting a suitable bit pattern over the sieve is much faster. Without presieving you'd use 0 as first_sieve_prime_index, naturally.
In any case this should make clear that the logic of a segmented sieve can be very simple indeed.
Using a delegate like process_sieve instead of fixed-function code has the added advantage that it simplifies unit testing. I.e. write a delegate that extracts and stores all primes, which can then be checked against any other prime source.
To put things into perspective: my SPOJ TDPRIMES submission from which the pseudo code was taken (printing every 100th prime up to 100,000,000) clocked 190 milliseconds on SPOJ (C#). The time for sieving the first 10,000 prime-indexed primes with this approach should be on the order of a millisecond.
