Introducción a la Factorización De Semiprimos (RSA)

[kub0x]

Normalmente, en este foro hablamos de la criptografía estándar, aquella que corre en nuestros dispositivos. Muchos sabreís que algunos de los esquemas criptográficos basados en clave pública o en una única clave (simétricos) sufrirán debilidades cuando el cuántico se torne realidad.

Esto es debido a que problemas como la factorización de enteros, el cálculo del logaritmo discreto, ya sea el aplicado a curvas elípticas o en finite fields podrán ser resueltos. Una respuesta obvia sería aumentar el tamaño de clave en la parametrización de estos esquemas, pero implicaría hacerlo cada vez que la computación cuántica avance.

Hoy vamos a explicar una forma de factorizar números compuestos por dos primos, es decir, semiprimos. Es un problema conocido pues nos daría la habilidad de romper esquemas como RSA. Nótese que RSA también presenta el RSA Problem y está basado en calcular raices e-ésimas mod n.

Antes de empezar, quiero decir que el lector tiene que tener un buen conocimiento en matemáticas: tanto en teoría de números como en algebra lineal.

Factorización De Fermat (Punto Medio)

Sea [latex]N=p\cdot q[/latex] tal que [latex]p,q \in \mathbb{P}[/latex] es un producto de dos primos, por lo tanto [latex]N[/latex] es semiprimo.

Queremos encontrar los factores [latex]p,q[/latex] que componen el módulo utilizado por una clave pública RSA. Esto nos permitiría generar firmas digitales válidas, además de intervenir en el TLS Handshake bien en falsificando valores cuando RSA se usa para firma, o bien, descifrando la PMK (Pre Master Key) lo cual nos permitiría hacer un Man In The Middle en plano, como antiguamente

Hace unos años, cuando empecé en esto, mi problema favorito era el de la factorización de enteros, curiosamente, la de semiprimos. Resulta que reinventé la rueda al descubrir que entre dos primos siempre existe un punto medio [latex]r[/latex] el cual equidista con misma distancia [latex]d[/latex] sobre los primos [latex]p,q[/latex]. Si traducimos a cristiano esto sería:

[latex](r %2B d) \cdot (r-d)=n[/latex] donde [latex]r=\frac{p %2B q}{2}[/latex] y [latex]d=\frac{p-q}{2}[/latex]. Si vemos más alla nos damos cuenta que [latex](x %2B y)\cdot (x-y)=x^2-y^2=n[/latex], lo que nos da la expresión [latex]r^2-d^2=p\cdot q=n[/latex]. ¿Bonito verdad? ¿Pero cómo aplicarlo?

Un ejemplo sería [latex]p=17,q=11,n=17\cdot 11=187[/latex] entonces [latex]r=\frac{17 %2B 11}{2}=14[/latex] y [latex]d=17-14=3[/latex]. Si nos imáginamos los enteros positivos alineados en una línea horizontal, vemos que [latex]14[/latex] equidista [latex]3[/latex] posiciones de [latex]17,11[/latex].

Para usarlo, al saber que [latex]r^2-d^2=n[/latex] hacemos [latex]\sqrt{n-r^2}=d^2[/latex] desde [latex]r=\sqrt n %2B 1[/latex] hasta dar con una solución entera. Pena que a Pierre de Fermat, se le ocurriera hace casi 400 años, conociéndose este método por Fermat Factorization.

En el mundo real nunca sabremos [latex]p,q[/latex] de antemano, sólo [latex]n[/latex] y nuestro interés se centra en encontrar los valores [latex]r,d[/latex] mediante una técnica llamada Quadratic Sieve o Criba Cuadrática.

Quadratic Sieve : Primera Fase

En teoría de números un número b-smooth es aquel que tiene factores primos menores o iguales que [latex]b[/latex], por ejemplo [latex]30=2\cdot 3 \cdot 5[/latex] es 5-smooth pero no 3-smooth. Establecemos [latex]b \in \mathbb{P}[/latex] como primo y seguimos nuestro análisis.

Si tomamos la congruencia [latex]x^2 \equiv y^2 \pmod{N}[/latex] entonces [latex]x^2 - y^2 \equiv 0 \pmod{N} \rightarrow x^2-y^2 \mid N[/latex] (divide a N) por lo tanto [latex](x %2B y)\cdot (x-y) \mid N[/latex] y tanto [latex]x %2B y[/latex] como [latex]x-y[/latex] son factores de [latex]N[/latex].

