## Large prime numbers

As of this writing, the largest known prime number is 277232917-1. “Not large enough!”, I say. Let’s write a computer program that will print out a larger prime number. No, not a program that will take billions of years to run. It should take only a few minutes.

There is a small catch. Hardly worth mentioning, really, but… Our program will print out more than one number. And it won’t tell which of them are prime. And I can’t absolutely guarantee that it will succeed. We’ll settle for a 99.9% chance of success, though we can easily make the probability much higher if desired.

Smaller numbers are more likely to be prime, so our hunting ground will be the numbers just slightly larger than 277232917-1. The probability that a randomly chosen number about as big as $N$ is prime is $\frac{1}{\ln(N)}$. To calculate $\ln(2^{77232917})$, multiply the exponent by $\ln(2)$, to get 53533778.66. So, of the integers in the vicinity of 277232917, about one in every 53.5 million is prime.

How many such integers would we have to look at to have a 99.9% of at least one being prime? This can be calculated directly, but instead let’s use the fact that, given $T$ trials, each with a $\frac{1}{T}$ chance of success, the probability of none of them succeeding is (for large $T$) very close to $\frac{1}{e}$: about 36.79%.

Each of our trials has a 1 in 53.5 million chance of success, so with 53.5 million trials, our chance of failure is 36.79%. How many of these sets of (53.5 million in this case) trials do we need, to get down to our desired 0.1% chance of failure? You can multiply 0.3679 by itself a few times with a calculator, to see that 7 times is sufficient. (The actual number is about 6.9.)

Multiply 53.5 million by 7 to get 374.5 million. That’s how many numbers we want to print. Here’s a Python program to do that:

#!/usr/bin/env python3
howmany = 374500000
for i in range(howmany):
print(f'2^77232917+{i}')


Its output starts like this:

2^77232917+0
2^77232917+1
2^77232917+2
2^77232917+3
2^77232917+4
...

That’s nice, but let’s go the extra mile, and be a little more efficient. There’s no point in printing out even numbers, for example. Hardly any of them are prime. Let’s use a sieve to filter out numbers that are a multiple of 2, 3, and/or 5.

In any block of 30 consecutive integers, 15 are divisible by 2. Another 5 are divisible by 3, and another 2 are divisible by 5. That’s 22 numbers out of 30 that are definitely composite, leaving 8 that might be prime. Let’s only print those 8 numbers. 374.5 million multiplied by $\frac{8}{30}$ leaves us with 99.866 million numbers to print, so let’s make it an even 100 million.

But it’s not that easy, because we still have to properly align our sieve. To do that, we need to know the remainder when 277232917 is divided by 30, and unfortunately most calculators can’t handle numbers with 23 million digits.

Fortunately, the powers of 2 (mod 30) quickly fall into a simple loop of length 4: 2, 4, 8, 16, 2, 4, 8, 16… If $N$ (mod 4) is 1, then $2^N$ (mod 30) is 2. 77232917 (mod 4) is 1, so our answer is 2.

Working mod 30, the numbers that could be prime are 1, 7, 11, 13, 17, 19, 23, 29. To program the sieve, we’ll make a list of the differences between consecutive numbers in that list: 6, 4, 2, 4, 2, 4, 6, 2. The first 6 is the difference between 1 and 7, and so on. The final 2 is because 29+2 wraps around to 1.

The very first number we’ll print is the first number greater than or equal to 277232917 that we haven’t ruled out as being prime: 277232917+5. That number is equivalent to 7 (mod 30). So the first item in our differences list that we’ll use is the second item in the list (index 1, value 4).

Here’s a program that implements the sieve:

#!/usr/bin/env python3
howmany = 100000000
d = [6,4,2,4,2,4,6,2]
di = 1
n = 5
for i in range(howmany):
print(f'2^77232917+{n}')
n += d[di]
di += 1
if di>=len(d):
di = 0


Its output starts and ends like this:

2^77232917+5
2^77232917+9
2^77232917+11
2^77232917+15
2^77232917+17
...
2^77232917+374999985
2^77232917+374999987
2^77232917+374999991
2^77232917+374999997
2^77232917+374999999

I ran the program, with its output directed to a file, and it took about 3 minutes to produce a 2 GB output file. A big file, but not extraordinarily big. On the off chance that all my math is correct, we can expect the file to have approximately 7 prime numbers in it larger than the largest known prime. There’s a 99.9% chance that it has at least 1.

If there’s a point to this post, it’s that prime numbers, even really big ones, are not as rare as one might think. The difficulty is in figuring out whether a really big number is or is not prime.