```/*********************************************************
/
/ 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);
}

```