Hi,

For the rest of my reply, I'll use ^ to mean exponentiation and let M(n) = 2^n - 1.

You are attempting to prove M(n) is prime by trial factoring. Unfortunately, trial factoring slows down very quickly. It is still useful for small number.

Your trial factoring code starts at 2 and then increments the trial divisor by 1 each time through the loop. This does work but performs many extra divisions. There is no need to try dividing by 4 since you would have already found 2 as a factor. You can double the speed just by dividing by 2 and then by dividing by all the odd numbers. With a little extra work, you can skip all the multiples of 3 and 5.

Here is an implementation that skips multiples of 2, 3, and 5. I included a limit to the range of factors to try - if you exceed this limit, you probably should be trying a different approach.

On my system, it verified that M(61) is prime in around 50 seconds.

- Code: Select all
`def trial_factor(n):`

'''Return the first non-trivial factor of n. Return n if n is

prime. Trial factoring is done up to 10,000,000,000. If no factor

is found below this limit, 0 is returned.'''

if n < 2:

raise ValueError('n must be greater than or equal to 2')

# Include a short list of primes so we can efficiently test

# for small factors. The primes must be listed in order and the

# last prime listed should be the last prime below a multiple of

# 30. For example, if you want to extend the list, the list could

# to to 59 (the last prime before 60) or to 89 (the last prime

# before 90).

primes = [2, 3, 5, 7, 9, 11, 13, 17, 19, 23, 29]

# Check small numbers very quickly.

if n in primes:

return n

# Test with all the values in our short list of primes.

for tf in primes:

if not n % tf:

return tf

# If a number is composite, the smallest factor must be less than

# the square root of number. Conversely, if we have tried trial

# factoring up to a limit K and n is less than K^2, then n must be

# prime.

if n < tf ** 2:

return n

# Begin trial factoring with divisors beyond our short list of

# primes. Since we require the last prime in a list be the last

# prime before a multiple of 30, we can start our trial factoring

# at the first odd number in the next multiple of 30.

tf = 30 * (tf // 30 + 1) - 1

# Skip multiples of 2, 3, and 5. Only 8 trial divisions are needed

# to increment the trial factor by 30.

skip30 = [2, 6, 4, 2, 4, 2, 4, 6]

while tf * tf < n:

for k in skip30:

tf = tf + k

if not n % tf:

return tf

# Enforce a reasonable upper limit.

if tf > 10000000000:

return 0

# Trial factoring is complete so n must be prime.

return n

The trial_factor() function is useful in general. It can be used to fully factor a reasonably sized number or verify that small numbers are prime. As an exercise, try writing a couple of functions:

1) use trial_factor() to fully factor a small number (it should return a list of all the prime factors)

2) use trial_factor() to return True (False) if a small number is prime (composite)

It can be proven that the exponent of a prime Mersenne number must also be a prime number. You can decrease the numbers you need to test by first checking if the exponent is prime and only then proceeding to trial factoring.

Can we make other improvements to trial factoring? If we know we are trying to factor a candidate for a Mersenne prime (the exponent p is prime), then any factor must be of the form 2*p*K + 1 for some K. Here is function that uses this fact to speed up trial division. Note that only the exponent is passed to the function.

- Code: Select all
`def mersenne_trial_test(p):`

'''Attempt to prove that 2**p - 1 is prime by trial factoring

and using the fact that any divisor must be of the form

2 * P * K + 1.'''

n = 2 ** p - 1

delta = 2 * p

tf = delta + 1

while tf * tf < n:

if not n % tf:

return False

tf = tf + delta

return True

mersenne_trial_test(61) completes in about 2 seconds on my computer. But mersenne_trial_test(89) takes a very long time. If it completes, I'll update the running time.

There are better algorithms for factoring numbers but factoring numbers larger than 200 digits in length is a very difficult problem. However, there are fast algorithms that can determine if a number is most likely a prime number. And for certain classes of numbers, the algorithms can give guaranteed results. The Lucas-Lehmer test is a very efficient test for Mersenne primes.

The Lucas-Lehmer test generates a sequence beginning with S_0 = 4. The next terms in the sequence are generated by the following rule: S_n = ((S_n-1 * S_n-1) - 2) mod M(p). Then M(p) is prime if S_p-2 is 0. Here is an implementation of Lucas-Lehmer test.

- Code: Select all
`def lucas_lehmer_test(p):`

'''Attempt to prove that 2**p - 1 is prime with Lucas-Lehmer

test. It is assumed that p is an odd prime.'''

s = 4

n = 2 ** p - 1

# Use range(p - 2) for Python 3.x

for i in xrange(p - 2):

s = (s * s - 2) % n

return not s

The Lucas-Lehmer test verifies M(61) is less than a thousandth of a second. M(89) only takes a couple thousandths of a second. M(44497) is verified prime in less than 3 minutes.

I hope you find these examples helpful.

casevh