What is an emirp? Can we find a new world record size prime of this form, using only MATLAB, and a home computer?
An emirp is a prime that is prime when viewed in in both directions. They are not too difficult to find at a lower level. For example...
isprime([199 991])
Gosh, that was easy. But what happens if the number is a bit larger? The problem is, primes themselves tend to be rare on the number line when you get into thousands or tens of thousands of decimal digits. And recently, I read that a world record size prime had been found in this form. You have probably all heard of Matt Parker and numberphile.
And so, I decided that MATLAB would be capable of doing better. Why not? After all, at the time, the record size emirp had only 10002 decimal digits.
How would I solve this problem? First, we can very simply write a potential emirp as
10^n + a
then we can form the flipped version as
ahat*10^(n-d) + 1
where ahat is the decimally flipped version of a, and d is chosen based on the number of decimal digits in the number a itself. Not all emirps will be of that form of course, but using all of those powers of 10 makes it easy to construct a large number and its reversed form. And that is a huge benefit in this. For example,
Pfor = sym(10)^101 + 943
Prev = 349*sym(10)^99 + 1
It is easier to view these numbers using a little code I wrote, one that redacts most of those boring zeros.
emirpdisplay(Pfor)
emirpdisplay(Prev)
And yes, they are both prime, and they both have 102 decimal digits.
isprime([Pfor,Prev])
Sadly, even numbers that large are very small potatoes, at least in the world of large primes. So how do we solve for a much larger prime pair using MATLAB?
The first thing I want to do is to employ roughness at a high level. If a number is prime, then it is maximally rough. (I posted a few discussions about roughness some time ago.)
https://www.mathworks.com/matlabcentral/discussions/tips/879745-primes-and-rough-numbers-basic-ideas
In this case, I'm going to look for serious roughness, thus 2e9-rough numbers. Again, a number is k-rough if its smallest prime factor is greater than k. There are roughly 98 million primes below 2e9.
The general idea is to compute the remainders of 10^12345, modulo every prime in that set of primes below 2e9. This MUST be done using int64 or uint64 arithmetic, as doubles will start to fail you above
format short g
sqrt(flintmax)
The sqrt is in there because we will be multiplying numbers together here, and we need always to stay below intmax for the integer format you are working with. However, if we work in an integer format, we can get as high as 2e9 easily enough, by going to int64 or uint64.
sqrt(double(intmax('int64')))
And, yes, this means I could have gone as high as primes(3e9), however, I stopped at 2e9 due to the amount of RAM on my computer. 98 million primes seemed enough for this task. And even then, I found myself working with all of the cores on my computer. (Note that I found int64 arithmetic will only fire up the performance cores on your Mac via automatic multi-threading. Mine has 12 performance cores, even though it has 16 total cores.)
I computed the remainders of 10^12345 with respect to each prime in that set using a variation of the powermod algorithm. (Not powermod itself, which was itself not sufficiently fast for my purposes.) Once I had those 98 millin remainders in a vector, then it became easy to use a variation of the sieve of Eratosthenes to identify 2e9-rough numbers.
For example, working at 101 decimal digits, if I search for primes of the form 10^101+a, with a in the interval [1,10000], there are 256 numbers of that form which are 2e9-rough. Roughness is a HUGE benefit, since as you can see here, I would not want to test for primality all 10000 possible integers from that interval.
Next, I flip those 256 rough numbers into their mirror image form. Which members of that set are also rough in the mirror image form? We would then see this further reduces the set to only 34 candidates we need test for primality which were rough in both directions. With now only a few direct tests for primality, we would find that pair of 102 digit primes shown above.
Of course, I'm still needing to work with primes in the regime of 10000 plus decimal digits, and that means I need to be smarter about how I test a number to be prime. The isprime test given by sym/isprime only survives out to around 1000 decimal digits before it starts to get too slow. That means I need to perform Fermat tests to screen numbers for primality. If that indicates potential primality, I currently use a Miller-Rabin code to verify that result, one based on the tool Java.Math.BigInteger.isProbablePrime.
And since Wikipedia tells me the current world record known emirp was
117,954,861 * 10^11111 + 1 discovered by Mykola Kamenyuk
that tells me I need to look further out yet. I chose an exponent of 12345, so starting at 10^12345. Last night I set my Mac to work, with all cores a-fumbling, a-rumbling at the task as I slept. Around 4 am this morning, it found this number:
emirp = @(N,a) sym(10)^N + a;
Pfor = emirp(12345,10519197);
Prev = sym(flip(char(Pfor)));
emirpdisplay(Pfor)
emirpdisplay(Prev)
isProbablePrimeFLT([Pfor,Prev],210)
I'm afraid you will need to take my word for it that both also satisfy a more robust test of primality, as even a Miller-Rabin test that will take more time than the MATLAB version we get for use in a discussion will allow. As far as a better test in the form of the MATLAB isprime utility to verify true primality, that test is still running on my computer. I'll check back in a few hours to see if it fininshed.
Anyway, the above numbers now form the new world record known emirp pair, at 12346 decimal digits. Yes, I do recognize this is still what I would call low hanging fruit, that having announced a largest prime of this form, someone else willl find one yet larger in a few weeks or months. But even so, for the moment, MATLAB owns the world record!
If anyone else wants a version of the codes I used for the search, I've attached a version (emirpsearchpar.m) that employs the parallel processing toolbox. I do have as well a serial version which is of course, much, much slower. It would be fun to crowd source a larger world record yet from the MATLAB community.
3 Comments
Time Descending
After having gotten this last prime pair to the point where I was satisfied all tests certified they were indeed prime, well, it was a day after Pi day. So I set all 16 cores to work once again. (And fired up my system fan too.) Now I was looking for primes of the form
10^14159 + a
and their mirror image counterparts. At this level, even a single Fermat test for primality requires around 20 seconds to execute. A couple of days later, my Mac finally gave a burp to tell me it had found an interesting pair, this time with over 14000 decimal digits.
a = 18010681;
ahat = 18601081;
Intriguingly, a and its mirror image, ahat look vaguely palindromish themselves, but nothing significant to be learned there I think.
Pa = sym(10)^14159 + a;
emirpdisplay(Pa)
Pa =
100000... (14141 zeros redacted) ...0000018010681
Pahat = ahat*sym(10)^14152 + 1;
emirpdisplay(Pahat)
Pahat =
1860108100000... (14141 zeros redacted) ...000001
Now I need to post-test them, to verify they are indeed prime. First, 9 independent Fermat tests, with 9 different witnesses. Square free witnesses are the best witnesses, and so prime witnesses make the most sense to me, even though composite witnesses are good too as long as they are square free. So 210 is an entirely good witness for a Fermat test. Regardless, this set is an excellent first test.
isProbablePrimeFLT(Pa,[2 3 5 7 11 13 17 19 23])
ans =
logical
1
isProbablePrimeFLT(Pahat,[2 3 5 7 11 13 17 19 23])
ans =
logical
1
Next, is a MIller-Rabin test, employing the Java.Math.BigInteger.isProbablePrime utility. isProbablePrimeJ is my wrapper for that call. This will be considerably slower than even 9 separate Fermat tests, (maybe 5 minutes apiece per test) and I need to do two of them to verify this new pair. Each of those tests is set to do 100 iterations.
isProbablePrimeJ(Pa)
ans =
logical
1
isProbablePrimeJ(Pahat)
ans =
logical
1
At this point, both have satisfied a fairly strict requirement for primality with Miller-Rabin, but no positive poof so far. My confidence in their primality as a pair is up way over 99.99%, but the number of nines in that confidence is not relevant.
The next step is a Pocklington certificate for Pahat. This can be applied as long as I have a complete factorization of Pahat-1. That part is trivially easy to obtain, since that number is a product of many powers of 10 with 18601081.
factor(18601081)
ans =
269 69149
And so the distinct prime factors of Pahat-1 are nothing more than [2, 5, 269, 69149]. This allows me to perform a true test for primality on Pahat that will require only 4 powermod calls. isPrimePocklington does that.
isPrimePocklington(Pahat,100)
ans =
logical
1
Ah, good. I now have a certificate of primality for Pahat. But a Pocklington certificate will not apply unless I can fully factor the number 1 less than the number.
So the next step is a Fibonacci test. Here I need to use the result that is Pa is prime, then it must divide either fibonacci(pa-1) or fibonacci(pa+1). Which one of those will depend upon the remainder of pa modulo 5.
mod(Pa,5)
ans =
1
This tels me that Pa must divide fibonacci(Pa-1), IF Pa is prime. However, there are also Fibonacci liars. The code isProbablePrimeFib calls my fibmod code. (fibmod computes the remainder of a Fibonacci number, modulo some given number. This is actually not that difficult to compute, even for huge Fibonacci arguments.)
Again, this test is a little slow. I'll just call fibmod, if this call returns zero, then Pa is at worst also a Fibonacci liar, as well as a multiple Fermat liar, and a Miller-Rabin liar.
tic,fibmod(Pa-1,Pa),toc
ans =
0
Elapsed time is 125.362510 seconds.
So 2 minutes later, we know that Pa is also now far more certainly prime, since no composite numbers that fail all of these tests are known. 2 minutes of CPU time is short enough that I'll try it on Pahat too, even though I'm not at all worried about that one after the positive result from the Pocklington certificate.
isProbablePrimeFib(Pahat)
ans =
logical
1
I'm actually quite confidant that both members of this new pair are prime. But the kicker is I'm not yet 100% certain about Pa. I've added a half dozen or so nines to my confidance level. But I'm still looking for
isprime(Pa)
...
Sadly, I expect this last call to take between 24 and 48 hours to finish. It was around 24 hours at 12346 decimal digits. So at 14160 decimal digits, I'm guessing 48 hours. And that call is not, and cannot be easily parallelized.
Anyway, it looks like I've now surpassed my previous record.
What kind of Mac do you have?
A truly pretty thing about this result is it employed two forms of parallel processing in MATLAB, initially massive multi-threading, working in int64 arithmetic, followed by the use of parfor to perform Fermat tests across all the cores of my computer.
I have since verified that
10^12338*79191501 + 1
is certainly prime, by way of a Pocklington certificate. This test requires the completely known factorization of X-1, which is easily done here.
The cousin of that number was more difficult to prove.
10^12345 + 10519197
Here I used a Fibonacci test, which did yield a positive result for primality. And while Fibonacci tests do admit pseudo-primes, they are somewhat rare, and in combination with a Fermat test, suggests the number is almost certainly prime.
I'll take another shot at it tonight, unless someone else has Mathematica and can put it into ProvablePrimeQ. (The MATLAB isprime has been running continuously on this for over 12 hours now.) Edit: After roughly 24 hours of effort, isprime came back with a finding of prime for the latter form. And isprime is a 100% test as I recall, so both halves of this emirp pair have been conclusively proven to be prime now.
I'll probably next start a search in the 20k-30k digit range, and hope my CPU can stand the heat.
Sign in to participate
