# Primes numbers with prime index

Write a program or function that outputs/returns the first 10000 prime-indexed prime numbers.

If we call the nth prime p(n), this list is

3, 5, 11, 17, 31, 41, 59 ... 1366661


because

p(p(1)) = p(2) = 3
p(p(2)) = p(3) = 5
p(p(3)) = p(5) = 11
p(p(4)) = p(7) = 17
...
p(p(10000)) = p(104729) = 1366661


Standard loopholes are forbidden, and standard output methods are allowed. You may answer with a full program, a named function, or an anonymous function.

• You should generally try to post challenges in the sandbox first (see the link on the right side) in order to work out the issues. – aditsu quit because SE is EVIL Oct 23 '15 at 17:43
• Optimizing for runtime is not what we do in a code-golf challenge; the shortest program always wins. – lirtosiast Oct 23 '15 at 18:53
• Primes with prime subscripts: A006450. – user12166 Oct 23 '15 at 21:55
• @bilbo Answers for code golf are usually accepted after a week, and should be accepted as the shortest successful code. If you wanted code speed, there is a tag for that. See this page about the tag code-golf. – Addison Crump Oct 24 '15 at 12:04
• All contests need an objective winning criterion; they are off topic otherwise. If you're going to judge answers by size and speed, you need to disclose a way to combine both. This should be done when the contest is posted, not 14 hours and 10 answers later. I have undone all speed-related edits, since the only other option would be to close this post for being off topic. – Dennis Oct 25 '15 at 14:10

# MATLAB/Octave, 25 bytes

p=primes(2e6)
p(p(1:1e4))


It doesn't get much more straightforward than this.

## Python, 72 bytes

P=p=1;l=[]
while p<82e5/6:l+=P%p*[p];P*=p*p;p+=1
for x in l:print l[x-1]


This terminates with a "list index out of range error" after printing the 10000 numbers, which is allowed by default.

Uses the Wilson's Theorem method to generate a list l of the primes up to the 10000th prime. Then, prints the primes with the positions in the list, shifted by 1 for zero-indexing, until we run out of bounds after the 10000th prime-th prime.

Conveniently, the upper bound of 1366661 can be estimated as 82e5/6 which is 1366666.6666666667, saving a char.

I'd like a single-loop method, printing prime-indexed primes as we add them, but it seems to be longer.

P=p=1;l=[]
while p<104730:
l+=P%p*[p]
if len(l)in P%p*l:print p
P*=p*p;p+=1

• This is way better than the garbage I was writing. +1 – Mego Oct 23 '15 at 17:51
• This only prints 1229 numbers – aditsu quit because SE is EVIL Oct 23 '15 at 19:03
• @aditsu I think I see my mistake. Are you able to run this code with the bigger bound? – xnor Oct 23 '15 at 19:07
• It will probably take a long time :p – aditsu quit because SE is EVIL Oct 23 '15 at 19:16
• I think it finished ＼（＠；◇；＠）／ , it seems correct – aditsu quit because SE is EVIL Oct 24 '15 at 15:59

## J, 11 bytes

p:<:p:i.1e4


Outputs the primes in the format

3 5 11 17 31 41 59 67 83 109 127 ...


## Explanation

        1e4  Fancy name for 10000
i.     Integers from 0 to 9999
p:       Index into primes: this gives 2 3 5 7 11 ...
<:         Decrement each prime (J arrays are 0-based)
p:           Index into primes again


# Mathematica, 2625 23 bytes

Prime@Prime@Range@1*^4&


Pure function returning the list.

• Prime is Listable so a simple Prime@Prime@Range@1*^4& will do – user46060 Oct 23 '15 at 22:34
• I know the feeling ... In any case, I think this is the prettiest Mathematica solution I have seen on here! – user46060 Oct 23 '15 at 22:36
• Let me guess, the @ operator has higher precedence than ^ when writing Range@10^4? That's classic Mathematica messing up your game of golf. Good trick! – user46060 Oct 23 '15 at 23:14

