MATLAB: Isprime function seems to have poor performance

isprimeMATLABperformanceprimeprime numberprime testingsymsymbolicSymbolic Math Toolbox

Why is MatLab's "isprime" function so much slower than Octave's "isprime" function?

I am using MatLab's "isprime" function to check whether a large number is a prime or not using the symbolic toolbox. I found that the performance of "isprime" in MatLab is much slower than in Octave. Why is this the case or what am I doing wrong with MatLab?

My tests with octave testing large Mersenne-primes produced the following runtimes for "isprime":

tested prime	    runtime in seconds
    2^607-1,          0.15724        
    2^1279-1,         0.18309      
    2^2203-1,         0.41784        
    2^2281-1,         0.70215        
    2^3217-1,         1.7013        
    2^4253-1,         2.4854        
    2^4423-1,         2.2523        
    2^9689-1,         25.7571        
    2^9941-1,         25.761        
    2^11213-1,        38.3376

and with MatLab's "isprime":

tested prime	    runtime in seconds
2^607-1                  31.930225            
2^1279-1                 547.414940            
2^2203-1                 5168.567632            
2^2281-1                 5578.169207            
2^3217-1                 461.535261            
2^4253-1                 739.918345             
2^4423-1                 3805.209681             
2^9689-1                 8954.457005             
2^9941-1                 10550.740359             
2^11213-1                11865.530147

MatLab's documentation about "isprime" says, that 10 random tests based on the Miller-Rabin method are done. I believe, Octave only does 4 random tests (I suppose also Miller-Rabin, but I am not sure). However this does by far not explain the huge difference in runtime.

In both cases no parallelisation was used and the program ran on one thread on the CPU.

This is the MatLab code I used to run the test. The Octave version is basically the same…

function speedtrace_isprime();
  
%  teste Dauer der Ausführung von "isprime" in Abhängigkeit von wachsenden Mersenne-Primzahlen 2^p-1
%  zu testende p's:  
  p = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941, 11213, 19937, 21701, 23209,   44497, 86243, 110503, 132049, 216091, 756839, 859433, ...
  1257787, 2976221, 3021377, 6972593, 13466917, 20996011, 24036583, 25964951, 30402457, 32582657, ...
  37156667, 42643801, 43112609, 57885161, 74207281, 77232917, 82589933];
  
  speedtrace = fopen('speedtrace_isprime.txt', 'w');                                %trace-file öffnen
  fprintf(speedtrace, "%s %s \n", "Start: ", string(datetime));                      % schreiben

  fprintf(speedtrace, "%s \n", "getestete Primzahl,   Zeit in Sekunden,   Uhrzeit/Datum");   % schreiben
  fclose(speedtrace);                                                               % file wieder zumachen
  
  base = sym("2");
  disp(["getestete Primzahl,   Zeit in Sekunden,    Uhrzeit/Datum"]);  
  for k = 1:1:numel(p);
    tic;
    isprime(base^p(k)-1);
    Zeit(k) = toc;
    fprintf("2^%i-1                 %f             %s  \n", p(k), Zeit(k), datetime);
    
    speedtrace = fopen('speedtrace_isprime.txt', 'a');
    fprintf(speedtrace,"2^%i-1                 %f            %s  \n", p(k), Zeit(k), datetime);
    fclose(speedtrace);
      
%   figure(1); plot(Zeit);
end; 
 
  speedtrace = fopen('speedtrace_isprime.txt', 'a');
  fprintf(speedtrace, "%s %s \n", ["Ende: ", string(datetime)]);
  fclose(speedtrace);  
    
end

Best Answer

First, using isprime to test for primality of a Mersenne number is a bad idea. The Lucas-Lehmer test, is as I recall, much more efficient here, and it gives a statement of primality, NOT a statement of isprobable primality. The Miller-Rabin test is an isprobable test, not a proof of primality! That is why multiple internal runs are employed. Do some reading about how these methods can fail.

But more imprtantly, read about Lucas-Lehmer.

https://en.wikipedia.org/wiki/Lucas–Lehmer_primality_test

