c - 使用 Lucas-Lehmer 素性检验的 Mersenne 素数

Question

这是代码，其中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;
}

输出:

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

它只返回第一个 5 而不是 8 我不明白为什么。

---

更新:

它跳过了从 13 开始的所有索引。我怀疑错误出在 isMersenneLucasLehmer(unsigned int) 的最后几行中的某处。找了很久都没找到。

Answer 1

改变这个:

unsigned int termN = 4;

为此:

unsigned long int termN = 4;

主要是因为你后来做了 termN * termN 这很可能在类型为 unsigned int 时导致溢出。

输出:

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

---

如果你应该打印你的类型会很好:

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

所以将 %lu 更改为 %u。

---

我是如何开始调试的？

在循环开始时使用打印语句，如下所示:

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

你会看到这个:

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

这意味着您看不到 8 个梅森数，因为您要结束循环，您的函数会找到其中的 8 个!