5
\$\begingroup\$

I have coded the Lucas-Lehmer primality test following Wikipedia's description. I used the mod 2^n - 1 suggested in the article but was wondering if there were any other improvements I could make. I am using the GMP library for arbitrary precision integers. The program asks an integer n as input, then checks every Mersenne number M(i) for all i less than n. n = 10000 takes around 25 seconds on my computer.

#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;
}
\$\endgroup\$
1
\$\begingroup\$

The GMP code on RosettaCode is about 2x faster if you're looking for performance improvements.

It has a reasonable number of comments explaining the various optimizations. Most of the help will be in the pre-tests, where we know easy ways to find out if the result will be composite. The test itself is not too dissimilar, with a few different choices made.

I don't particularly like your for loop, which you wrote more like a while loop -- mixing the iteration and test together.

\$\endgroup\$
  • \$\begingroup\$ While this answer is welcome, it would be preferable to include the essential parts of the linked resource, so that the answer is complete in itself, and still valuable when the link target changes or disappears. I suggest you edit to incorporate relevant suggestions (with samples if you can). Also, could you you be more specific with your final comment? Perhaps you can show a clearer way to write the loop. \$\endgroup\$ – Toby Speight Aug 9 '17 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.