/****************************************************************************/
/* Primos de Wilson (p-1)! = -1 (mod p^2)                                   */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://elparaiso.mat.uned.es                                            */
/****************************************************************************/

#include <stdio.h>
#include "lip.h"

void uso(char *progname);
void facmod(verylong p, verylong m, verylong *res);
int esprimoWilson(verylong p);

int main(int argc, char *argv[])
{
	verylong min=0, max=0, p=0, tmp=0, uno=0;

	if (argc==1 || argc>3) {
		uso(argv[0]); 
		return 1;
	}
	zsread(argv[1],&min);
	if (argc==3) zsread(argv[2],&max);
	
	zcopy(min,&p);
	zone(&uno);
	if (zcompare(max,0)) {
		while (zcompare(p,max)<0) {
			if (esprimoWilson(p)) zwriteln(p);
			zadd(p,uno,&tmp);			
			zcopy(tmp,&p);
		}
	}
	else {
		while (1) {
			if (esprimoWilson(p)) {zwriteln(p); break;}
			zadd(p,uno,&tmp);			
			zcopy(tmp,&p);
		}
	}

	return 0;
}

/*
 * 	Devuelve 1 si p es un primo de Wilson, 0 en otro caso
 */

int esprimoWilson(verylong p)
{
	verylong pcuad=0,uno=0,muno=0,tmp=0,res=0;

	zone(&uno);
	if (zprobprime(p,100)==0) return 0;	/* no es primo */
	zsq(p,&pcuad);
	zsub(pcuad,uno,&muno);
	zsub(p,uno,&tmp);
	facmod(tmp,pcuad,&res);			/* res= (p-1)! mod p^2 */
	if (zcompare(res,muno)==0) return 1;	/* ¿res=-1 mod p^2? */
	return 0;
}

/*
 * 	Calcula  res = n! mod m 
 */

void facmod(verylong n, verylong m, verylong *res)
{
	verylong uno=0,tmp=0;

	zone(res); 
	zone(&uno);       	   	/* res=1 uno=1 */
	while(zcompare(n,uno)) {
		zmul(*res,n,&tmp);	/* res= res*n */
		zcopy(tmp,res);
		zsub(n,uno,&tmp);	/* n=n-1        */
		zcopy(tmp,&n);
	}
	zmod(*res,m,&tmp);
	zcopy(tmp,res);
	
	return;
}

/*
 * 	Instrucciones de uso
 */

void uso(char *progname)
{
	fprintf(stderr,"-------\n");
	fprintf(stderr,
	 "%s <n> busca el primer primo de Wilson > n.\n",progname);
	fprintf(stderr,
	 "%s <n> <m> todos los primos de Wilson entre n y m.\n",progname);
	fprintf(stderr,"-------\n");

	return;
}