Mersenne primes using Lucas-Lehmer primality test

Question

Here is the code, where limit = 8:

#include <stdio.h>
#include <math.h> // pow(x, exp)

//----------------------------------------------------------

char isMersenneLucasLehmer(unsigned int prime)
{
    unsigned int i, termN = 4;
    unsigned long mersenne;
    unsigned int limit;
    int res;

    mersenne = (unsigned long) pow(2, (double)prime) - 1;
    if (prime % 2 == 0)
    {
        return prime == 2;
    }
    else 
    {
        res = (int) sqrt((double) prime);
        for (i = 3; i <= res; i += 2)
        {
            if (prime % i == 0)
            {
                return 0;  
            }
        }

        limit = prime - 2;
        for (i = 1; i <= limit; ++i)
        {
            termN = (termN * termN - 2) % mersenne;
        }
    }
    return termN == 0;
}
//----------------------------------------------------------

/*
    Function: findMersenneLucasLehmer()

*/
void findMersenneLucasLehmer(unsigned int limit)
{
    unsigned int i, current = 0;
    unsigned long mersenne, bitsInLong = 64;

    for (i = 2; i <= bitsInLong; i++)
    {
        if (current >= limit)
        {
            break;
        }

        if (isMersenneLucasLehmer(i))   
        {
            mersenne = (unsigned long) pow(2, (double)i) - 1;
            printf("current = %lu, mersenne = %lu, index = %u\n", current, mersenne, i);
            ++current;
        } 
    }
}
//----------------------------------------------------------

int main()
{
    unsigned int limit = 8;
    findMersenneLucasLehmer(limit);
    return 0;
}

Output:

current = 0, mersenne = 3, index = 2
current = 1, mersenne = 7, index = 3
current = 2, mersenne = 31, index = 5
current = 3, mersenne = 127, index = 7
current = 4, mersenne = 8191, index = 13

It is only returning the first 5 instead of 8 and I can't figure out why.

---

Update:

it is skipping all the indexes from 13 and on. I suspect that the error is in somewhere in the last lines of isMersenneLucasLehmer(unsigned int). I've been staring for too long and couldn't find it.

Answer 1

Likely integer overflow at termN * termN. In general you should represent values that could be very large numbers as doubles, and avoid casting between different numeric types, especially between integers and floats, whenever possible.

Answer 2

Change this:

unsigned int termN = 4;

to this:

unsigned long int termN = 4;

mostly because you later do termN * termN which is likely to cause an overflow when a type of unsigned int.

Output:

current = 0, mersenne = 3, index = 2
current = 1, mersenne = 7, index = 3
current = 2, mersenne = 31, index = 5
current = 3, mersenne = 127, index = 7
current = 4, mersenne = 8191, index = 13
current = 5, mersenne = 131071, index = 17
current = 6, mersenne = 524287, index = 19
current = 7, mersenne = 2147483647, index = 31

---

It would be nice to print your types as you ought to:

C02QT2UBFVH6-lm:~ gsamaras$ gcc -Wall main.c
main.c:58:67: warning: format specifies type 'unsigned long' but the argument has type 'unsigned int' [-Wformat]
            printf("current = %lu, mersenne = %lu, index = %u\n", current, mersenne, i);
                              ~~~                                 ^~~~~~~
                              %u

So change %lu to %u.

---

How did I start debugging?

By using a print statement in the start of your loop, like this:

for (i = 2; i <= bitsInLong; i++)
{
    printf("Loop i = %u, current = %u\n", i, current);
    ...

You will see this:

current = 4, mersenne = 8191, index = 13
Loop i = 14, current = 5
...
Loop i = 63, current = 5
Loop i = 64, current = 5

which means that you don't see 8 Mersenne numbers, because you are ending your loop, before your function fins 8 of them!

Answer 3

The issue is in the line:

termN = (termN * termN - 2) % mersenne;

You declared termN as an unsigned int (32 bit, in your environment) but that product can become so big not to be representable by this type and the resulting overflow cause the loop to diverge.

The solution is to use a type with a bigger range, like unsigned long long int (64 bit).

See the example at Ideone.com.

Answer 4

Let's push this one number further and toss all that floating point code. Although the next mersenne prime will fit into 64 bits, the issue is the termN * termN expression which will overflow before the modulus can reign it in. If we had true modular exponentiation, we might avoid this problem. Instead, we'll use the emulated 128 bit type in GCC/clang while computing that value:

#include <stdio.h>
#include <stdbool.h>

typedef unsigned __int128 uint128_t;

bool isPrime(unsigned number)
{
    if (number % 2 == 0)
    {
        return number == 2;
    }

    for (unsigned i = 3; i * i <= number; i += 2)
    {
        if (number % i == 0)
        {
            return false;
        }
    }

    return true;
}

bool isMersenneExponent(unsigned exponent)
{
    if (exponent == 2)
    {
        return true;
    }

    if (!isPrime(exponent))
    {
        return false;
    }

    unsigned long termN = 4, mersenne = (1L << exponent) - 1;

    unsigned limit = exponent - 1;

    for (unsigned i = 1; i < limit; i++)
    {
        termN = (((uint128_t) termN * termN) % mersenne) - 2;
    }

    return termN == 0;
}

void findMersennePrime(unsigned limit)
{
    unsigned bit_limit = sizeof(unsigned long) * 8;

    for (unsigned current = 0, i = 2; current < limit && i < bit_limit; i++)
    {
        if (isMersenneExponent(i)) 
        {
            unsigned long mersenne = (1L << i) - 1;
            printf("current = %u, mersenne = %lu, index = %u\n", current++, mersenne, i);
        }
    }
}

int main()
{
    unsigned limit = 9;
    findMersennePrime(limit);
    return 0;
}

The i * i in isPrime() is a slight inefficiency but since the exponent primes are tiny, it hardly matters.

OUTPUT

current = 0, mersenne = 3, index = 2
current = 1, mersenne = 7, index = 3
current = 2, mersenne = 31, index = 5
current = 3, mersenne = 127, index = 7
current = 4, mersenne = 8191, index = 13
current = 5, mersenne = 131071, index = 17
current = 6, mersenne = 524287, index = 19
current = 7, mersenne = 2147483647, index = 31
current = 8, mersenne = 2305843009213693951, index = 61