I can comment on multiple things here. But you ask about making the code look cleaner, so I'll start with that
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)