Post

Motivation

I recently found a rather nice prime number, belonging to Trinity Hall in Cambridge. The “Trinity Hall Prime” was gifted to the college by a junior research fellow (Professor McKee) and bears a likeness of the Trinity Hall coat of arms.

Since I liked how cute the number was, I decided I would give it a go for the generic case of any input image being converted into a prime number bearing its likeness. I decided on the use of Python 3 for this project.

In this short post I will show how I managed to solve this problem. My code for this project (still being worked on) is over on my GitHub.

Pictures to Numbers

To convert a picture to a number, we have to take a few simple intermediate steps.

Step 1: Resizing

By using the scipy.misc module, we have access to some basic image processing tools.

We can load an image (as a grey-scale) into memory with f = scipy.misc.imread(path, mode='L') (path being a string containing the path of the target image).

From there we can resample the grey-scale image f (represented as a numpy ndarray) with f = scipy.misc.imresize(f, (size_y, size_x), interp='biliniear'). This function takes as arguments a tuple of the y and x dimensions in pixels of our new image, which as default I set as(25,50). The y dimension is half the size of x as monospaced font is twice as tall as it is wide - together when rendered this gives me a square image.

Step 2: Assigning Numbers

The current version of my program only has two states a pixel can be in: on (8) or off (1). These two numbers were chosen as 1 takes up the least space, whereas 8 takes up the most.

To use all 10 numbers, we can define a list from darkest to lightest: pixelset = [8,9,6,5,4,3,0,2,7,1] which I decided on through inspection.

The principle is the same if we use two or ten symbols, we use threshholding to decide if a number is to be assigned that value.

I am going to use this opportunity to define an array named bounds whose elements bounds[i] and bounds[i+1] define the range of allowable values for pixelset[i].

For each row in the ndarray, we look at each pixel individually, and then consider all the possible numbers we can assign it through a for and an if statement as follows:

numimage = []

for row in image:
  for pixel in row:
    for i, pixelset_item in enumerate(pixelset):
      if (bounds[i] <= pixel) & (pixel < bounds[i+1]):
        numimage.append(pixelset_item)

We can then concatenate the list by doing a string.join and then feed the result to the integer constructor.

The full expression looks like this:

s = ''.join(map(str,numimage))
imageNumber = int(s)

Note that constructors like ‘str’ can be called without arguments - this will create the empty string (calling int yields an integer zero).

And that’s it - we should have our number to start playing around with!

Making it Prime

N.B. This section will change in the future as I’m moving towards changing digits in the mass of the image similar and not having an obvious block that was iterated over.

I had a very basic workflow for my MVP (Minimum Viable Product) release of the program. It checked only odd numbers (as all primes apart from 2 are odd) which speeds up compute times.

1. Check if the number is odd.
   • If odd, goto Step 2.
   • if not odd, increment by 1 and goto Step 2.
2. Check if the number is prime
   • If not prime, increment number by 2 and goto Step 1.
   • If prime, return the number

The hardest part of this workflow is Step 2 - checking if the number is prime! If we recall the standard size I created is 25 x 50 pixels - 1250 digits long! A large prime indeed.

Primality Testing

We can define a very basic prime checking function, is_prime_simple(n):

def is_prime_simple(n):
 """Returns 'True' if n is prime, otherwise 'False'"""
 if n == 1:
   return False
 if n > 2 and n % 2 == 0:
   return False
 maxval = (int_sqrt(n)) + 1
 for d in range(3,maxval,2):
   if n % d == 0:
     return False
 return True

The tests the number has to jump through to reach the return True at the end makes sense:

1. It must not be 1
2. It must be odd
3. For every integer up to sqrt(n) (we use integer square root¹ in this case), it must not be divisible by that number

This looks great! Computers are screamingly fast in this day and age, and Cambridge was finding numbers bigger than our poxy little 1250 digit long prime, way back when in the 1960s.

Let’s give it a go!