p=[x|x<-[2..],all((>0).mod x)[2..x-1]]
f=take 10000$map((0:p)!!)p  Outputs: [3,5,11,17,31,41,59,67,83,109,127.....<five hours later>...,1366661] Not very fast. How it works: p is the infinite list of primes (naively checking all mod x ys for y in [2..x-1]) . Take the first 10000 elements of the list you get when 0:p!! (get nth element of p) is mapped on p. I have to adjust the list of primes where I take the elements from by prepending one number (-> 0:), because the index function (!!) is zero based. # PARI/GP, 25 bytes apply(prime,primes(10^4))  # AWK - 129 bytes ...oookay... too long to win points for compactness... but maybe it can gain some honor for the speed? The x file: BEGIN{n=2;i=0;while(n<1366662){if(n in L){p=L[n];del L[n]}else{P[p=n]=++i;if(i in P)print n}j=n+p;while(j in L)j=j+p;L[j]=p;n++}}  Running: $ awk -f x | nl | tail
9991  1365913
9992  1365983
9993  1366019
9994  1366187
9995  1366327
9996  1366433
9997  1366483
9998  1366531
9999  1366609
10000  1366661


BEGIN {
n=2
i=0
while( n<1366662 ) {
if( n in L ) {
p=L[n]
del L[n]
} else {
P[p=n]=++i
if( i in P ) print n
}
j=n+p
while( j in L ) j=j+p
L[j]=p
n++
}
}


The program computes a stream of primes using L as "tape of numbers" holding found primes jumping around on L to flag the nearby numbers already known to have a divisor. These jumping primes will advance while the "tape of numbers" L is chopped off number by number from its beginning.

While chopping off the tape head L[n] being empty means there is no known (prime) divisor.

L[n] holding a value means, this value is a prime and known to divide n.

So either we have found a prime divisor or a new prime. Then ths prime will be advanced to the next L[n+m*p] on the tape found being empty.

This is like the Sieve of Eratosthenes "pulled thru a Klein's bottle". You always act on the tape start. Instead of firing multiples of primes thru the tape, you use the primes being found already as cursors jumping away from the tape start by multiple distances of their own value until a free position is found.

While the outer loop generates one prime or not prime decission per loop, the primes found get counted and stored in P as key, the value of this (key,value) pair is not relevant for the program flow.

If their key i happens to be in P already (i in P), we have a prime of the p(p(i)) breed.

Running:

$time awk -f x.awk | wc -l 10000 real 0m3.675s user 0m3.612s sys 0m0.052s  Take into account that this code does not use external precalculated prime tables. Time taken on my good old Thinkpad T60, so I think it deserves to be called fast. Tested with mawk and gawk on Debian8/AMD64 • good 129 bytes in gawk: now with Debian10/AMD64 on my corei7-i870@3.6Ghz: real 0m2,417s user 0m2,205s sys 0m0,042s – JeanClaudeDaudin Jan 25 at 14:42 • you can save one byte with: BEGIN{n=2;i=0;while(n<1366662){if(n in L){p=L[n];del L[n]}else{P[p=n]=++i;if(i in P)print n}j=n+p;while(j in L)j+=p;L[j]=p;n++}} – JeanClaudeDaudin Jan 25 at 15:08 # CJam, 19 3D#{mp},_1e4<:(\f=p  You can try it online, but you'll need a little patience :p For the record, the last number is 1366661. # Perl, 55 bytes use ntheory':all';forprimes{print nth_prime$_,\$/}104729


Uses @DanaJ's Math::Prime::Util module for perl (loaded with the pragma ntheory). Get it with:

cpan install Math::Prime::Util
cpan install Math::Prime::Util::GMP


# 05AB1E, 7 bytes (non-competing)

Code:

4°L<Ø<Ø


Try it online!, note that I have changed the 4 into a 2. If you have a lot of time, you can change the 2 back to 4, but this will take a lot of time. I need to fasten the algorithm for this.

Explanation:

4°       # Push 10000 (10 ^ 4)
L      # Create the list [1 ... 10000]
<     # Decrement on every element, [0 ... 9999]
Ø    # Compute the nth prime
<   # Decrement on every element
Ø  # Compute the nth prime