/****************************************************************************/
/* Solución de un sistema de congruencias x= ai (mod mi) i=1...r            */
/* por medio de la función fi de Euler.                                     */
/*                                                                          */
/* Se demuestra que las soluciones vienen dadas por                         */
/* x = a1.M1^fi(m1)+...+ ar.Mr^fi(mr) (mod M)                               */
/* siendo M= m1.m2...mr , Mi=M /mi                                          */
/*                                                                          */
/* 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 <stdlib.h>


#define PAR(x) !(x&1)

#define NOMEM 1

#define NOPRIM 2


typedef unsigned long ulong;

int scngeuler(int n, ulong *a, ulong *m, ulong *x);
int primos_a_pares(int n,ulong *m);
ulong mcd(ulong a, ulong b);
ulong fi(ulong n);

int main()
{
	ulong *m,*a,x,M;
	int i,n,ret;

	puts("Este programa resuelve sistemas de congruencias.");
	puts("del tipo x=a1 mod m1;... x=an mod mn.");
	puts("siendo los mi primos entre si dos a dos.\n");
	
	printf("Número de ecuaciones n= "); scanf("%d",&n);
	m=(ulong *)malloc(n*sizeof(ulong));
	a=(ulong *)malloc(n*sizeof(ulong));
	if (a==NULL) {
		puts("No hay suficiente memoria.");
		return NOMEM;
	}
	for (i=0; i<n; i++) {
		printf("\na%d  =",i+1); scanf("%ld",a+i);
		printf("m%d  =",i+1);   scanf("%ld",m+i);
	}
	ret=scngeuler(n,a,m,&x);
	
	if (ret==NOPRIM) {
		puts("\nLos modulos no son primos entre si dos a dos.");
		return NOPRIM;
	} else if (ret==NOMEM) {
		puts("\nNo hay suficiente memoria.\n");
		return NOMEM;
	}
	
	for (M=1L,i=0; i<n; i++) M*= *(m+i);
	printf("\nLa solución es x= %ld (mod %ld)\n",x,M);
	free(m); free(a);
	return 0;
}

/*
 * Funcion scngeuler()
 * Descrip: resuelve un sistema de ecuaciones x=a1 mod m1... x=an mod mn
 *          siempre y cuando los modulos sean primos entre si dos a dos.
 * Entrada: n numero de ecuaciones,a coeficientes de las ecuaciones
 *          m modulos de las ecuaciones, x resultado del sistema
 * Salida : O      --> Ok 
 * 	    NOMEM  --> No hay suficiente memoria
 * 	    NOPRIM --> Los mi no son primos entre si dos a dos
 * Efecto : si devuelve 0 la solucion del sistema esta en *x,
 *          en otro caso el valor de x no se modificara.
 *          
 * ¡ATENCION! no se hacen comprobaciones de overflow, que efectivamente 
 * se produce si hay suficiente número de ecuaciones o los módulos "mi"
 * son suficientemente grandes. En estos casos la solución sera erronea.
 * 
 */

int scngeuler(int n, ulong *a, ulong *m, ulong *x)
{
	ulong *M,*y,K=1;
	int i,j;

	M=(ulong *)malloc(n*sizeof(ulong));
	y=(ulong *)malloc(n*sizeof(ulong));
	if (y==NULL) return NOMEM;
	if (primos_a_pares(n,m) != 1) return NOPRIM;
	
	/* K=m1.m2...mn */
	for (i=0; i<n; i++) K*= *(m+i);
	/* Mi=K/mi   yi= Mi^fi(mi)  */
	for (i=0; i<n; i++) {
		*(M+i)= K / *(m+i);
		for (*(y+i)=1, j=0; j<fi(*(m+i)); j++) 
			*(y+i)= (*(y+i)* *(M+i))% K;
	}
	
	/* x=a1.y1+a2.y2+...+an.yn (mod K) */
	for (*x=0,i=0; i<n; i++)
		*x +=  (*(a+i) * *(y+i))%K;
	*x %= K;
	free(M);
	free(y);
	return 0;
}

/*
 * funcion primos_a_pares()
 * Averigua si una serie de enteros son primos dos a dos
 */

int primos_a_pares(int n, ulong *m)
{
	int i,j;

	for (i=0; i<n-1; i++)
		for (j=i+1; j<n; j++) {
			if (mcd(*(m+i),*(m+j))!=1) return 0;
		}
	return 1;
}

/*
 * funcion mcd()
 * Descrip: devuelve el mcd de dos enteros sin usar divisiones,
 *          solo restas y desplazamientos.
 */

ulong mcd(ulong a, ulong b)
{
	if (a==b) return a;
	else if (a<b) return mcd(b,a);
	else if ( PAR(a) &&  PAR(b))  return 2*mcd(a>>1,b>>1);
	else if ( PAR(a) && !PAR(b))  return mcd(a>>1,b);
	else if (!PAR(a) &&  PAR(b))  return mcd(a,b>>1);
	else                          return mcd(a-b,b);
}

/*
 * Función fi de Euler, número de enteros menores que uno dado
 * y primos con él. Suponiendo que n = p1^e1 . p2^e2 ... pk^ek 
 * entonces   fi(n)= n.(1-1/p1).(1-1/p2)...(1-1/pk)
 */

ulong fi(ulong n)
{
	ulong p;
	float f;

	f=n;
	if ((float)n/2 == n/2) {
		f *= (1-1.0/2);
		while ((float)n/2 == n/2) n/=2;
	}
	
	p=3;
	while (p*p <= n) {
		if ((float)n/p == n/p) {
			f *= (1-1.0/p);
			while ((float)n/p == n/p) n/=p;
		}
		p += 2;
	}
	if (n!=1) f *= (1-1.0/n);

	return (ulong)f;
}