Empezamos tomando valores [latex]y[/latex] de la congruencia [latex]x^2 \equiv y \pmod{N}[/latex]. Para eso vamos generando valores [latex]x[/latex].  Expresamos [latex]y[/latex] en la base prima b-smooth [latex]y=p_1^{e_1} \ldots,p_b^{e_b}=\prod_{i=1}^{b}p_i^{e_i}[/latex].

Guardamos y reducimos mod 2 los exponentes de la base prima de cada residuo [latex]y[/latex] en una matriz de exponentes [latex]M_e \in F_2^{b \times k}[/latex], donde [latex]k[/latex] indica el número de congruencias [latex]y[/latex] que hemos obtenido de la criba. En cada columna guardamos los exponentes [latex]e_i[/latex] de cada congruencia [latex]x_i^2 \equiv y_i \pmod{N}[/latex]. Entonces vemos [latex]M_e[/latex] como:

[latex]M_e=\begin{bmatrix}e_{1,1}, \ldots, e_{k,1}\\ \vdots \ldots \vdots \\ e_{1,b} \ldots e_{k,b} \end{bmatrix}[/latex]

Segunda Fase

¿Por qué hemos reducidos los exponentes [latex]e_i[/latex] modulo 2 y encima en columnas? Tened en cuenta que si multiplicamos dos números, sus exponentes se suman. Si tras realizar la multiplicación todos los exponentes son pares/even, entonces [latex]a \cdot b[/latex] tiene raíz cuadrada entera. Si los exponentes de [latex]a,b[/latex] son [latex]\alpha=(\alpha_1,\ldots,\alpha_b), \beta=(\beta_1,\ldots,\beta_b)[/latex] entonces la multiplicación resulta en [latex]a \cdot b = \prod_{i=1}^{b}p_i^{\alpha_i %2B \beta_i}[/latex], si cada sumando [latex]\alpha_i %2B \beta_i[/latex] es par, la tupla modulo 2 resultado de [latex]\alpha %2B \beta[/latex] daría [latex]0[/latex] en cada posición.

De esta forma, encontramos residuos b-smooth cuyo producto resulta en una tupla de exponentes enteramente par, y como estamos en [latex]\mathbb{F}_2[/latex], obtenemos una tupla que consiste de 0's y de longitud [latex]b[/latex] por lo tanto  [latex]\alpha %2B \beta=(0,\ldots,0)_b[/latex].

Si revisaís los apuntes de Álgebra Lineal, el rango de una matriz nos indica la dimensión de la imagen. Una matriz se toma como una transformación lineal, cuenta con dominio y codominio. Como nuestra matriz tiene dimensión [latex]b \times k[/latex] entonces el sistema [latex]M_e \overrightarrow{x} = \overrightarrow{0}[/latex] nos daría en el vector [latex]\overrightarrow{x}=(x_1,\ldots,x_k)[/latex] la ansiada combinación lineal, puesto que una matriz por un vector, no es más que la suma de una combinación lineal de sus columnas. Las ecuaciones de tipo [latex]A \cdot \overrightarrow{x} = \overrightarrow{0}[/latex] se resuelven calculando el kernel o nullspace de la matriz [latex]A[/latex]. Si tenemos una matriz [latex]m \times n[/latex], con [latex]m>n[/latex] entonces [latex]rk(A) \leq n[/latex] por lo tanto [latex]Dim(A)=n[/latex] y [latex]Dim(A)-rk(A)=n-r=Dim(Ker(A))[/latex].

