Calculate Mersenne prime numbers in MATLAB

A Mersenne prime is equal to (2^n - 1), where n is an integer. 31 is a Mersenne prime since 31 = 2^5 -1. I need to write a code that calculates Mersenne primes between 1 and 10,000. Lastly I'm not allowed to use the build in function isprime.

I have used the formula to calculate till an arbitrary n then used a for loop to verify if the numbers I got are prime numbers.

Running the script I get 3, 7, 31, 127, 8191.

How could I make the code look cleaner? Do you think my code is giving me the right answer?

    for k = 1:50000 %calculates mersenne primes using mersenne formula

p(k) = 2^k -1;

end
p = p(p<10000 & p>1);
psize =size(p);
x = psize(2);
mersenneprime = 0; k = 2; f = 0;

for m = 1:x %verifies if mersenne are prime numbers
k = k +1;
a = p(m);
b = p(m)-1;
if all((rem(a,2:b)) ~= 0)
flag = 1;
else
flag = 0;
end
if flag == 1
f = f+1;
mersenneprime(f) = p(m);
end
end %it eliminates 7 numbers from original array
mersenneprime


Indenting

Be consistent with indenting. Indents show the structure of your code. Inconsistent indenting makes it hard to follow the structure.

Consistent spacing around operators also helps readability.

Variable names

It's common in MATLAB code to use single-letter variable names. I do this too. Try to use variable names that make the purpose of the variable clear. For example, p could be called candidates, since they are the possible Mersenne prime candidates.

Variable usage

psize = size(p);
x = psize(2);
% <snip>
for m = 1:x


You define a variable x here, but then use it only once. The name x is meaningless, so it takes a bit of time to find out what you are looping over. Consider instead:

for m = 1:length(p)


That one statement is immediately clear.

But since you are only using m to index into p (a = p(m)), you can also do this:

for a = p


Now m will have one of the values of p in each loop iteration (you don't need to do the indexing). Note that this form of for iterates over the columns of p, so if p were a column vector instead of a row vector, then the loop would only run once, with a containing all values in p.

You actually index p(m) three times within this loop. Two of those you can replace with a, which you have assigned the value of p(m).

Finally, k is initialized, and incremented each loop iteration, but not used at all. Unnecessary assignments and operations add confusion to your code, try to avoid these.

Vectorize

There are several places where you can simplify your code.

for k = 1:50000 %calculates mersenne primes using mersenne formula
p(k) = 2^k -1;
end


This is identical to

p = 2.^(1:50000) - 1;


(note the .^ element-wise operator). Here you generate 50 thousand powers of two, but next you keep only 12:

p = p(p<10000 & p>1);


Why not directly generate the values you intend to keep? For the condition 2.^k-1 < 10000, you need k < log2(10000-1), and for 2.^k-1 > 1 you need k > 1:

p = 2.^(2:log2(10000-1))-1;


Testing for primeness

Next, you loop over each of the candidates to see if it's a prime:

a = p(m);
b = p(m)-1;
if all((rem(a,2:b)) ~= 0)


Note that to test for primeness, you need to test division only up to sqrt(a), any larger value is guaranteed not to divide the value evenly. And you don't need to test with any power of 2 except 2 itself:

if all(rem(a,[2,3:2:sqrt(a)]) ~= 0)


You could consider making this test into a function, which will improve readability (then you don't need to add the comment):

isprime = @(a)all(rem(a,[2,3:2:sqrt(a)]) ~= 0);
% ...
if isprime(a)


(now you're using the isprime function, but not the built-in one!)

Simplify logic

Within this test you set flag = 1 or flag = 0. Considering true has a value of 1, and false a value of 0, you might as well do:

flag = isprime(a);


But next you test flag and add the number to the list of Mersenne primes. So you don't need the flag at all:

if isprime(a)
f = f+1;
mersenneprime(f) = a;
end


Appending to an array

Adding a value to the mersenneprime list can be simplified too, using the end keyword in indexing, do you don't need f:

mersenneprime(end+1) = a;


This is actually fairly efficient. They have tried to make people preallocate arrays like these that grow within a loop for decades, and finally decided to optimize their allocation strategy for this use case. Also, since your list will always be very small, no matter how large your upper limit is, you don't need to worry about efficiency here at all. You do need to initialize the mersenneprime to an empty array.

Putting it all together

So putting this all together I get to the following bit of code:

isprime = @(a)all(rem(a,[2,3:2:sqrt(a)]) ~= 0);
limit = 10000;
candidates = 2.^(2:log2(limit-1))-1;
mersenneprime = [];
for c = candidates
if isprime(c)
mersenneprime(end+1) = c;
end
end
disp(mersenneprime)

• Great, comprehensive answer! I feel like the isprime function must be fully vectorisable, i.e. you should be able to bin off the candidates loop? – Wolfie Oct 13 '18 at 15:29
• Thanks @Wolfie! It does look like it should be vectorizable, but I don’t immediately see how. Because the upper limit is different for each input, bsxfun would not work (or do lots of useless computation). – Cris Luengo Oct 13 '18 at 16:05