/****************************************************************************/
/* Orden de un entero modulo otro.                                          */
/*                                                                          */
/* Si a y m son enteros primos entre si, llamamos orden de a modulo m       */
/* al menor entero x tal que a^x = 1 (mod m)                                */
/*                                                                          */
/* Jaime Suarez <mcripto@bigfoot.com> 2003                                  */
/* en http://elparaiso.mat.uned.es                                            */
/****************************************************************************/

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


#include <limits.h>

#include <stdio.h>


typedef unsigned long ulong;

ulong mcd(ulong a, ulong b);
ulong exp_mod(ulong x, ulong y, ulong m);
ulong ord(ulong a, ulong m);

int main(void)
{
	ulong a,m;
	
	printf("Orden de a modulo m\n");
	printf("a= "); scanf("%ld",&a);
	printf("m= "); scanf("%ld",&m);

	if (mcd(a,m)!=1) {
		printf("a y m deben ser primos entre sí.\n");
		return 1;
	}

	printf("Ord(%d, %d)= %d\n", a, m, ord(a,m));
	return 0;
}

/* Maximo comun divisor de dos enteros positivos */

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

/* 
 * Orden de un entero modulo otro 
 * Observese que siempre tiene solucion pues al menos a^fi(m)=1 mod m
 * por el teorema de Euler
 */

ulong ord(ulong a, ulong m)
{
	ulong x;

	for (x=1; x<m; x++) 
		if (exp_mod(a,x,m) == 1) return x;
  
}


/* 
 * Calcula x^y mod n 
 * Devuelve 0 y un mensaje de error si se produce desbordamiento.
 */

unsigned long exp_mod(unsigned long x, unsigned long y, unsigned long n)
{
	unsigned long s,t,u;
	int i;

	s=1; t=x; u=y;
	while(u) {
		if (s>ULONG_MAX/t || t>ULONG_MAX/t) {
			fprintf(stderr,"Overflow en exp_mod.");
			return 0;
		}
		if (u&1) s=(s*t)%n;
		u>>=1;
		t=(t*t)%n;
	}
	return s;
}