isprime function seems to have poor performance

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

[James Tursa]

Maybe submit an enhancement request so that TMW can investigate this.

[John D'Errico]

TMW did apparently make enhancements to isprime as far as I can tell at some time in the past. It used to be slower than my VPIJ code for isprime (which is seriously faster than my VPI toolbox ) but the current release is faster in my tests.

Note that assuming isprime is using Miller_Rabin, doing more iterations will be a stronger test for primality, but also takes more time. There is a tradeoff there. (I've not checked the docs to see if that is true, but have just assumed it to be so, as Miller-Rabin seems to be the best general test. It has been a while since I was doing active reading one these things though.)

[Walter Roberson]

Octave sym appears to use Python sympy and the isprime for that is documented :

https://docs.sympy.org/latest/modules/ntheory.html#sympy.ntheory.primetest.isprime

Your inputs are over 2^64 so BPSW would be used, not Miller Rabin.

[Thomas]

Thank you very much for this information! - I'll definitely have a closer look into this!

Answer 1

[John D'Errico]

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.

[Thomas]

Hi John, Thanks for this great answer and the code...! I know, that Lucas-Lehmer is much better than Miller-Rabin for Meresenne primes. I just wanted to get prime numbers of different orders of magnitude to compare computation times. A list of Mersenne primes, was the easiest to get.

Do you think it makes sense to submit an enhancement request so that TMW can investigate in this, as was mentioned in one answer above?

[John D'Errico]

Feel free to submit the enhancement request. They may or may not decide to invest programmer time on something like this, as there are always things they need to do.

Part of the problem may be that what you see is the non-multi-threaded nature of the symbolic TB calls. Check an activity monitor to learn if Octave runs all 4 cores on those operations.

By the way, in some tests I recall that Python is considerably faster in this respect too. Subtle changes in algorithms can be a huge difference.

[Emerson Dutra]

Hi, I really like your algorithm, can you explain it to me?

[Emerson Dutra]

n = 10000;
primes = [];
if n>=2
    primes = [2];
end
for value = 3:2:n
    isprime = true;
    for divisor = 2:floor(value/2)
        if mod(value, divisor) == 0
            isprime = false;
            break;
        end
    end
    if isprime
        primes (end +1) = value;
    end
end
primes

Answer 2

[Chris]

>> tic, isprime(2^607-1), toc
ans =
  logical
   0
Elapsed time is 0.453590 seconds.

>> tic, isprime(uint64(2^1279-1)),toc
ans =
  logical
   0
Elapsed time is 26.890225 seconds.

Cant talk about the speed difference with octave but ML seemes to have trouble with the large numbers; without the uint64 ML took the second input as inf. I am using 19a and have a few year old computer.

[Thomas]

I can understand why isprime does not work with uint64. 2^607-1 is far to large to be correctly represented. This is why I use symbolic math.

According to MatLab's doc isprime is supposed to work with symbolic expressions, and it does return correct results for the prime-testings. With octave I succeeded to check primes with more than 65.000 digits (took about 57 hours). With MatLab this would run forever...

[John D'Errico]

If it took 57 hours to use isprime to test a number with 65000 digits, you need to get a seriously faster computer. As well, make sure your MATLAB release is current. Anyway, it is true that isprime on a double precison number cannot return a valid result above a relatively small number.

[Thomas]

I have the latest release (2019a). My computer is about 5 years old, 4 CPU cores, with 8 threads and I think about 3 or 4 GHz. A "normal" home-PC I would say....

How can you be so much faster? I assume, you did not use any parallelisation or GPU?

[Thomas]

Btw.: In octave this test ran in about 10 or 15 minutes. MatLab took 57 hours. I cannot really believe that the reason for this difference is only my old computer...

[John D'Errico]

Um, first of all, 4 cores is 4 cores. 8 threads is not relevant. That is just your computer making believe it has 8 cores. Hyperthreading is a fallacy in terms of CPU performance, sold to you by the manufacturers, as a way of letting you think you have real power there. If one core is running flat out, you don't gain by having two hyperthreaded cores interleaved, both running flat out.

Hyperthreading was a gain when you had only a limited set of cores (1 is the real issue here) and you want to do several things at once. It allows you to browse the web, or check your mail, while MATLAB is at work in the background. And since your computer is always doing things in the background, hyper threading was arguably a good idea in the past. Now? Just a mirage.

Next, you will probably find that a CPU monitor tells you that MATLAB is using only ONE of those cores anyway, in a call to isprime. symbolic toolbox tools are not multi-threaded at all. At least they are not so in release R2019a. (This is one thing I'd like to see change in the future for MATLAB symbolic toolbox calls.)

So you might find that part of the speed bump found by Octave is it might be multi-threaded for symbolic problems. And if they do fewer iterations of Miller-Rabin than does MATLAB, that might explain as much as an 8-1 speed bump for this specific problem.

[Thomas]

Octave is not multithreaded, the CPU behaves exactly the same as MatLab (as far as I can see this in the Windows task manager). I have exactly the same CPU usage with Octave and MatLab (= one logical porcessor) when running isprime.

When running performance tests with other operations than the isprime function - for example calculating eigenvalues of large random matrices - MatLab is a factor 2 to 10 faster than Octave. Therefor I am even more confused, why "isprime" is so much slower in MatLab...