…and nope, this is definitely not going to work, this is far too slow! We have to ask ourselves, why is this so?

Modulo uses an integer division operation, which is very computationally expensive for chips - especially if the number is very large (greater than 64 bits - we in fact need around 2000 bits to resemble this 1250 digit long integer).

So dividing big numbers is very slow, is there a way we can turn this around and multiply small numbers? It turns out yes! The Miller-Rabin primality test is a Monte-Carlo method (numerical experiments are run and there is a statistical chance of success), which works extremely well for large numbers.

You can find examples of its implementation here and a brief description here.

Miller-Rabin is an example of a function that uses modular exponentiation to gain a significant speed advantage - it’s $\mathcal{O}({\log^3(n)})$, which is very nice.

How long until I get a prime?

I honestly can’t tell you. Statistically, the density of primes $\rho (x) \approx \frac{1}{\log (x)}$, meaning for $x \approx 10^{1250}, \quad \rho(x) \approx \frac{1}{1250}$. This implies we can expect a prime to show up once every 1250 numbers, since we only are testing odd numbers, we can half it to every 625 numbers!

If we assume on average we start halfway between the last and the next prime, we expect to have to (on average) test 313 numbers before we have our prime numbere-d image!

Of course - this is only statistical. There are arbitary gaps in prime numbers that could continue for thousands of unsuccessful trials, so don’t be upset if you don’t find a prime within the first 1000 candidates.

Wrapping up - Rastering the text.

Since I wanted this project to just spit out a single text file containing the text image you so desire, I made a small script to raster the output number into the original defined size.

The script worked as follows:

try:
  fh = open("primeout.txt","w")
  try:
    # split prime by row and save to file
    sy, sx = size
    for i in range(sy):
      fh.write(str(maybeprime)[(i*sx):((i+1)*sx)] + "\n")
  finally:
    fh.close()
except IOError:
  print("Error, Can't find file or write data")

1. We try to open the file. If this fails, an exception is thrown and a print statement is shown to the user.
2. We define the size of our image
3. We loop over a list from 0 to the size_y-1 (the number of lines):
   1. We interpret maybeprime as a string, and take the elements between (i*sx) and ((i+1)*sx).
   2. We append a newline character (\n) on the end of the string we just made.
   3. We attempt to write - if we throw an exception here we proceed to the finally: block before going to the print statement to the user
4. When finished writing, we proceed to the finally: block, and close the file. If you don’t close the file, will puppies die? Every time.

Done!

And that should be all you need to tackle the problem of generating a prime number that looks like an image! Remember the code is on the GitHub, if you want to play around with it some more, and make edits to your personal taste.

Thanks for your attention, and I’ll leave you with a prime number, 1249 digits long (the year of University College’s founding), in the shape of its crest!

11111111111111111111111188111111111111111111111111
11111111111111111111881888818811111111111111111111
11111111111111111111118888881111111111111111111111
11111888811111111111118888881118888111111111111111
11111188881111111111118888881111888811111111111111
11111188888888111111118888881111888888811111111111
11111118811188881111118888881111881118888111111111
11111111188888888881118888881111118888888888111111
11111111111888111881118888881111111188811188111111
11111111111111111111118888881111111111111111111111
11188111111111111111118888881111111111111111188111
11118888888888888888888888888888888888888888881111
11888888888888888888888888888888888888888888888811
11118888888888888888888888888888888888888888881111
11188111111111111111118888881111111111111111188111
11111111111111111111118888881111111111111111111111
11111188881111111111118888881118888111111111111111
11111118888111111111118888881111888811111111111111
11111118888888111111118888881111888888811111111111
11111118811188881111118888881111881118888111111111
11111111188888888881118888881111118888888888111111
11111111111888111881118888881111111188811188111111
11111111111111111111118888881111111111111111111111
11111111111111111111881888818811111111111111111111
1111111111111111111111118811111111111111111112257