Code Golf Stack Exchange is a question and answer site for programming puzzle enthusiasts and code golfers. It only takes a minute to sign up.
Sign up to join this community
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.
p=primes(2e6)
p(p(1:1e4))
It doesn't get much more straightforward than this.
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
p:<:p:i.1e4
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
Prime@Prime@Range@1*^4&
Pure function returning the list.
Listable so a simple Prime@Prime@Range@1*^4& will do
\$\endgroup\$
– user46060
Oct 23 '15 at 22:34
@ 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 '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.
apply(prime,primes(10^4))
...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
Readable:
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
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.
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
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
