卢卡斯·莱默素数检验

Question

我按照维基百科的描述对卢卡斯-莱默素数检验进行了编码。我使用了文章中建议的mod 2^n - 1，但我想知道是否还有其他改进可以改进。我正在使用GMP库来处理任意精度的整数。该程序要求一个整数n作为输入，然后检查每个Mersenne数字M(i)的所有i小于n。n = 10000在我的电脑上大约需要25秒。

#include <gmp.h>
#include <stdio.h>

_Bool mersenne_prime(unsigned long exponent) {
    mpz_t sequence, number, temp;
    // Intialize variables
    mpz_init_set_ui(sequence, 4);
    mpz_init(temp);
    // Set number to 2^n - 1
    mpz_init_set_ui(number, 1);
    mpz_ui_pow_ui(number, 2, exponent);
    mpz_sub_ui(number, number, 1);
    // Repeat exponent-2 times
    for (unsigned long counter = exponent; --counter - 1;) {
        mpz_mul(sequence, sequence, sequence);
        // Modulus suggested by wikipedia
        while (mpz_cmp(sequence, number) > 0) {
            // Most significant bits of sequence
            mpz_div_2exp(temp, sequence, exponent);
            // Least significant bits of sequence
            mpz_mod_2exp(sequence, sequence, exponent);
            mpz_add(sequence, sequence, temp);
        }
        // Erratic case
        if (mpz_cmp(number, sequence) == 0) {
          mpz_set_ui(sequence, 0);
        }
        mpz_sub_ui(sequence, sequence, 2);
    }
    // sequence == 0 means prime
    _Bool result = mpz_sgn(sequence) == 0;
    // Clear variables
    mpz_clears(sequence, number, temp, NULL);
    return result;
}

int main() {
    mpz_t num;
    unsigned long limit;
    mpz_init_set_ui(num, 3);
    // Get limit
    scanf("%lu", &limit);
    printf("Searching for Mersenne primes...\nM2 is prime!\n");
    while (mpz_cmp_ui(num, limit) < 0) {
        unsigned long num_ui = mpz_get_ui(num);
        if (mersenne_prime(num_ui)) {
            printf("M%lu is prime!\n", num_ui);
        }
        mpz_nextprime(num, num);
    }
    // Clear num
    mpz_clear(num);
    return 0;
}

Answer 1

如果您正在寻找性能改进，RosettaCode上的GMP代码大约快2倍。

它有合理数量的注释解释了各种优化。大部分的帮助将在预测试中，在那里我们知道简单的方法来找出结果是否是合成的。测试本身并没有太大的不同，做了一些不同的选择。

我不特别喜欢您的for循环，您编写的这个循环更像是while循环--将迭代和测试混合在一起。