Mersenne Prime Finder Help

This is the place for queries that don't fit in any of the other categories.

Mersenne Prime Finder Help

Postby joelsjet » Sat Sep 21, 2013 2:37 am

Hi All,

I have some very messy and un commented code here which checks every number form zero to see if it is a Mi Prime

Here is the code:
Code: Select all
check=0
n=0
n2=0
per2=0
while n>=0:
    n=n+1


    import math
    a=2
    p=(pow(2,n)-1)
    count=math.sqrt(p)

    while a<count:

        per=((a/count)*100)
        per=math.ceil(per*10)/10
        if per2!=per:
            print per,"%"

            f2=open("C:/Mi Primes Log.txt", "a");
            string3=str(per)
            f2.write(string3)
            f2.write("%")
            f2.write("\n");

        per2=per

        if p % a==0:
            check=0
            break

        else:
            a=a+1
            check=1
            if n2!=n:
                print "finding prime", n

                string2=str(n)
                f2.write("finding prime ")
                f2.write(string2)
                f2.write("\n")
            n2=n
            continue


    if check==1:
        print n,"is a Mi Prime!!"
        print "which becomes", p
        f=open("C:/Mi Primes2.txt", "a");
        string=str(p)
        string6=str(n)
        f.write(string6)
        f.write("  ")
        f.write(string)
        f.write("\n");
        f.close

        f2.write(string6)
        f2.write(" is a Mi Prime!!")
        f2.write("\n");
        f2.write("which becomes ")
        f2.write(string)
        f2.write("\n");
        f2.close


I know i know its an absolute mess.

If you run it you will see the txt files appear at your drives start point and they work fairly well, the log looks messy but does still work.

The problem is if i stop the script, it then restarts from zero, which is hopeless, so i need some help working out a way to make it restart were it left of.

There is also a bug that it doesn't find the first mi prime being 2 which i thick has to do with the variable count but i'm not sure.

Also just in-case anyone was interested this type of brute force attack is hugely inefficient and slow, for instance the last mi prime i found was the 9th which is described 2^61-1, or 2,305,843,009,213,693,951. Which took roughly 1,518,500,250 division operations to check it was truly a prime!
joelsjet
 
Posts: 2
Joined: Sat Sep 21, 2013 1:56 am

Re: Mersenne Prime Finder Help

Postby casevh » Sat Sep 21, 2013 5:37 am

joelsjet wrote:Hi All,

I have some very messy and un commented code here which checks every number form zero to see if it is a Mi Prime

Here is the code:
Code: Select all
check=0
n=0
n2=0
per2=0
while n>=0:
    n=n+1


    import math
    a=2
    p=(pow(2,n)-1)
    count=math.sqrt(p)

    while a<count:

        per=((a/count)*100)
        per=math.ceil(per*10)/10
        if per2!=per:
            print per,"%"

            f2=open("C:/Mi Primes Log.txt", "a");
            string3=str(per)
            f2.write(string3)
            f2.write("%")
            f2.write("\n");

        per2=per

        if p % a==0:
            check=0
            break

        else:
            a=a+1
            check=1
            if n2!=n:
                print "finding prime", n

                string2=str(n)
                f2.write("finding prime ")
                f2.write(string2)
                f2.write("\n")
            n2=n
            continue


    if check==1:
        print n,"is a Mi Prime!!"
        print "which becomes", p
        f=open("C:/Mi Primes2.txt", "a");
        string=str(p)
        string6=str(n)
        f.write(string6)
        f.write("  ")
        f.write(string)
        f.write("\n");
        f.close

        f2.write(string6)
        f2.write(" is a Mi Prime!!")
        f2.write("\n");
        f2.write("which becomes ")
        f2.write(string)
        f2.write("\n");
        f2.close


I know i know its an absolute mess.

If you run it you will see the txt files appear at your drives start point and they work fairly well, the log looks messy but does still work.

The problem is if i stop the script, it then restarts from zero, which is hopeless, so i need some help working out a way to make it restart were it left of.

There is also a bug that it doesn't find the first mi prime being 2 which i thick has to do with the variable count but i'm not sure.

Also just in-case anyone was interested this type of brute force attack is hugely inefficient and slow, for instance the last mi prime i found was the 9th which is described 2^61-1, or 2,305,843,009,213,693,951. Which took roughly 1,518,500,250 division operations to check it was truly a prime!


Hi,

To restart where you left off, I would write the variable n to disk each time it changed. Then at the start of your program, read the value from the file. If the file doesn't exist, start with 0. Below is some sample code that should illustrate the idea. It is not a complete program.

Code: Select all
# Try to read a value for n from a file. If the file
# doesn't exist, set n to 0.

try:
    n = int(open("checking.txt", "w").readline())
except FileNotFoundError:
    n = 0
   
# Other code...

# Whenever you start processing a new value for n,
# save it.

open("checking.txt", "r").write(str(n))


The reason 2 is not detected is due to a combination of the test while a<count: and that check is still 0 when you reach if check==1:. I hope that is enough of a hint.

Are you interested in suggestions for testing Mersenne primes faster? If so, let me know.

casevh
casevh
 
Posts: 70
Joined: Sat Feb 09, 2013 7:35 am

Re: Mersenne Prime Finder Help

Postby joelsjet » Sat Sep 21, 2013 6:50 am

Hi casevh,

Thanks, i don't know why i didn't think to do it like that, and i could even further that by saving 'a' so it starts dividing from that point.

Sure it would be great to get some help learning how to improve the speed. I looked at things like Lucas–Lehmer test however i decided i was in way over my head.
I realize now after mucking around that the currant method is much to slow i took out the math.ceil for a quick test run at prime 89 and it was increasing at the 7th decimal place! .000000x super slow.

As for the bug i'm not sure. i picked up what your talking about while debugging and looking at value of the variables, by the time you square root p to get count its to small and the while function at line 14 is skipped, but i'm not sure how to get around that. I might leave that alone for now if you have input on making it faster code as everything is probably going to change
joelsjet
 
Posts: 2
Joined: Sat Sep 21, 2013 1:56 am

Re: Mersenne Prime Finder Help

Postby casevh » Sun Sep 22, 2013 5:22 am

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
casevh
 
Posts: 70
Joined: Sat Feb 09, 2013 7:35 am


Return to General Coding Help

Who is online

Users browsing this forum: Baidu [Spider], snippsat, W3C [Linkcheck] and 4 guests