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

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!