Guest User

mersenne.cpp

a guest
Feb 9th, 2018
50
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <limits.h>
  4. #include <gmp.h>
  5. /*g++ mersenne.c -o mersenne.o -L/gmp_install/lib -lgmp */
  6. int lucas_lehmer(unsigned long p)
  7. {
  8.   mpz_t V, mp, t;
  9.   unsigned long k, tlim;
  10.   int res;
  11.  
  12.   if (p == 2) return 1;
  13.   if (!(p&1)) return 0;
  14.  
  15.   mpz_init_set_ui(t, p);
  16.   if (!mpz_probab_prime_p(t, 25)) /* if p is composite, 2^p-1 is not prime */
  17.     { mpz_clear(t); return 0; }
  18.  
  19.   if (p < 23)                     /* trust the PRP test for these values */
  20.     { mpz_clear(t); return (p != 11); }
  21.  
  22.   mpz_init(mp);
  23.   mpz_setbit(mp, p);
  24.   mpz_sub_ui(mp, mp, 1);
  25.  
  26.   /* If p=3 mod 4 and p,2p+1 both prime, then 2p+1 | 2^p-1.  Cheap test. */
  27.   if (p > 3 && p % 4 == 3) {
  28.     mpz_mul_ui(t, t, 2);
  29.     mpz_add_ui(t, t, 1);
  30.     if (mpz_probab_prime_p(t,25) && mpz_divisible_p(mp, t))
  31.       { mpz_clear(mp); mpz_clear(t); return 0; }
  32.   }
  33.  
  34.   /* Do a little trial division first.  Saves quite a bit of time. */
  35.   tlim = p/2;
  36.   if (tlim > (ULONG_MAX/(2*p)))
  37.     tlim = ULONG_MAX/(2*p);
  38.   for (k = 1; k < tlim; k++) {
  39.     unsigned long q = 2*p*k+1;
  40.     /* factor must be 1 or 7 mod 8 and a prime */
  41.     if ( (q%8==1 || q%8==7) &&
  42.          q % 3 && q % 5 && q % 7 &&
  43.          mpz_divisible_ui_p(mp, q) )
  44.       { mpz_clear(mp); mpz_clear(t); return 0; }
  45.   }
  46.  
  47.   mpz_init_set_ui(V, 4);
  48.   for (k = 3; k <= p; k++) {
  49.     mpz_mul(V, V, V);
  50.     mpz_sub_ui(V, V, 2);
  51.     /* mpz_mod(V, V, mp) but more efficiently done given mod 2^p-1 */
  52.     if (mpz_sgn(V) < 0) mpz_add(V, V, mp);
  53.     /* while (n > mp) { n = (n >> p) + (n & mp) } if (n==mp) n=0 */
  54.     /* but in this case we can have at most one loop plus a carry */
  55.     mpz_tdiv_r_2exp(t, V, p);
  56.     mpz_tdiv_q_2exp(V, V, p);
  57.     mpz_add(V, V, t);
  58.     while (mpz_cmp(V, mp) >= 0) mpz_sub(V, V, mp);
  59.   }
  60.   res = !mpz_sgn(V);
  61.   mpz_clear(t); mpz_clear(mp); mpz_clear(V);
  62.   return res;
  63. }
  64.  
  65. int main(int argc, char* argv[]) {
  66.   unsigned long i, n = 2147483647;
  67.   if (argc >= 2) n = strtoul(argv[1], 0, 10);
  68.   for (i = 2 ; i <= n; i++) {
  69.     if (lucas_lehmer(i)) {
  70.       printf("M%lu ", i);
  71.       fflush(stdout);
  72.     }
  73.   }
  74.   printf("\n");
  75.   return 0;
  76. }
RAW Paste Data