Code for finding Mersenne Primes

I am trying to answer the following question and am failing miserably at it.

A Mersenne prime is a prime number that is equal to 2^n-1, where n is an integer. For exampke, 31 is a Mersenne prime since 31=2^5-1. Write a program that finds all the Mersenne primes between 1 and 10,000. Do not use MATLAB's built-in function isprime.

Here's what I have so far:

clc,clear
n=10000;
prime=[1 2 3];
s=0;
for i=4:n
    isprime=1;
    for j=2:i-1 %(or i/2) because i/2 is optimal
        if rem(i,j)==0
            isprime=0;
        end
    end
    if isprime==1
        prime(end+1)=i;
    end
end
k=length(prime);
n=1;
l=0;
MersennePrimeIndex=0;
for i=2:k
    if (2^n)-1==find(prime(i));
        MersennePrimeIndex=MersennePrimeIndex+1;
        l(MersennePrimeIndex)=prime(i);
        n=n+1;
    end
end

Answer 1

[Walter Roberson]

Your line

      if (2^n)-1==find(prime(i));

is wrong. prime(i) is going to be a single non-zero value, and find() on a scalar non-zero value is going to return 1.

I suggest you change your test: take each prime, add 1, and determine whether the result is a power of 2.

Alternately, just go through the powers of 2 that are in the designated range, subtract 1, and test to see if the values are present in the list of primes you constructed.

[Walter Roberson]

Hint: factor() helps.

[John D'Errico]

Funny. I realized why I mistook the goal of this question. 31=2^5-1 is a Mersenne prime. But Mersenne primes are usually defined simply by the exponent of 2. And since 2^31-1 is also a Mersenne prime, I saw the 31 there, and jumped to the wrong conclusion. So I should have more carefully read the question.