In fact, sym/isprime has greatly improved over previous ersions. I recall testing it some years ago, and it was pretty slow. But that is not so true now. This means if you have an old release, it might gain you to upgrade. By the way, if I do a test for primality of a number with on the order of 70000 digits, it is a matter typically of 20 minutes, as I recall. (Its been a month or so since I was running massive overnight jobs, searching for prime members of a specific family of primes, thus quasi-modified Woodal primes of a specific class. They can be exceedingly rare if you choose the right family.)

Anyway, you don't want to use isprime to test for a Mersenne prime. Use Lucas-Lehmer. It is REALLY fast in comparison. For a quick implementation...

function s = Lucas_Lehmer(p)
  % Returns true for prime 2^p-1
  
  % don't even bother if p is not prime
  if ~isprime(p)
    s = false;
    return
  end
  
  s = 4;
  M = sym(2)^p - 1;
  for i =  1:p-2
    s = mod(s*s - 2,M);
  end
  s = logical(s == 0);

Seriously, you can't do something much simpler.

Lucas_Lehmer(607)
ans =
  logical
   1
   
isprime(sym(2)^607-1)
ans =
  logical
   1
   
timeit(@() Lucas_Lehmer(607))
ans =
           0.4162580338225
tic,isprime(sym(2)^607-1),toc
ans =
  logical
   1
Elapsed time is 9.295827 seconds.
Lucas_Lehmer(613)
ans =
  logical
   0
% A quick check of the list of known Mersenne primes tells me that 21701 is...
tic,Lucas_Lehmer(21701),toc
ans =
  logical
   1
Elapsed time is 17.560138 seconds.

2^21701-1 is not really that large of a number, but I expect that isprime will take considerably longer. (My own test right now is still running after about 10 minutes.)

tic,isprime(sym(2)^21701-1),toc
???

Finally, if your goal is to find seriously large Mersenne primes, MATLAB is probably not the tool to use if you are looking for primes with millions of digits at some point. But I have found it to work reasonably well in my own play time, searching for primes with 50000-100000 decimal digits. For any seriously large investigation, I'd suggest looking at GIMPS, which seems to be the state of the art.

If your goal is to search for primes in some other family, then there are sometimes tricks to make things more efficient, but they are often strongly dependent on the number family itself. For example, Cullen and Woodall number families have some neat things you can do, to avoid calling isprime too often.

Related Solutions

MATLAB: Is the function isprime deterministic

Your problem is that number is NOT represented in MATLAB as you think it is. MATLAB CANNOT represent a number as a double exactly when the integer is larger than 2^53-1. Double precision is limited in that respect. And yes, I know you THINK you created the number 2^64-59 as a uint64. You did not.

So, YES, I believe it can be shown that isprime is deterministic within the limits of double precision. However, NO, if you throw numbers at it that are too large, then it will fail. The biggest reason why isprime failed for you was that you created it as a double precision number, AND only then converted to uint64. 2^64 is NOT a number you can trust, because it is created as a double.

There is another problem, even more serious. Can you even compute 2^64, and store it in a uint64 variable? NO!

uint64(2)^64
ans =
  uint64
   18446744073709551615
>> intmax('uint64')
ans =
  uint64
   18446744073709551615

Should 2^64 be an odd number? Of course not, yet as you see, 2^64 overflows even a uint64 variable.

Instead, use sym, which her will be quite efficient.

isprime(sym(2)^64-59)
ans =
  logical
   1

See that I created 2 as a symbolic number, only then raising it to a power.

If I did not have the symbolic toolbox, I could have done this - subtly avoiding any overflows, because I recognize that intmax('uint64') is just 2^64-1.

isprime(intmax('uint64') - 58)
ans =
  logical
   1

Now it works. the symbolic toolbox isprime is MUCH faster, although I expect it is really best called isprobable prime for large integers, since it will be based on a multiple set of calls to a MIller-Rabin test, and possibly other tricks. I recall looking up the basis of the mathematica isprime test, which is similar. A true test for primality that cannot be fooled is more difficult to achieve, and will be MUCH slower running most of the time.

By the way, what number did you test to see if it was prime?

uint64(2^64) - 59
ans =
  uint64
   18446744073709551556
   
intmax('uint64') - 58
ans =
  uint64
   18446744073709551557

