Test de primalidad Miller-Rabin con GMP

[soplo]

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.

Código [Seleccionar] Expandir

# 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.

Código [Seleccionar] Expandir

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

Código [Seleccionar] Expandir

  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)

Código [Seleccionar] Expandir

gmp_randstate_t semilla;
gmp_randinit_default(semilla);
gmp_randseed_ui (semilla, time(NULL));

Inicializo algunos valores que necesito

Código [Seleccionar] Expandir

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.

Código [Seleccionar] Expandir

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 while
Ahora al final solo queda devolver el resultado EsPrimo que será verdadero o falso y limpiar memoria

 

Código [Seleccionar] Expandir

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.

[APOKLIPTICO]

Wii! Copado, genial, realmente muy bien, como siempre, un par de críticas constructivas, si me permitis...
1) Comprobar la divisibilidad con los primeros 10.000 primos, te va a tomar bastante tiempo, que realmente no es necesario si utilizas miller rabin. Las comprobaciones que deberías hacer para eliminar rápidamente números compuestos triviales son:
Si es par (num%2 == 0) y si es "2" o menor a "2" (1, 0 o negativos).
Si es par, la respuesta seria obviamente que no es primo, la comprobación de paridad se debe hacer despues de comprobar si es igual a "2" ya que "2" es par y primo (el unico caso). "0" y "1" son valores inválidos o no primos (ya que las condiciones necesarias para que un número sea primo es que sea divisible por 1 y por si mismo, en caso de "1" si mismo es igual a "1" asi que no tiene sentido y "0" no es divisible por si mismo). Si es "2" entonces es primo.

Probar con los primeros 10000 primos, va a disminuir mucho tu cálculo de primos por segundo.

2) Quizas convendría que el número aleatorio, no sea mayor a 1000 o a 500, para disminuir tiempo de cálculo. Total, solo necesitas a lo sumo 10 números.
3) No entendí muy bien algunas partes del código, las tendría que ver con mayor detenimiento, pero quiero ver si entendiste bien el algoritmo, porque me suena que hay algo que no está bien:

a) Si a^d mod n != 1 y != -1 (esto significa n -1). Entonces n es probablemente primo.
b) Hay que probar todos los a^(d^(2*s)) mod n. Siendo "s" la variable que va desde 1 hasta las veces que lo dividiste. Vas incrementando "s" hasta que de -1 o hasta que llegues a las veces que lo dividiste. Si te dio -1 alguna de las potencias, es probablemente primo, sino, es compuesto.

PD: La próxima utilizá los tags [code=cpp

[soplo]

a mi también me parece que al final he acabado mareando la perdiz y poniendo algo mal. Por eso no me da bien. Tengo en mente hacer un manual de operaciones GPM porque es bastante antipático tal como viene en el manual de gnu mp 5.1.

En cualquier caso si tienes alguna duda con alguna línea no dudes en decirlo a ver si encontramos donde está el error (que seguro que es de bulto porque ya no se cuantas veces lo he tocado)

je je je

[APOKLIPTICO]

Estoy luchando con gmp para instalarlo bajo mingw, pero me sigue tirando undefined references...
Ahora cuando lo haga andar, me parece que ya se donde está el problema, en un lugar pones que si valoraux == 0 es compuesto, en realidad, si valoraux = -1 entonces es primo... Se tienen que probar todos los "s" para que sea compuesto.

[WestOn]

Muy bueno, al final te salió bien me alegro

Saludos

[APOKLIPTICO]

Estuve haciendo pruebas, lo depuré a mas no poder y encontré un error en la función mpz_powm, que en vez de dar "n-1" da "1"...
Ese error ocurre cuando utilizo numeros un poco grandes. El problema debe de estar en la función, aunque me parece extraño.

[soplo]

Para facilitar la comprensión he modificado el código anterior (solo para este hilo) a fin de ver si encontramos donde está el error. He omitido las funciones GMP y las he escrito en C clásico con el único interes de ver si hay algún error lógico.

A mi me parece que hay algo del algoritmo que no he puesto bien y por eso no me funciona adecuadamente, pero es que no lo veo.

Tal como está no funciona porque solo he cambiado textos, De hecho sé que faltan puntos y comas y tampoco están ni stdlib.h ni math.h pero da igual. Solo se trata de ver donde está el error para que no calcule bien.

Ahi va

Código [Seleccionar] Expandir

int miller_rabin ( long  Num, int Comprobaciones )
{
  int r, s, EsPrimo=ES_PRIMO;
  long 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

Código [Seleccionar] Expandir

  NumDec=Num-1; //NumDec=Num-1
  d=NumDec

  s = 0;
  while ( mpz_even_p(d) ) {  //mientras d sea par
    mpz_fdiv_q_2exp(d, d, 1); // d=d/2
    s++;
  }

// Obtener la semilla para números aleatorios

Código [Seleccionar] Expandir

srand(time(NULL))
ValorTemp=Num-4 //este valor se utiliza para el límite superior de número aleatorio a obtener

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.

Código [Seleccionar] Expandir

while ( Comprobaciones-- > 0 && EsPrimo) {
    //Obtener un valor aleatorio en ValAleatorio entre 2 y Num-2
   
    ValorAleat=rand() * (ValorTemp)
    ValorAleat+=2 //obtener un valor aleatorio entre 2 y n-2

// Calcular ValorAux=ValorAleat^d mod Num
    ValorAux=pow(ValorAux,d) % Num

    if ( ValorAux == 1 ) { //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 && ValorAux == NumDec ) { //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++) {     
          ValorAux=pow(ValorAux,2) % Num
          if ( ValorAux == 1 ) {  // si ValorAux=1 es compuesto
             EsPrimo=NO_ES_PRIMO;
             break;
          }
       }
       if ( ValorAux<>NumDec ) { //si ValorAux<>Num-1 no es primo
             EsPrimo=NO_ES_PRIMO;
       }
       
       if ( !EsPrimo ) {
          break;
       }
    }   
  } //fin del while

¿donde demonios está el error?

[APOKLIPTICO]

Yo lo compile y todo, el problema esta en GMP, por alguna razon la funcion del modulo exponencial, devuelve 1 cuando deberia devolver n-1.