/****************************************************************************/
/* Teorema chino del resto                                                  */
/* Ref: Rosen "Elementary number theory..."                                 */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://www.matematicas.net                                            */
/****************************************************************************/

#include <stdio.h>

#include <stdlib.h>

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

int tchino(int n, long *a, long *m, long *x);
int primos_a_pares(int n, long *m);
long mcd(long a, long b);
int invmod(long a, long m, long *inv);
void mcd_aux(long a, long b, long *alfa, long *beta, long *m);

int main()
{
	long *m,*a,x;
	int i,n;

	puts("Este programa resuelve sistemas de congruencias.");
	puts("del tipo x=a1 mod m1;... x=an mod mn.\n");
	printf("Número de ecuaciones n= "); scanf("%d",&n);
	m=(long *)malloc(n*sizeof(long));
	a=(long *)malloc(n*sizeof(long));
	if (a==NULL) {
		puts("No hay suficiente memoria.");
		return 1;
	}
	for (i=0; i<n; i++) {
		printf("\na%d  =",i+1); scanf("%ld",a+i);
		printf("m%d  =",i+1);   scanf("%ld",m+i);
	}
	if (tchino(n,a,m,&x)) {
	    puts("\nNo puedo solucionar ese sistema.");
	    return 1;
	}
	printf("\nLa solución es %ld.\n",x);
	return 0;
}

/*
 * Funcion tchino()
 * 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 1 --> Error, no se pudo solucionar el sistema
 * Efecto : si devuelve 0 la solucion del sistema esta en *x,
 *          en otro caso el valor de x no se modificara.
 */

int tchino(int n, long *a, long *m, long *x)
{
	long *M,*y,K=1;
	int i;

	M=(long *)malloc(n*sizeof(long));
	y=(long *)malloc(n*sizeof(long));
	if (y==NULL) return 1;
	if (primos_a_pares(n,m) != 1) return 1;
	for (i=0; i<n; i++) K*= *(m+i);
	for (i=0; i<n; i++) {
		*(M+i)= K / *(m+i);
		invmod(*(M+i),*(m+i),y+i);
	}
	*x= 0;
	for (i=0; i<n; i++)
		*x +=  *(a+i) * *(M+i) * *(y+i);
	*x %= K;
	return 0;
}

/*
 * funcion primos_a_pares()
 * Descrip: averigua si una serie de enteros son primos dos a dos
 * Entrada: n numero de enteros, m puntero a la lista de enteros
 * Salida : 0 no son primos entre si, 1 primos entre si de dos en dos.
 */

int primos_a_pares(int n, long *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.
 */

long mcd(long a, long 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);
}

/*
 * funcion invmod()
 * Descrip: devuelve el inverso de un entero modulo m si es que existe
 * Entrada: a un entero, m el modulo, inv contendra el resultado
 * Salida : 0-->Ok, 1-->Error, no inversible
 * Efecto : *inv contiene el inverso de a modulo m cuando sea posible
 */

int invmod(long a, long m, long *inv)
{
	long alfa, beta, r;

	mcd_aux(a,m,&alfa,&beta,&r);
	if (r!=1) return 1;
	*inv = alfa;
	if (*inv<0) *inv += m;
	return 0;
}

/* 
 * funcion mcd_aux()
 * Descrip: calcula el mcd de dos enteros y los coeficientes que
 *          permiten expresar dicho mcd como combinacion lineal
 *          de dichos enteros.
 * Entrada: a y b enteros long, alfa, beta y m punteros a long
 * Salida : ninguna
 * Efecto : *m contendra el mcd de a y b
 *          *alfa y *beta los enteros tales que m=alfa*a+beta*b
 */

void mcd_aux(long a, long b, long *alfa, long *beta, long *m)
{
	long s0,t0,s1,t1,s,t,q,r=1;
	
	if (a<b) return mcd_aux(b,a,beta,alfa,m);
	s0=1; t0=0;
	s1=0; t1=1;
	while (r != 0) {
		q = a/b;
		r = a%b;
		s = s0 - q*s1;
		t = t0 - q*t1;
		s0=s1; t0=t1;
		s1=s ; t1=t ;
		a =b ; b =r ;
	}
	*m=a;
	*alfa=s0;
	*beta=t0;
	return;
}