/****************************************************************************/
/* Busca factores de un número de Mersenne:                                 */
/*  Si n es  un entero  positivo, entonces Mn=2^n -1  se                    */
/*  llama el n-simo número de Mersenne. Si p es un primo                    */
/*  impar,  cualquier divisor del  número Mp debe ser de                    */
/*  la forma 2kp+1 donde k es un entero positivo.                           */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://elparaiso.mat.uned.es                                            */
/****************************************************************************/

#include <gmp.h>
#include <stdio.h>
#include <stdlib.h>
typedef unsigned long ulong;

int esPar(ulong);
int esPrimo(ulong);
int divisorMersenne(mpz_t Mp, ulong p, mpz_t div);

int main(int argc, char *argv[])
{
	int i;
	ulong p;
	mpz_t Mp, dos, tmp, div;

	if (argc<2) {
		printf("Uso: %s p1 {p2} {p3} ... \n",argv[0]);
		printf("Busca divisores de los números de Mersenne\n");
		printf("Mp1, Mp2,... donde pi es un primo impar.\n");
		return 1;
	}
	mpz_init(Mp);
	mpz_init(dos);
	mpz_init(tmp);
	mpz_init(div);
		
	mpz_set_ui(dos,2);

	for (i=1; i<argc; i++) {
		p=atol(argv[i]);
		if (esPar(p) || !esPrimo(p)) {
			printf("* %ld no es un primo impar.\n\n",p);
			continue;
		}
		mpz_pow_ui(tmp,dos,p);
		mpz_sub_ui(Mp,tmp,1);
		printf("M%ld = ",p);
		mpz_out_str(stdout,10,Mp);
		if (divisorMersenne(Mp,p,div)) {
			printf("\n es divisible por ");
			mpz_out_str(stdout,10,div);
		}
		else printf(" es primo.");
		
		printf("\n\n");
	}
	return 0;
}

int esPar(ulong n) { return !(n&1); }
int esPrimo(ulong n)
{
	ulong div;

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

/* Prueba con todos los divisores del tipo 2pk+1 */
/* para encontrar divisores de Mp                */
/* Devuelve 1 si tiene divisores, 0 en otro caso */
/* el divisor,cuando existe, se coloca en div    */
int divisorMersenne(mpz_t Mp, ulong p, mpz_t div)
{
	mpz_t d,tmp,k,dosp,raiz,resto;

	mpz_init(d);
	mpz_init(dosp);
	mpz_init(k);
	mpz_init(raiz);
	mpz_init(resto);
	mpz_init(tmp);

	mpz_set_ui(tmp,p);
	mpz_mul_ui(dosp,tmp,2);		/* dosp= 2.p	*/
	mpz_set_ui(k,1);		/* k= 1		*/
	mpz_sqrt(raiz,Mp);		/* raiz=[raiz cuadrada de Mp] */
	
	do {
		mpz_mul(tmp,dosp,k);	/* tmp=2p.k     	*/
		mpz_add_ui(d,tmp,1);	/* d= 2p.k + 1  	*/
		mpz_mod(resto,Mp,d);	/* resto = Mp mod d 	*/
		if (mpz_sgn(resto)==0) {
			mpz_set(div,d); /* cuando el resto es 0 */
			return 1;
		}
		mpz_add_ui(tmp,k,1);
		mpz_set(k,tmp);		/* k=k+1		*/
	} while (mpz_cmp(raiz,d)>0);

	return 0;
}