https://projecteuler.net/problem=7
By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10,001st prime number?
Solution 1
Use Eratosthenes’ Sieve to generate enough primes. Then pick the 10,001st.
The Sieve: start with an initial state of known primes (2) and a vector flagging the (odd) numbers still to test.
Here it is in slo-mo.
And in code:
q)kc:{(1#2;0b,1_x#10b)} / (known primes; flag candidates)
q)kc 20
,2
00101010101010101010b
The step function: sieve the next prime. It’s the first number flagged. Join it to the known primes and remove its multiples from the flag vector.
q)snp:{(x,n;y&count[y]#i<>til n:1+i:y?1b)}. / sieve next prime
q)snp kc 20
2 3
00001010001010001010b
Sieve the initial state to get all the primes under a million. (Surely enough?)
q)es:{{x>last first y}[floor sqrt x] snp/kc x} / Eratosthenes' Sieve
q)es 50
2 3 5 7
00000000001010001010001000001010000010001010001000b
q)rslv:raze @[;1;1+where@]@ / resolve result
q)rslv es 50
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
q)pt:rslv es@ / primes to
q)pt 50
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
Solution 2
The While form of the Over iterator here uses a test function that stops when candidates exceed the square rot of N. All the remaining flags then mark primes, and we know that we have sieved out all the primes below N.
But until then we do not know how many of the remaining flags mark primes. The first item of the state lists the primes we have found.
If we examine the intermediate result of sieving for primes below a million, we see there have been 168 iterations.
q)count first es 1000000 / # iterations
169
Sieving for primes below a million was a lucky guess. Is there a lower N we could use?
By convention π(x) denotes the number of primes below x. Its values have been calculated for powers of 10.
PI:0 4 25 168 1229 9592 78498 664579 5761455 50487534 455052511 4118054813 37607912018 346065536839
This tells us that to get 10,001 primes the lowest power of 10 to use for N is 106.
q)sum PI<10001
6i
q)count pt "j"$10 xexp sum PI<10001
78498
We have replaced a lucky guess with a reasoned approximation. But it’s still the same number of iterations, and seven times as many primes as we need.
Solution 3
The Prime Number Theorem tells us that
π(x)≈x÷ln(x)
What value of x gives us π(x)>10000?
Let’s look. Start with 1000 and, instead of going up in orders of magnitude, keep doubling.
pi:{x%log x} / π(x) first approximation
q)(10000>pi@)(2*)/1000
128000
q)pi 128000
10884.55
q)count pt (10000>pi@)(2*)/1000
11987
We could refine the search to find a lower value of N but it is unlikely to reduce the number of iterations by much.
q)count first es 128000 / # iterations
72
Comments
This is a nice example of using q in the REPL not only to find a solution but to explore how to improve it.
kc:{(1#2;0b,1_x#10b)} / (known primes; flag candidates)
snp:{(x,n;y&count[y]#i<>til n:1+i:y?1b)}. / sieve next prime
es:{{x>last first y}[floor sqrt x] snp/kc x} / Eratosthenes' Sieve
rslv:raze @[;1;1+where@]@ / resolve result
pt:rslv es@ / primes to
pi:{x%log x} / π(x) first approximation
p:pt (10000>pi@)(2*)/1000 / first 10000 or so primes
p@10000 / 10001st prime
Contributors
- Stephen Taylor