You should see clearly why MATLAB thought the former version was not prime. After all, most even numbers are not prime. Almost all of them, in fact.

MATLAB: I wrote a code to produce twin primes, however it is pretty strenuous for the computer and I need some high numbers to check a conjecture. Is there a better way to get twin primes. Note: these twin primes are represented by 6m’s where 6m+-1 are pri

Ok, it looks like you are still interested in this, so here is an approach targetted more at what I hear from your last comment.

Simple would be, if you just want to find all twin primes less than some value, just use primes to get all primes, then find all primes that are exactly 2 larger than the one before.

N = 1e7;
P = primes(N);
ind = find(diff(P) == 2);
TP = [P(ind)',P(ind+1)'];
size(TP)
ans =
       58980           2
TP(1:10,:)
ans =
     3     5
     5     7
    11    13
    17    19
    29    31
    41    43
    59    61
    71    73
   101   103
   107   109

But there is no formula to find all such twin primes, nor to predict when a twin prime pair will occur, except that it must be a pair of the form 6*m+[-1 1].

Yes, you could use isprime to test each pair of numbers of that form, but that would be extremely inefficient, far less good than just calling primes directly.

So, is there something slightly better? For example, could we write a variation of a sieve, targetted specifically at twin primes? For example, suppose I wanted to find all twin primes less than 1e9? primes(1e9) is probably too big for me to do on my computer. But I know that twin primes must be of the form 6*m +/-1.

tic
% find all twin prime pairs where both are less than N
N = 1e9;
Plist = primes(sqrt(N));
% Keep the small twin primes from Plist
twin_primes = Plist(diff(Plist) == 2);
twin_primes = [twin_primes;twin_primes+2]';
% We don't need to check 6*m-1 for divisibility by 2 or 3.
Plist(1:2) = [];
 
% m will correspond to the index of a potential twin prime, 6*m-1.
% At the end of this comutation, if m(k) is true, then
% 6*m(k) +/- 1 will be a twin prime. Note that all of this will miss the
% very first twin prime pair, [3 5]. But all other twin prime pairs are of
% the form 6*m +/- 1.
m = true(1,floor((N-2)/6));
 
for p_i = Plist
  % which elements of m correspond to 6*m-1 that are divisible by p_i?
  % Find them using the modular inverse of 6, modulo p_i.
  p_i_inv = minv(6,p_i);
  
  m1 = p_i_inv;
  if p_i == (6*m1 - 1)
    m1 = m1 + p_i;
  end
  m(m1:p_i:end) = false;
  
  m2 = mod((p_i - 1)*p_i_inv,p_i);
  
  if p_i == (6*m2 - 1)
    m2 = m2 + p_i;
  end
  m(m2:p_i:end) = false;
  
end
m = find(m);
 
twin_primes =[twin_primes;[6*m-1;6*m+1]'];
toc

So reasonably fast.

Elapsed time is 9.773307 seconds.

This ran pretty quickly on my computer, since the only RAM it required was a vector of size N/6 bytes. So for N=1e9, m only required 0.17 gigabytes. And the loop was also pretty efficient.

Note that the above code uses my minv tool, which is found in my vpi toolbox, although when used on doubles, it returns double precision output.

whos m twin_primes
  Name                   Size                    Bytes  Class     Attributes
  m                      1x3424019            27392152  double              
  twin_primes      3424506x2                  54792096  double

So it found every one of the 3.4 million twin prime pairs for which neither of them exceed 1e9. But since it did not need to explicitly test any of them explicitly for primality, nor did it need to generate the entire list of primes below 1e9, the computation was actually pretty fast. While my computer is moderately old, I'd guess I could easily push N out further out than 1e10. (I just tried it for N=1e10, and it took 110 seconds to find all 27412679 twin prime pairs less than 1e10. Not bad for a 10 year old computer.)

I'm not sure if this will be of any real value to you, but it will indeed generate all such twin primes less than N.

Best Answer

Related Solutions

MATLAB: Is the function isprime deterministic

MATLAB: I wrote a code to produce twin primes, however it is pretty strenuous for the computer and I need some high numbers to check a conjecture. Is there a better way to get twin primes. Note: these twin primes are represented by 6m’s where 6m+-1 are pri

Related Question