使用 Lucas-Lehmer 素数检验的梅森素数

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

让我们把这个数字推得更远，扔掉所有的浮点代码。虽然下一个 mersenne 素数将适合 64 位，但问题是 termN * termN 表达式会在模数控制它之前溢出。如果我们有真正的模幂运算，我们可能会避免这个问题。相反，我们将在计算该值时使用 GCC/clang 中的模拟 128 位类型：

#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;
}

isPrime() 中的 i * i 效率稍低，但由于指数素数很小，所以这并不重要。

输出

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

Answer 2

问题出在一行：

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

您将 termN 声明为 unsigned int（在您的环境中为 32 位），但该产品可能变得如此之大而无法用这种类型表示，并且由此产生的溢出导致循环发散。

解决方案是使用范围更大的类型，例如unsigned long long int（64 位）。

请参阅Ideone.com 的示例。

Answer 3

改变这个：

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。

---

我是如何开始调试的？

通过在循环开始时使用 print 语句，如下所示：

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 个！

Answer 4

termN * termN 可能整数溢出。一般来说，您应该将可能是非常大的数字的值表示为双精度值，并尽可能避免在不同的数字类型之间进行转换，尤其是在整数和浮点数之间。