Hola.
Estoy teniendo muchos problemas con la dichosa rutinita esta y aún no funciona. Para saberlo le meto un número que sé que es primo (porque aparece en alguna lista de números primos en internet) y me dice que no es primo así que no está bien.
Os pongo lo que tengo hecho a ver si encontramos la causa:
Puede salir un tocho porque voy a poner muchas anotaciones para que no os perdais con GMP pero bueno ...
Aquí empiezo. Incluyo gmp para poder trabajar con números grandes, defino dos funciones y creo una variable 104729 que es un número primo conocido para pasar el test. El resultado debe ser que ese número es primo.
Miller Rabin tiene un factor de error de 1/4 así que siempre se debe repetir el test varias veces para asegurarse que siempre que ese número pasa el test se devuelve que es primo. Si una sola vez saliera que es compuesto ya sabríamos que el númeor es compuesto.
# include <gmp.h>
# define ES_PRIMO 1
# define NO_ES_PRIMO 0
int miller_rabin_test ( mpz_t n, int k );
int ComprobarDivisibilidad( mpz_t n );
int main(void)
{
int i;
char *c="104729";
mpz_t n;
mpz_init(n); //inicio la variable n
mpz_set_str(n,c,10); //meto el número a probar en n
if ( miller_rabin(n, 10) == ES_PRIMO ) { // 10 es el número de comprobaciones que haré
printf("%s es primol\n",c);
} else {
printf("%s no es primol\n",c);
}
}
Ahora el test de Miller Rabin en sí mismo. Para facilitar la rapidez lo primero que hago es ver si es divisible por los primeros diez mil números primos conocidos. Si es divisible por alguno es que el número es compuesto y no hay que perder mas tiempo.
int miller_rabin ( mpz_t Num, int Comprobaciones )
{
int r, s, EsPrimo=ES_PRIMO;
mpz_t NumDec, ValorAleat, ValorAux, d, ValorTemp;
if (ComprobarDivisibilidad(Num) == 0) { // comprobar si es divisible entre los primeros
EsPrimo=NO_ES_PRIMO; // 10.000 números primos
}
Omito la función ComprobarDivisibilidad para no hacer mas tocho y me centro en Miller Rabin.
Miller Rabin requiere que se encuentre un par de números R y S tal que satisfagan que (N-1)=R * 2^S. Además R debe ser impar.
Num es el número a comprobar.
NumDec es Num-1
mpz_init(d);
mpz_init(NumDec);
mpz_sub_ui(NumDec,Num,1); //NumDec=Num-1
mpz_set(d,NumDec); // d=Num-1
s = 0;
while ( mpz_even_p(d) ) { //mientras d sea par
mpz_fdiv_q_2exp(d, d, 1); // d=d/2
s++;
}
Hasta aquí todo bien. Si llego hasta aquí es que el número a comprobar no es divisible entre los primeros 10.000 primos y he obtenido un valor R y S que requiere MillerRabin
Inicializo valores aleatorios (la semilla)
gmp_randstate_t semilla;
gmp_randinit_default(semilla);
gmp_randseed_ui (semilla, time(NULL));
Inicializo algunos valores que necesito
mpz_init(ValorAleat);
mpz_init(ValorAux);
mpz_init(ValorTemp);
mpz_sub_ui(ValorTemp,Num,4); //ValorTemp=num-4
Aquí comienza el bucle de comprobación. Como dije MillerRabin tiene un factor de error y para asegurarnos que la respuesta es correcta en vez de hacer el test una vez lo hacemos varias. En mi caso 10 porque llamé a la función MillerRabin con el parámetro Comprobaciones=10.
Lo de EsPrimo es un valor boleano. Si antes resultó que es divisible entre algún número primo que no haga el bucle y se vaya al final.
while ( Comprobaciones-- > 0 && EsPrimo) {
//Obtener un valor aleatorio en ValAleatorio entre 2 y Num-2
mpz_urandomm(ValorAleat,semilla,ValorTemp); //Obtener un valor aleatorio entre 0 y num-4
mpz_add_ui(ValorAleat,ValorAleat,2); // ValorAleat=ValorAleat+2
// Calcular ValorAux=ValorAleat^d mod Num
mpz_powm(ValorAux,ValorAleat,d,Num);
if ( mpz_cmp_ui(ValorAux, 1) == 0 ) { //si valoraux=1 es probable primo
EsPrimo=ES_PRIMO; // es primo pero hay que comprobar mas
} else {
EsPrimo=NO_ES_PRIMO; //de momento no es primo.
}
if ( !EsPrimo && mpz_cmp(ValorAux, NumDec) == 0 ) { //si valoraux= Num-1 probable primo
EsPrimo=ES_PRIMO; //de momento es primo pero hay que comprobar mas
} else {
EsPrimo=NO_ES_PRIMO; //de momento no es primo.
}
if ( !EsPrimo ) { //si antes fue visto como no primo
for(r=0; r < s-1; r++) {
mpz_powm_ui (ValorAux, ValorAux, 2, Num); //ValorAux=ValorAux^2 mod Num
if ( mpz_cmp_ui(ValorAux,1) == 0) { // si ValorAux=1 es compuesto
EsPrimo=NO_ES_PRIMO;
break;
}
}
if ( mpz_cmp(ValorAux,NumDec) != 0 ) { //si ValorAux<>Num-1 no es primo
EsPrimo=NO_ES_PRIMO;
}
if ( !EsPrimo ) {
break;
}
}
} //fin del whileAhora al final solo queda devolver el resultado EsPrimo que será verdadero o falso y limpiar memoria
mpz_clear(ValorAux);
mpz_clear(ValorTemp);
mpz_clear(d);
mpz_clear(NumDec);
mpz_clear(ValorAleat);
return (EsPrimo);
}
Pues no me funciona y no veo la razón. Se que estoy muy cerca pero hay algo que se me escapa.
Ese código se puede optimizar mucho. Lo he puesto así porque se ven mas claros los pasos que doy para seguir a MillerRabin.
Mostrar Temas
