Results for
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.
