View Single Post
 2019-01-04, 09:29 #38 GP2     Sep 2003 50338 Posts 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 #include #include 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