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


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.

  • 3
    \$\begingroup\$ 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. \$\endgroup\$ Oct 23, 2015 at 17:43
  • 6
    \$\begingroup\$ Optimizing for runtime is not what we do in a code-golf challenge; the shortest program always wins. \$\endgroup\$
    – lirtosiast
    Oct 23, 2015 at 18:53
  • 1
    \$\begingroup\$ Primes with prime subscripts: A006450. \$\endgroup\$
    – user12166
    Oct 23, 2015 at 21:55
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ Oct 24, 2015 at 12:04
  • 2
    \$\begingroup\$ 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. \$\endgroup\$
    – Dennis
    Oct 25, 2015 at 14:10

16 Answers 16


MATLAB/Octave, 25 bytes


It doesn't get much more straightforward than this.


Python, 72 bytes

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.

while p<104730:
 if len(l)in P%p*l:print p
  • 1
    \$\begingroup\$ This is way better than the garbage I was writing. +1 \$\endgroup\$
    – user45941
    Oct 23, 2015 at 17:51
  • \$\begingroup\$ This only prints 1229 numbers \$\endgroup\$ Oct 23, 2015 at 19:03
  • \$\begingroup\$ @aditsu I think I see my mistake. Are you able to run this code with the bigger bound? \$\endgroup\$
    – xnor
    Oct 23, 2015 at 19:07
  • \$\begingroup\$ It will probably take a long time :p \$\endgroup\$ Oct 23, 2015 at 19:16
  • \$\begingroup\$ I think it finished \(@;◇;@)/ , it seems correct \$\endgroup\$ Oct 24, 2015 at 15:59

J, 11 bytes


Outputs the primes in the format

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


        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, 26 25 23 bytes


Pure function returning the list.

  • 1
    \$\begingroup\$ Prime is Listable so a simple Prime@Prime@Range@1*^4& will do \$\endgroup\$
    – user46060
    Oct 23, 2015 at 22:34
  • \$\begingroup\$ I know the feeling ... In any case, I think this is the prettiest Mathematica solution I have seen on here! \$\endgroup\$
    – user46060
    Oct 23, 2015 at 22:36
  • \$\begingroup\$ 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! \$\endgroup\$
    – user46060
    Oct 23, 2015 at 23:14

Haskell, 65 bytes

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.


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++}}


$ 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


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

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.


$ time awk -f x.awk | wc -l

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

  • \$\begingroup\$ good 129 bytes in gawk: now with Debian10/AMD64 on my [email protected]: real 0m2,417s user 0m2,205s sys 0m0,042s \$\endgroup\$ Jan 25, 2020 at 14:42
  • \$\begingroup\$ 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++}} \$\endgroup\$ Jan 25, 2020 at 15:08

PARI/GP, 25 bytes


CJam, 19


You can try it online, but you'll need a little patience :p

For the record, the last number is 1366661.


Husk, 11 10 bytes

(or 6 5 bytes to list the first 10,000 prime-indexed prime numbers, followed by [the infinite number of] all the rest of them)

Thanks to user and Razetime for suggesting the 6 5-byte version!


Try it online! (unfortunately times-out on TIO before printing all 10,000 prime-indexed primes, but at least outputs what it's found so far...)

↑               # print the first
 !4             # 4th index of
   İ⁰           # powers of 10 (so, the first 10,000), of
     Ṡ          # hook: apply function to its own argument
      m!        # map index function to
        İp      # prime numbers
                # (implied by hook) to prime numbers
  • \$\begingroup\$ Nice answer. The 10000 primes restriction is a shame - too many bytes are wasted on that. \$\endgroup\$
    – user
    Nov 18, 2020 at 15:22
  • 1
    \$\begingroup\$ If you output an infinite list, it will still have the first 10000 items ;) \$\endgroup\$
    – Razetime
    Nov 18, 2020 at 15:27
  • \$\begingroup\$ @user & Razetime - thanks for the suggestion from both of you - I've incorporated it (& credited you both)! \$\endgroup\$ Nov 18, 2020 at 15:45
  • \$\begingroup\$ It didn't actually occur to me to print an infinite list - that's all Razetime, but thanks anyway. \$\endgroup\$
    – user
    Nov 18, 2020 at 16:51

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

Japt, 11 bytes

Èj ©°Tj}jL²

Test it


Python 3, 76 109 bytes - disqualified - requalified?

edit: I am sorry, I was going for primes but missed the putting, didn´t reach the goal of "primes with primeindex". Your are free to grab and enhance. I'll be glad to be learning from you.

Based on Wilson Theorem

def p(x):
 for i in range(3,x):
  if c%i!=0:

without index, 62 bytes

def p(x):
 for i in range(3,x):
  if c%i!=0:print(i)

Anyone can give me an example or good link, how to make a lambda function from it?

edit: just to have it done right: 109bytes

def p(x):
 for i in range(2,x):
  if c%i!=0:
   if n in p: print(p[-1])
  • \$\begingroup\$ I think you've misunderstood the question. Your functions seem to output all prime numbers (and the first one precedes them with an index). The question defines 'prime-indexed prime numbers' as those in the list of all primes that have an index value that is also a prime. For instance, 7 - the 4th prime - isn't prime-indexed because 4 isn't a prime... \$\endgroup\$ Nov 21, 2020 at 16:25
  • \$\begingroup\$ I am sorry. You are right, I misunderstood it :( Thank you very much for spending your time on clearing this up for me! \$\endgroup\$
    – Bambeleme
    Nov 22, 2020 at 0:59

Vyxal H, 8 7 6 5 bytes

Saved 1 + 1 bytes thanks to Aaron Miller


Try it Online!

H flag - Push 100 to the stack
²        Squared (10,000)
 ʀ       Range [0, 10,000)
  ǎ      Map each a in that range to the ath prime
   ‹     Decrement each, because the challenge uses 1-indexing
    ǎ    Get the nth prime for each of those
  • 1
    \$\begingroup\$ 7 bytes \$\endgroup\$ May 3, 2021 at 2:26
  • \$\begingroup\$ @AaronMiller Thanks, forgot about auto-closing! \$\endgroup\$
    – user
    May 3, 2021 at 2:27
  • 1
    \$\begingroup\$ I completely forgot about the H tag, which presets the stack to 100, allowing for 5 bytes \$\endgroup\$ May 3, 2021 at 15:31

05AB1E, 7 bytes (non-competing)



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.


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

Jelly, 6 bytes


Try it online!

Obviously, this times out on TIO. Here is a version that prints as it goes, which reaches \$57493\$ (the \$765\$th term) before timing out

How it works

ȷ4RÆN⁺ - Main link. Takes no arguments
ȷ4     - 10000
  R    - [1, 2, 3, ..., 9999, 10000]
     ⁺ - Do twice:
   ÆN  -   n'th prime of each

05AB1E, 6 bytes


Try it online!

     Ø  # prime numbers with...
    <   # 1-based...
     Ø  # indices in...
  Åp    # first...
т       # 100...
 n      # ^ 2...
  Åp    # primes
        # implicit output

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.