Por lo tanto la matriz que representa el kernel de [latex]M_e[/latex] tiene dimension [latex]k \times (k-rk(A)[/latex]. Si el rango de [latex]M_e[/latex] no baja de [latex]\min(b,k)[/latex] la única solución es [latex]\overrightarrow{x}=(0,\ldots,0)_k[/latex] llamada solución trivial. Para hallar soluciones no triviales el rango ha de ser inferior a [latex]\min(b,k)[/latex] y para esto, como las columnas de [latex]M_e[/latex] corresponden a los exponentes de cada residuo, si el rango baja sabemos que existe un producto de residuos que tiene raíz entera. De esta forma, encontramos una combinación lineal [latex]\sigma=(\sigma_1,\ldots,\sigma_k)[/latex] sobre las columnas que representan los exponentes, es decir, una suma finita de estas columnas (modulo 2) que da como solución una tupla de longitud [latex]b[/latex] con todos los elementos a [latex]0[/latex]. Entonces [latex]M_e \cdot \sigma = \overrightarrow{0}_b[/latex] tal y como deseabamos.

Tercera Fase

Ahora, vamos a encontrar los valores [latex]a,b[/latex] tal que [latex]a^2-b^2=(a %2B b)(a-b)=N[/latex]. Como hemos guardado todos los valores [latex]x_i,y_i[/latex] tal que [latex]x_i^2 \equiv y \pmod{N}[/latex], haciendo uso de la tupla [latex]\sigma[/latex] calculamos el producto de los valores [latex]x_i=\sigma_i[/latex] en [latex]a[/latex]. Lo mismo para los valores de [latex]y_i[/latex] y lo hacemos en [latex]b^2[/latex]. Entonces [latex]b=\sqrt{b^2}[/latex] y [latex]\gcd(a \pm b,N) \mid N[/latex].

La fase de criba finalizaría al dar con el kernel de [latex]M_e[/latex]. Como podeís ver es una fase muy amplía, donde necesitamos muchas congruencias para encontrar una combinación lineal de los exponentes que sea par. Además es necesario en todo momento almacenar las listas de las congruencias, tanto los [latex]x_i[/latex] como los [latex]y_i[/latex]. Todo ello mediante una representación matricial, donde cada columna corresponde a los exponentes de cada termino [latex]y_i[/latex] en la base b-smooth. Para concluir, vemos como el kernel nos permite extraer valores para los cuales existe una combinación lineas de las tuplas de exponentes que da cero. De esta forma sabemos que valores [latex]y_i[/latex] forman un cuadrado mediante su producto o multiplicación.

Ejemplo Quadratic Sieve

Vamos a poner un pequeño ejemplo para ilustrar la teoŕia anteriormente expuesta.

Tomamos [latex]p=4663, q=3581, N=16698203[/latex] la base está compuesta de los 50 primeros primos entonces es 229-smooth es decir [latex]B=\{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229\}[/latex]. Todos los valores[latex]y_i[/latex] que tratemos tienen que ser 229-smooth, es decir, factorizables entre productos de los primos en la base B.

Si empezamos con [latex]x=\lfloor \sqrt{N} \rfloor %2B 1[/latex], tras [latex]11[/latex] iteracciones, tenemos la lista [latex]X=(4107, 4114, 4125, 4371, 4401, 4657, 4745, 4836, 5217, 5221, 5355)[/latex] y la de residuos cuadraticos [latex]Y=(169246,226793,317422,2407438,2670598,4989446,5816822,6688693,10518886,10560638,11977822)[/latex].

La matriz de exponentes [latex]M_e[/latex] tiene dimensión [latex]50 \times 11[/latex] porque cada columna tiene longitud [latex]b=50[/latex] por lo tanto 50 filas y al haber 11 congruencias, de ahí sale. Nos enfrentamos al sistema [latex]M_e \cdot \overrightarrow{x} = \overrightarrow{0}_{50}[/latex]. Queremos una suma de columnas en [latex]M_e[/latex] que nos de completamente cero, y de esta forma sabremos que productos [latex]y_i y_j[/latex] dan un cuadrado mod N. La matriz [latex]M_e[/latex] se muestra a continuación:

\
La matriz [latex]M_e[/latex] tiene rango [latex]10[/latex] modulo 2. Por lo tanto, su kernel o nullspace tiene dimensión 1 y es representable mediante una matriz [latex]11 \times 1[/latex] en este caso [latex]nsM = \begin{bmatrix}0 \\ 0 \\ 1 \\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\1 \end{bmatrix}[/latex]. Entonces, la ecuación [latex]M_e \cdot x = 0 \rightarrow A \cdot (nsM \cdot c) = 0[/latex] tiene dos soluciones (una de ellas la trivial) porque la dimensión es 1, por lo tanto una tupla en [latex]\mathbb{F}_2[/latex] tiene [latex]2[/latex] valores, en este caso 0 ó 1. Por lo tanto fijamos [latex]\sigma=(0,0,1,0,0,0,0,0,0,0,1)[/latex]  y [latex]A \sigma = 0[/latex]. Comprobad vosotros mismos que la tercera columna y la onceava columna de la matriz [latex]M_e[/latex] suman 0 mod 2 en cada coordenada/posición.

Por lo tanto, sabemos que de la lista X multiplicamos los valores 3º y el 11º, esto es, [latex]a=4125 \cdot 5355=22089375[/latex] y lo mismo con la lista Y, nos queda [latex]b^2=317422 \cdot 11977822=3802024214884[/latex]. Calculamos [latex]\sqrt{b^2}=\sqrt{3802024214884}=1949878[/latex]. Entonces [latex]\gcd(a %2B b,N)=\gcd(22089375 %2B 1949878,N)=3581[/latex] y [latex]\gcd(a-b,N)=\gcd(22089375-1949878,N)=4633[/latex].

Aquí os dejo el código del algoritmo en Mathematica. Sé que lo puedo mejorar, lo he hecho en un par de horas para colgarlo junto a este tuto. En C++ utilizando la lib de GMP es igual de sencillo, mientras tengas las matemáticas de tu lado, todo va bien. El code es lenguaje de Mathematica aunque ponga el GeShI en MatLab, recuérdenlo.

Código (matlab) [Seleccionar] Expandir


bsmooth[x_] := (
   exps = Array[0 &, Length[plist]];
   found = False;
   i = 1;
   While[found == False && i <= Length[x],
    pos = Position[plist, x[[i]][[1]]] ;
    If[pos == {},
     found = True;
     ,
     pos = pos[[1]];
     exps[[pos]] = Mod[x[[i]][[2]], 2];
     ];
    i++;
    ];
   If[found == True,
    Return[{}];
    ,
    Return[exps];
    ];
   );

ComputeFactors[list_] := (
  k = 1;
  Do[If[list[[i]] != 0, k = k*list[[i]]], {i, 1, Length[list]}];
  Return[k];
  )

QSieve[n_, blen_, limit_] := (
  ctr = 1;
  xlist = Array[0 &, limit];
  ylist = Array[0 &, limit];
  plist = Table[Prime[i], {i, 1, blen}];
  A = ConstantArray[0, {blen, limit}];
  k = 1;
  While[ctr <= limit && k <= n,
   x = Floor[Sqrt[n] + k++];
   x2 = PowerMod[x, 2, n];
   bx2 = FactorInteger[x2];
   bx2 = bsmooth[bx2];
   If[bx2 != {},
    xlist[[ctr]] = x;
    ylist[[ctr]] = x2;
    A[[All, ctr]] = bx2;
    ctr++;
    ];
   ];
  If[MatrixRank[A, Modulus -> 2] < Min[limit, blen],
   nsM = NullSpace[A, Modulus -> 2] // Transpose;
   inNS = Mod[nsM.RandomInteger[{0, 1}, Dimensions[nsM][[2]]], 2];
   factorsX = inNS*xlist;
   factorsY = inNS*ylist;
   a = ComputeFactors[factorsX];
   b = Sqrt[ComputeFactors[factorsY]];
   Return[{GCD[a + b, n], GCD[a - b, n]}];
   ,
   Print[-1];
   ];
  )

Para llamar a la función, se necesita el semiprimo n, un número blen que indica la longitud de la base de primos b-smooth y un límite, cuantas congruencias ha de recolectar. Después de recolectar, intenta obtener una solución en la transpuesta del nullspace de la matriz de exponentes. El resto es componer los factores y hacer el GCD.

He capturado un .GIF para mostrar un ejemplo del algoritmo corriendo en Mathematica:

Vemos que el nullspace aquí ya no es [latex]11 \times 1[/latex] sino [latex]1000 \times 869[/latex]. Obviamente, el ejemplo que he puesto, era facilito pero el código trabaja para valores más altos.

Y nada, ha quedado demostrado que la criba cuadrática encuentra una solución mientras que seamos capaces de factorizar números en una base b-smooth. Además de necesitar una combinación lineal válida que nos asegure la obtención de un cuadrado mediante producto de residuos cuadráticos.

[KMALIGNO]

Muy buena la informacion, soy nuevo por aqui, gracias por el aporte, espero  poder tambien dar mi aporte