View Single Post
Old 2019-01-04, 09:29   #38
GP2
 
GP2's Avatar
 
Sep 2003

A1416 Posts
Default

Here is a Mersenne prime test. I didn't write it. It uses the GMP library.

It prints a 64-bit hexadecimal residue: "0" for a Mersenne prime, non-zero otherwise.

I compiled on Linux and ran on a single Skylake core on AWS.

For 86243 it takes 11 seconds; for 216091 it takes 1 minute 38.5 seconds.

And yet this is much slower than mprime/Prime95.

Code:
/* Adapted from http://rosettacode.org/wiki/Lucas-Lehmer_test#GMP */

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

void lucas_lehmer(unsigned long p, mpz_t *V)
{
  mpz_t mp, t;
  unsigned long k;

  mpz_init_set_ui(t, p);

  mpz_init(mp);
  mpz_setbit(mp, p);
  mpz_sub_ui(mp, mp, 1);

  mpz_init_set_ui(*V, 4);
  for (k = 3; k <= p; k++) {
    mpz_mul(*V, *V, *V);
    mpz_sub_ui(*V, *V, 2);
    /* mpz_mod(*V, *V, mp) but more efficiently done given mod 2^p-1 */
    if (mpz_sgn(*V) < 0) mpz_add(*V, *V, mp);
    /* while (n > mp) { n = (n >> p) + (n & mp) } if (n==mp) n=0 */
    /* but in this case we can have at most one loop plus a carry */
    mpz_tdiv_r_2exp(t, *V, p);
    mpz_tdiv_q_2exp(*V, *V, p);
    mpz_add(*V, *V, t);
    while (mpz_cmp(*V, mp) >= 0) mpz_sub(*V, *V, mp);
  }

  mpz_clear(t); mpz_clear(mp);

  /* Residue is returned in V */
}

int main(int argc, char* argv[]) {
  mpz_t m64, V;

  /* We want to print a 64-bit residue */
  mpz_init(m64);
  mpz_setbit(m64, 64);
  mpz_sub_ui(m64, m64, 1);

  unsigned long p;
  if (argc >= 2) {
      p = strtoul(argv[1], 0, 10);
  }
  else {
      fprintf(stderr, "Provide one argument, a Mersenne exponent, for instance: 11213\n");
      return 1;
  }

  lucas_lehmer(p, &V);

  mpz_and(V, V, m64);
  mpz_out_str(stdout, 16, V);
  fprintf(stdout, "\n");

  mpz_clear(m64); mpz_clear(V);

  return 0;
}

Last fiddled with by GP2 on 2019-01-04 at 09:38
GP2 is offline   Reply With Quote