/********************************************************* 
/ 
/ mersenne.c 
/ 
/ A program to implement a Lucas-Lehmer primality test 
/ for Mersenne primes, to find all Mersenne primes up to 
/ (2^31)-1, and their associated perfect numbers. 
/ 
/ Paul Soper 
/ 
/ June, 2014 
/ 
/********************************************************/ 
 
#include "stdio.h" 
#include "stdlib.h"  /* because we use malloc */ 
 
/*  long longpow(long x, long y) does integer  
    exponentiation, returning x^y */ 
 
long longpow(long x, long y) { 
  long i = 0;  /* a counter */ 
  long product = 1; 
  for (i = 0; i < y; ++i) { 
    product = product * x; 
  } 
  return (product); 
} 
 
int main() { 
  /* Mersenne primes take the form M = 2^p - 1, where 
     p is a prime number.  So to find all M up to p = 31, 
     we need a list of primes up to 31.  We could generate 
     the list by trial division or by a sieve, but for 
     simplicity here I am going to populate the list by 
     hard coding it */ 
 
  long* pr = (long*)malloc(sizeof(long)*12); 
  *pr = 2; 
  *(pr+1) = 3; 
  *(pr+2) = 5; 
  *(pr+3) = 7; 
  *(pr+4) = 11; 
  *(pr+5) = 13; 
  *(pr+6) = 17; 
  *(pr+7) = 19; 
  *(pr+8) = 23; 
  *(pr+9) = 29; 
  *(pr+10) = 31; 
 
  /* Now we will step through the list of primes, and  
     determine for which primes 2^p-1 is also prime */ 
 
  int i = 0; /* a counter */ 
  for (i = 0; i <= 10; ++i) { 
    /* The Lucas-Lehmer primality test for Mersenne 
       primes works like this: 
       Establish a sequence:  s(0) = 4, and 
       s(n) = s(n-1)^2 - 2 
       M=2^p-1 is prime only if s(p-2) mod M = 0 
 
       Note that, to keep s from growing too quickly, 
       we can apply the modular division at each 
       step rather than waiting for s(p-2)  
 
       The test only works for odd p, so we have to  
       put in a work around for p = 2 */ 
     
 
    int p = *(pr + i);  /* the prime number for the  
			   exponent in the candidate M */ 
 
    long M = longpow(2, p) - 1;  /* the candidate M */ 
    
    long s = 4; /* the seed for the sequence */ 
 
    int j; /* a counter */ 
     
    /* apply the Lucas-Lehmer test */ 
 
 
    for (j = 0; j < (p-2); ++j) { 
      s = ((s*s) - 2) % M; 
    } 
 
     /* If s == 0, M is prime.  M = 2^p - 1.  If M is 
       prime, then (2^p - 1)*(2^(p-1)) is a perfect 
       number.  Print p, M, and the associated perfect 
       number */ 
     
    /* This is the work around for p = 2.  The LLT  
       does not work unless p is odd.  The only even 
       prime is p = 2, so we'll hard code that option */ 
    
    if ((s == 0) || (p == 2)) { 
      long perf = M * longpow(2, p-1); 
      printf ("p = %d, M = %ld is prime, %ld is perfect\n",  
	      p, M, perf); 
    } 
  } 
 
  /* free the allocated array and end */ 
 
  free (pr); 
 
  return (0); 
}