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.
To convert a picture to a number, we have to take a few simple intermediate steps.
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.
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!
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.
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.
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:
sqrt(n)
(we use integer square root^{1} in this case), it must not be divisible by that numberThis 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.
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.
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")
print
statement is shown to the user.size_y-1
(the number of lines):
maybeprime
as a string, and take the elements between (i*sx)
and ((i+1)*sx)
.\n
) on the end of the string we just made.finally:
block before going to the print
statement to the userfinally:
block, and close the file. If you don’t close the file, will puppies die? Every time.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
An example of integer square root in Python, as the primes we are playing with are too big for floating point:
def fisqrt(n):
"""Fast Integer Square Root"""
xn = int(1)
xn1 = (xn + n // xn) // 2
while abs(xn1 - xn) > 1:
xn = xn1
xn1 = (xn + n // xn) // 2
while xn1 * xn1 > n:
xn1 -= 1
return xn1