/****************************************************************************/
/*  Calculo de raices primitivas                                            */
/*  Si r y n son enteros primos entre si con n>0 y ord(r,n)=fi(n)           */
/*  se dice que r es un raiz primitiva modulo n                             */
/*                                                                          */
/* 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);
ulong fi(ulong n);
void raices_prim(ulong n);

int main(int argc, char *argv[])
{
	ulong n;
	int i;
	char *bp;
	
	if (argc<2) {
		printf("Uso: %s <n1> <n2> ... \n",argv[0]);
		printf("Calcula las raices primitivas modulo n1, n2...\n");
		return 1;
	}

	for (i=1; i<argc; i++) {
		n = strtoul(argv[i], &bp, 0);
		raices_prim(n);
	}
}

/*
 * Escribir las ra¡ces primitivas de un entero n
 */

void raices_prim(ulong n)
{
	ulong phi_de_n, i;

	phi_de_n=fi(n);
	printf("Raices primitivas de %lu { ",n);
	for (i=1; i<n; i++) {
		if (mcd(i,n) != 1) continue;
		if (ord(i,n) == phi_de_n) printf("%lu ",i);
	}
	printf("}\n");
	return;
}

/* 
 * maximo comun divisor de dos enteros positivos 
 * sin utilizar divisiones, solo 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 mcd(a>>1,b>>1) <<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.
 *  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) Observese
 *  que siempre tiene solucion pues al menos a^fi(m)=1 mod m por
 *  el teorema de Euler. Ademas necesariamente x debe dividir a fi(m)
 */

ulong ord(ulong a, ulong m)
{
	ulong x, phi_de_m;
	
	phi_de_m= fi(m);
	for (x=1; x<=phi_de_m; x++) 
		if (phi_de_m%x == 0 && exp_mod(a,x,m) == 1) 
			return x;
}

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

ulong exp_mod(unsigned long x, unsigned long y, unsigned long n)
{
	ulong 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;
}

/*
 * Funcion fi de Euler, numero de enteros menores que uno dado
 * y primos con el. 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 (n%2 == 0) {
		f *= (1-1.0/2);
		while (n%2 == 0) n/=2;
	}
	
	p=3;
	while (p*p <= n) {
		if (n%p == 0) {
			f *= (1-1.0/p);
			while (n%p == 0) n/=p;
		}
		p += 2;
	}
	if (n!=1) f *= (1-1.0/n);

	return (ulong)f;
}