使用 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

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

Answer 2

改变这个：

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

Answer 3

问题在行中：

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

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

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

参见 Ideone.com 中的示例。

Answer 4

让我们进一步推动这个数字并抛出所有浮点代码。虽然下一个 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

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

Mersenne primes using Lucas-Lehmer primality test

c

algorithm

primes

scientific-computing

primality-test

更新：