/****************************************************************************/
/* Resuelve una congruencia lineal aplicando el teorema de Euler.           */
/* Teorema de Euler: Si m es un entero positivo y a es un entero primo con  */
/* m entonces a^fi(m) = 1 (mod m) donde aparece la funcion fi de Euler,     */
/* número de enteros menores que uno dado y primos con el.                  */
/*                                                                          */
/* Por tanto si mcd(a,m)=1 entonces a.a^(fi(m)-1)= a^fi(m) = 1 (mod  m) y   */
/* dada la congruencia ax=b (mod m) tendremos que a^(fi(m)-1).ax = x =      */
/* a^(fi(m)-1).b (mod m) que es la solucion de la ecuación.                 */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://elparaiso.mat.uned.es                                            */
/****************************************************************************/

#include <stdio.h>

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

typedef unsigned long ulong; ulong fi(ulong n);
ulong fi(ulong);
long mcd(long, long);

int main()
{
	long a,b,c,i;
	ulong m;
	
	printf("Este programa resuelve congruencias del tipo\n");
	printf("ax = b (mod m) siendo a y m primos entre si.\n");
	printf("a= "); scanf("%ld",&a);
	printf("b= "); scanf("%ld",&b);
	printf("m= "); scanf("%ld",&m);

	if (mcd(a,(long)m)!=1) {
		printf("\na y m deben ser primos entre si.\n");
		return 1;
	}
	
	for (c=1,i=0; i<fi(m)-1; i++) c= (c*a)%m;
	printf("\n%ldx = %ld mod %ld ; ",a,b,m);
        printf("x= %ld (mod %ld)\n",(c*b)%m,m);

        return 0;
}
	
/*
 * Funcion 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;
	}
	
	for (p=3; p*p <= n; p+=2) {
		if ((float)n/p == n/p) {
			f *= (1-1.0/p);
			while ((float)n/p == n/p) n/=p;
		}
	}
	if (n!=1) f *= (1-1.0/n);

	return (ulong)f;
}

/*
 *  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);
}