/****************************************************************************/
/* Test de Lucas-Lehmer, comprueba la primalidad de los números de Mersenne.*/
/* Rosen, Elementary number theory and its applications ISBN 0-201-11958-7  */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://elparaiso.mat.uned.es                                            */
/****************************************************************************/

#include <stdio.h>

#include "gmp.h"


int LucasLehmer(unsigned long p);
int esPrimo(unsigned long n);

int main()
{
	unsigned long p;
	
	printf("Decide si Mp, el p-simo número de Mersenne es primo.\n");
	printf("p= ");
	scanf("%ld",&p);
	
	printf("M%ld ",p);
	if (LucasLehmer(p)) printf("es primo.\n");
	else printf("es compuesto.\n");
	return 0;
}

/*
 * Test de Lucas-Lehmer: devuelve 1 si Mp= 2^p-1 es primo y 0 en otro caso
 * Si R(1)=4 y R(k) =R(k-1)^2 -2 (mod Mp) entonces 
 * Mp es primo si y solo si R(p-1)=0 (mod Mp)
 * - Coste del algoritmo: O(p^3)
 */
int LucasLehmer(unsigned long p)
{
	mpz_t r,Mp;
	int escero;
	unsigned long i;
	
	/* Si p no es primo Mp no puede serlo */
	if (!esPrimo(p)) return 0;
	
	mpz_init(r); mpz_init(Mp);

	/* Mp= 2^p -1 */
	mpz_ui_pow_ui(Mp,2,p);
	mpz_sub_ui(Mp,Mp,1);
	
	/* Comienzo del ciclo principal */
	mpz_set_ui(r,(unsigned long)4);		/* r = 4              */
	for (i=2; i<p; i++) {
		mpz_mul(r, r,r);		/* r' = r*r           */
		mpz_sub_ui(r,r,2);		/* r' = r*r -2        */ 
		mpz_mod(r,r,Mp);		/* r' = r*r -2 mod Mp */
	}

	escero=(mpz_sgn(r)==0);			/* ¿ r= 0 ? */
	mpz_clear(r); mpz_clear(Mp);
	return escero;
}


/* Un test por divisiones para ver si un entero sin signo es primo */
int esPrimo(unsigned long n)
{
	unsigned long d;

	if (n%2 ==0) return 0;
	for (d=3; d*d<=n; d+=2)
		if (n%d ==0) return 0;
	return 1;
}