El Descubrimiento de los Números Primos: Usar un ODROID-C2 para hacer un historial matemático

“Se sabe que el problema de distinguir números primos de números compuestos y de resolver estos últimos en sus factores primarios es conocido como uno de los más importantes y útiles en aritmética.” − Carl Friedrich Gauss

En este artículo, proporcionare una introducción sobre algunos de los aspectos algorítmicos de las pruebas de primalidad, los ejemplificaré utilizando la utilidad de Linux bc y describiré algunos de los algoritmos más avanzados utilizados en la famosa prueba de primalidad Lucas-Lehmer (LL) para los números de Mersenne y la propia implementación del autor en su programa Mlucas. Este software tiene una versión optimizada para la funcionalidad de hardware vector-aritmética disponible en la familia de procesadores ARMv8, específicamente en el SBC ODROID-C2. Sin embargo, ten en cuenta que el software también se puede compilar en SBC ODROID que no son v8, aunque no utilizan las instrucciones vectoriales. Puesto que la página LEAME de Mlucas (vinculada más abajo) proporciona instrucciones detalladas de compilación para una amplia variedad de plataformas de hardware, incluidas las SBC ODROID, me voy a centrar en las matemáticas y los algoritmos, a un nivel que debería ser comprensible para cualquier persona con conocimientos básicos de álgebra y programación de computadoras.

Raíces primitivas y primalidad

LL es un ejemplo de lo que se conoce como prueba de primalidad no factorial, que hace referencia al hecho de que no se necesita conocimiento alguno sobre la factorización en primos del número de entrada N, o módulo, aunque normalmente realizamos una "factorización de prueba" pre-tamizada a fin de verificar los números para los pequeños factores primos antes de recurrir a la prueba LL. Tales pruebas se basan en propiedades algebraicas muy intensas de campos numéricos que están fuera del alcance de este artículo, pero en esencia equivalen a probar si existe una raíz primitiva  del grupo matemático definido por el módulo de multiplicación N, lo que significa una raíz del máximo orden posible N−1. (Esto suena a matemáticas complejas, pero en realidad se reduce a términos muy simples, como ya veremos más adelante). Tal raíz existe si y solo si (es decir, lo contrario) N es primo. Por ejemplo, si calculamos las potencias sucesivas de 2 módulo N = 7, obtenemos la secuencia 2,4,1, ... que se repite con la longitud 3, lo que significa que 7 es compuesto o 2 no es una raíz primitiva (mod 7). Si en cambio probamos potencias de 3 obtenemos la secuencia 3,2,6,4,5,1 que es de la máxima longitud posible N−1 = 6 para el grupo multiplicativo (mod 7), por lo tanto, concluimos que 7 es primo a modo de haber encontrado que 3 es una raíz primitiva. Si en cambio probamos las secuencias de potencia resultantes de los mismos módulos bases del compuesto 15 obtenemos 2^k (mod 15) = 2,4,8,1, ... y 3^k (mod 15) = 3,9, -3, -9, ..., para el índice k = 1,2,3, .... Observamos que ambas secuencias se repiten con periodicidad 4, que en primer lugar es menor que N−1 y en segundo lugar no divide N−1 (es decir, no hay posibilidad de que una repita en la secuencia de 2^k que acaba en la ranura N-1), y que la secuencia 3^k ni siquiera contenga un 1, mientras que ambas secuencias (mod 7) contienen uno o más unos, en particular ambos tienen 1 en el intervalo (N−1) Es decir, para N = 7 tenemos a^N-1 ≡ 1 (mod N) para todas las bases, no un múltiplo de 7, donde la triple-barra-igual es el símbolo para la equivalencia modular, es decir, los 2 lados de la relación son iguales cuando se reduce el módulo N. Para las secuencias (mod 15), por otro lado, la no ocurrencia de 1 como potencia N−1 en cualquiera de las dos es suficiente para demostrar el compuesto 15.

La computación de estas secuencias de mod-potencia para módulos grandes no es práctica para grandes N, de modo que la idea de calcular solo la potencia (N-1) de la base es crucial, ya que eso requiere el cálculo de unos meros intermedios O( lg N), donde lg es el logaritmo binario. La prueba resultante resulta ser rigurosa solo en el sentido de identificar correctamente los números compuestos, siendo meramente probabilística para los primos porque arroja un resultado de falso positivo '1' para un pequeño porcentaje de módulos compuestos además de los primos, pero puede ser modificada para obtener una prueba determinante (rigurosa) para una fracción convenientemente grande de estos números primos falsos problemáticos.

El "pequeño" teorema de Fermat y la prueba de probable-primalidad

Tanto las rigurosas pruebas de primalidad como las llamadas pruebas de probable-primalidad se basan de una manera u otra en la existencia de una raíz primitiva, y es muy útil utilizar la variedad probable-primo para ejemplificar el tipo de aritmética necesaria para implementar esta prueba, especialmente para módulos grandes. El "padre de la teoría de números", Pierre de Fermat, a principios del siglo XVII ya había observado y formulado en un teorema lo que hemos señalado anteriormente, que para cualquier primo p, si a es coprimo para (no tiene factores en común con) p - puesto que p primo significa que a no debe ser un múltiplo de p – entonces

a^(p−1) ≡ 1 (mod p).  [*]

Esto se conoce actualmente7 como el "pequeño teorema" de Fermat para diferenciarlo de su "último teorema" más famoso (pero de menor importancia práctica), acerca de las soluciones sobre los enteros de la ecuación aⁿ+bⁿ = cⁿ; la prueba resultante aplicada a los números de caracteres desconocidos se denomina, indistintamente, como una prueba de compostura probada de Fermat o de probabilidad probable de Fermat. La primera denominación se refiere al hecho de que un entero N satisfactorio a^N−1 ≡ 1 (mod N) para una base coprimo a N es muy probable que sea primo, para una gran muestra estadísticamente hablando, el segundo al hecho de que los números que fallan con este criterio son ciertamente compuestos, incluso si no se ha encontrado un factor explícito de N.

Tenga en cuenta que el inverso de [*] no es válido, es decir, hay pares enteros coprimos a,N para los cuales a^N−1 ≡ 1 (mod N) pero donde N no es un primo; por ejemplo, usando la utilidad de linux 'bc' se puede ver que el compuesto N = 341 = 11 × 31 satisface el teorema para la base a = 2, invocando bc en un intérprete de comandos y simplemente escribiendo '2^340%341'. Esto subraya la importancia de probar bases múltiples si el primero genera  a^N−1 ≡ 1 (mod N): para N primo, cada entero en 2,...,N-2† es coprimo para N, y por lo tanto todas estas bases producen 1 en el resultado, mientras que para el compuesto N, casi siempre basta con probar un pequeño número de bases para revelar N como compuesto. Decimos "casi siempre" porque existe una clase especial de números enteros conocidos como números de Carmichael, que pasan la prueba de Fermat para todas las bases que son coprimo a N; el más pequeño es  561 = 3×11×17. Los números de Carmichael solo revelan su carácter compuesto si elegimos una base a para la prueba de Fermat que no es coprimo para N, i.e. es decir, si encontramos un factor de N. En la práctica siempre se debe verificar primero N para pequeños factores primos hasta algún limite (fijado por el coste de dicha factorización de prueba) antes de someter a N a una prueba de primalidad probable al estilo de Fermat.

† Saltamos tanto 1 como N-1 como bases potenciales porque, siendo g ≡±1 (mod N),estos dos valores "bookend" de ambos producen trivialmente a^N−1 ≡ 1 para todos los N impares.

Para la base a = 2 hay precisamente 10403 de estos pseudoprimes de Fermat, un número minúsculo comparado con los cerca de 200 millones de números primos por debajo de ese límite, de modo que incluso usar simplemente esta base única nos proporciona una manera notablemente efectiva de determinar si un número es susceptible de ser primo. Por ejemplo, se puede combinar una prueba de pseudoprimo Fermat base-2 con una tabla preconfigurada de los compuestos conocidos inferior a 2³² que pase la prueba para producir un algoritmo de primalidad determinista muy eficiente para números por debajo de ese límite. Sin embargo, en un contexto más general de probar números de tamaño arbitrario, es importante darse cuenta que hay ciertas clases de números, los cuales son pseudoprimos Fermat base-2, independientemente de si son primos o compuestos. Las dos clases más conocidas son, en primer lugar, los números de Mersenne M(p) = 2^p−1 (para los cuales restringimos la definición a los exponentes principales ya que se requiere que un número de esta forma tenga la oportunidad de ser primo); por ejemplo,, 2¹¹−1 pasa la prueba aunque tiene un factor de 23 × 89.

La segunda clase de estos números son los números de Fermat F_m = 2^2m+1. Parece que los primeros cinco números de ese tipo, F₀ - F₄ = 3,5,17,257,65537 son lo suficientemente pequeños como para ser susceptibles a la división de prueba de lápiz y papel, y se pone de manifiesto fácilmente que son primos de esta manera, junto con el hecho que todos los que hayan pasado la prueba nombrada puede haber llevado a Fermat a hacer su famosa conjetura incorrecta de que todos estos números son primos. Euler refuta esta conjetura décadas más tarde, al mostrar que 641 es un factor primo de F5, y el trabajo posterior ha llevado a la creencia general de que con toda probabilidad los cinco números más pequeños son de hecho los únicos primos en la secuencia. Simplemente cambiando la base de la prueba de Fermat a, digamos, 3, es suficiente para distinguir los primos de los compuestos en esta secuencia de números, pero parece que esta idea nunca se le ocurrió a Fermat.

Exponenciación modular eficiente

Para llevar a cabo la exponenciación modular, usamos una técnica general, o mejor, un conjunto relacionado de técnicas, conocido como jerarquía de alimentación binaria. El nombre se refiere al hecho de que los diversos enfoques se basan en analizar los bits del exponente representado en forma binaria. Para animar al lector a calcular junto con la lectura, utilizamos la calculadora de precisión arbitraria integrada POSIX, bc, que es relativamente lenta en comparación con los programas gratuitos de teoría de números de gama superior como Pari/GP y la librería multiprecisión Gnu, GMP, pero nos puede ser extremadamente útil para este tipo de 'prototipo rápido' algorítmico básico. Invocamos la calculadora en modo de número entero por defecto simplemente a través de 'bc' en un terminal; 'bc -l' invoca la calculadora en modo de punto flotante, en la cual la precisión puede ajustarse para adaptarse usando el valor del parámetro 'scale' (cuyo valor predeterminado es 20 dígitos decimales), y '-l' define la librería matemática estándar, que contiene un puñado de funciones útiles que incluyen el logaritmo y el exponente natural, el seno trigonométrico, el coseno y la arcotangente (a partir de los cuales se pueden desarrollar otras cosas, por ejemplo '4*a(1)' calcula π con la precisión fijada por el valor actual de escala usando ese arctan (1) = π/4), y la función de Bessel del primer tipo. La base aritmética de las entradas y salidas de bc se puede controlar modificando los valores de ibase y obase de sus valores por defecto del 10, por ejemplo, para ver 23 en binario, escriba 'obase = 2; 23 'que fielmente produce 10111; Resetea la escala a su valor decimal por defecto mediante 'scale = 10'.

Obsérva que para los módulos en general, y en particular los muy grandes, no podemos simplemente calcular la potencia en el lado izquierdo de [*] y reducir el módulo resultante N, ya que los números son lo suficientemente grandes como para saturar incluso nuestro extenso almacenamiento informático. En general, las potencias se duplican en longitud para cada bit del exponente, por lo que elevar la base 2 a un exponente de 64 bits da un resultado del orden de 264 o 18,446,744,073,709,551,616 bits, o más de 2000 petabytes, ejemplificado muy bien el famoso problema del tablero de ajedrez y el trigo en las matemáticas de las series geométricas. Para probar un número del tamaño del primo Mersenne, más recientemente descubierto, tendríamos que hacer decenas de millones de este tipo de reiteraciones de duplicado de tamaño, entonces, ¿cómo es posible qué esto esté disponible en el hardware de cálculo? Volvemos a las propiedades de la aritmética modular, una de las más importantes es que al calcular un gran 'powermod' de este tipo, podemos hacer una reducción modular en cada paso del camino, siempre que sea conveniente hacerlo. Así, en la práctica, se usa una jerarquía binaria de potencia para dividir la exponenciación en una serie de escuadrados y multiplicaciones, a cada paso le sigue una reducción modular del intermedio resultante.

Compararemos y contrastaremos dos planteamientos básicos para el cálculo de a^b (mod n), de potencia binaria, que recorre los bits del exponente b en direcciones opuestas. Ambos hacen una cuadratura por bit de b, así como algunas multiplicaciones más generales cuyo conteo preciso depende del patrón de bits de b, pero nunca puede exceder el de los escuadrados. El método de derecha a izquierda activa un acumulador  y = 1 y un cuadrado actual z = a, luego para cada bit en n empezando con el bit más a la derecha (unos), si el bit actual = 1, multiplicamos el acumulador por el cuadrado actual z, luego de nuevo el cuadrado z para prepararse para el próximo bit a la izquierda. Aquí hay una simple función bc definida por el usuario que ilustra esto: por razones de simplicidad, hemos omitido algunos pre-procesamientos básicos que se incluirían para verificar la cordura de las entradas tales como las pruebas de módulo cero y de exponente no negativo:

/*
 * right-to-left-to-right binary modpow, a^b (mod n):
 */

define modpow_rl(a,b,n) {
  
  auto y,z;
  y = 1; 
  z = a%n;

  while (b) {

     if(b%2) y = (y*z)%n;

     z = (z*z)%n;

     b /= 2;

  }

  return (y);

}

Lo dejamos como un ejercicio para que el lector implemente una simple optimización que agregue una anticipación dentro del ciclo de forma que la actualización del cuadrado actual solo se realice si hay un bit próximo a la izquierda para procesar, lo cual es útil si la base a es grande pero el exponente es pequeño. Tras pegar el código anterior en tu intérprete de comandos bc, intenta hacer una prueba de pseudoprimo de Fermat de base 2 del primo Mersenne conocido M(2203): 'n=2^2203-1; modpow_rl(2,n-1,n)' (esto llevará unos segundos). Ahora vuelve a intentarlo con un exponente compuesto de tamaño similar, digamos 2205, y observa el resultado no unitario que indica que n = 2²²⁰⁵−1 también es compuesto. Como 2205 tiene pequeños factores primos 3,5 y 7, podemos utilizar aún más la función de módulo bc '%' para verificar que este número sea exactamente divisible entre 7,31 y 127.

Por supuesto, ya tenemos presente que una prueba de pseudoprimo de Fermat de base 2 devuelve 'primo probable' para todos los números de Mersenne, es decir, para todos los  2^p−1 con p primo, podemos usar la función modpow anterior para verificar esto, ahora modificando el valor de retorno a un 0 binario (compuesto) o 1 (pseudoprimo a la base dada):  'n=2^2207-1; modpow_rl(2,n-1,n) == 1' devuelve 1 a pesar de que el número de Mersenne en cuestión tiene un factor pequeño, 123593 = 56 × 2207 + 1. Repitiendo la misma prueba de pseudoprimo de Fermat pero ahora con la base 3 correctamente, revela que este número es compuesto. Observa el salto asociado en el tiempo de ejecución de bc. Esto parece reflejar algunas optimizaciones de casos especiales en la lógica interna de bc relacionadas con el reconocimiento de argumentos que son potencias de 2. Vemos una discrepancia de velocidad similar al repetir la prueba de pseudoprimo usando las bases 4 y 5, por ejemplo.

Nuestro próximo algoritmo de potencia procesa los bits en la dirección opuesta, de izquierda a derecha. Este método activa el acumulador y = a, que corresponde al bit fijado más a la izquierda, luego para cada bit hacia la derecha cuadramos y, y si el bit actual = 1, multiplicamos el acumulador por la base de potencia a. En un lenguaje de codificación como C podríamos, ya sea a través del compilador o a través de una pequeña macro de código de ensamblaje, implementar las funciones de bits mediante el metodo binario de dividir y vencer o accediendo de manera eficiente a cualquier instrucción de hardware disponible para liderar el recuento de ceros, y nuestra función de reversión de bits reverse () podría implementarse de manera eficiente utilizando una pequeña tabla de 256 bytes precomputados con bit invertidos, un ciclo para hacer intercambios de byte a los extremos izquierdo y derecho de nuestro exponente, y un paso final para desplazar hacia la derecha el resultado de 0-7 bits, dependiendo de dónde se encuentre el bit establecido más a la izquierda en el byte configurado más a la izquierda. En BC no tenemos esa funcionalidad bit a bit y debemos desplegar nuestras propias tareas de emulación ineficientes, pero como nuestro foco está en modpow de grandes operandos, el coste de tiempo de tales operaciones bit a bit es mínimo en comparación con las multiplicaciones modulares:

define bits(n) {
  auto ssave, r;
  ssave = scale;  scale = 0;  /* In case we're in floating-point mode */
  r = length(n)*3321928095/1000000000;
  while ( 2^r > n ) { r -= 1; }
  scale = ssave;
  return(r+1);
}

define reverse(n,nbits) {
  auto tmp;
  tmp = 0;
  while(nbits) {
    tmp = 2*tmp + (n % 2);
    n /= 2;
    nbits -= 1;
  }
  return(tmp);
}

/* left-to-right binary modpow, a^b (mod n): */
define modpow_lr(a,b,n) {
  auto y,len;
  len = bits(b);  b = reverse(b,len);
  y = a%n;  b /= 2;
  while(--len) {
    y = (y*y)%n;
    if(b%2) y = (a*y)%n;
    b /= 2;
  }
  return(y);
}

La necesidad de inversión de bits también se puede evitar implementando el algoritmo de forma recursiva, pero como bc es, digamos, menos que espectacular cuando se trata de soporte eficiente para la recursión, preferimos un algoritmo no recursivo en el caso de izquierda a derecha. Instamos a los lectores a pegar lo anterior en su intérprete de comandos bc, utilízalo para realizar nuevamente las pruebas de Fermat-pseudoprimo base-2 y base-3 en 2²²⁰⁷−1, y compáralas con las del algoritmo de derecha a izquierda. En mi interprete de comandos bc, el método de izquierda a derecha se ejecuta en casi la mitad del tiempo en el número compuesto de Mersenne antes mencionado, a pesar de que las entrañas de los bucles en nuestras funciones de alimentación RL y LR son bastante similares, cada una con un mod-cuadrado y un mod-multiplicador.

El motivo de la aceleración en el método LR se vuelve claro cuando examinamos con precisión qué operandos están involucrados en las dos operaciones respectivas de multiplicación de mod. En la potencia RL, multiplicamos el acumulador de corriente y y el mod-square actual z, que son del tamaño del módulo n. En la alimentación LR, multiplicamos el acumulador de potencia actual y la base a, siendo este último generalmente el orden de la unidad, u O (1) en la notación de orden asintótica estándar.

De hecho, para bases pequeñas, podemos reemplazar el producto a×y por una serie de cambiar y acumular bits hacia la izquierda de y si resulta ventajoso, aunque la línea inferior - con vistas a la implementación de hardware subyacente de estas operaciones de precisión arbitraria a través de la aritmética de múltiples palabras – que es la del algoritmo LR multiplicamos el vector y por la base escalar a, que es lineal en la longitud del vector en términos de coste. Para los exponentes generales cuyos bits se dividen por igual entre 0 y 1, solo nos damos cuenta de esta reducción de costes para los bits 1, pero nuestro ejemplo numérico particular implica un número de Mersenne, cuyos bits son 1, así el exponente de potencia binaria n- 1 tiene solo un bit 0 en la posición más baja, y si vector-times-vector mod-square y mod-multiply tienen aproximadamente el mismo coste (como en el caso de bc), reemplazamos este último por un vector-times-scalar reduciendo el tiempo de ejecución aproximadamente a la mitad, como se ha observado.

Prueba determinista de la primalidad

Si bien la prueba de Fermat-pseudoprimo es una forma eficiente de identificar si un número dado es probablemente primo, nuestro interés real es establecer rigurosamente el carácter, primo o compuesto del mismo. Por lo tanto, es importante complementar estas pruebas de primalidad probable con alternativas deterministas siempre que sea posible. Si dicha alternativa puede realizarse con un coste computacional similar, sería lo ideal, ya que la calculo a^N−1 (mod N) utilizando el enfoque modular de alimentación de binario de izquierda a derecha requiere del cálculo de O(lg N) escuadrados modulares de intermedios de tamaño N y un número similar de multiplicaciones modulares de tales intermedios por la base a, que generalmente es mucho menor que N, lo que significa que estas multiplicaciones escalares suplementarias no tienen significado en términos de coste algorítmico global asintótico de N grande. Hay razones para creer, aunque esto no se ha demostrado, que el coste de una cuadratura modular por bit de N es de hecho el óptimo que se puede lograr para una prueba de primalidad no factorial.

Por desgracia, parece que solo hay una clase limitada de números de formas muy especiales para las cuales existe una prueba determinista de la primalidad con un coste computacional similar; las más famosas son las dos clases mencionadas anteriormente. Para los números de Fermat utilizamos una generalización de la prueba de pseudoprimos de Fermat gracias a Euler, en la que uno calcula un a^(N−1)/2(mod N) y compara el resultado con ± 1, con el signo apropiado dependiendo de una particular propiedad algebraica de N. Para los números de Fermat, basta con tomar a = 3 y verificar si  3^(F_m−1)/2 = 3^22m−1 ≡ −1 (mod F_m), que requiere precisamente 2^m−1 mod-squarings de la simiente inicial 3. Lo suficiencite para determinar la primalidad de los números de Fermat se conoce como el teorema de Pépin. El lector puede usar cualquiera de las funciones modpow anteriores de bc para realizar la prueba de primalidad Pépin en un número de Fermat de hasta varios kilobits de tamaño; por ejemplo, para probar el F11 de 2 kilobits, 'n = 2^(2^11)+1; modpow_lr(3,(n-1)/2,n) == n-1'. Ten en cuenta que los algoritmos LR y RL se ejecutan aproximadamente al mismo tiempo en los números de Fermat, ya que la potencia calculada en la exponenciación modular de la prueba Pépin es una potencia de 2, por lo tanto solo tiene 1 bit en la posición más a la izquierda.

Para los números de Mersenne M(p) = 2^p−1, la prueba de primalidad se basa en las propiedades algebraicas de las llamadas secuencias de Lucas según el matemático francés Édouard Lucas, pero luego fue refinada por el teórico de los números Derrick Lehmer en una forma algorítmica simple conocida como la prueba de primalidad de Lucas-Lehmer: comenzando con cualquiera de los valores iniciales algebraicamente permitidos x (el más utilizado de los cuales es 4), realizamos precisamente actualizacones repetitivas p-2 de la forma x = x²−2 (mod M(p)); el número de Mersenne en cuestión es primordial si y solo si el resultado es 0. Por ejemplo, para p = 5, tenemos iteraciones no modificadas 4,14,194,37634; en forma modificada estas son 4,14,8,0, lo que indica que la iteración final 37634 es divisible por el módulo 31 y por lo tanto que este módulo es primo. Al igual que con nuestras pruebas de primalidad probable y la prueba de Pépin, tenemos una de esas iteraciones por cada bit del módulo, dar o recibir uno o dos..

Para dar una idea de las eficiencias relativas de estas pruebas de módulo especializado en comparación con las pruebas de primalidad determinística para los módulos generales, los más rápidos conocidos se basan en la aritmética de las curvas elípticas y se han utilizado para probar la primalidad de números que tienen alrededor de 20,000 dígitos decimales, mientras que las pruebas de Lucas-Lehmer y Pépin se han realizado hasta la fecha en números que se acercan a los 200 millones de dígitos. Puesto que, como mostraremos a continuación, el coste total de las pruebas especializadas es ligeramente mayor que la cuadrática en la longitud de la entrada, esta disparidad de tamaño del factor de 10 se traduce en una diferencia efectiva proporcionalmente mayor en la eficacia de la prueba. En términos de las pruebas de módulo especializado, la diferencia de velocidad es más o menos equivalente a una disparidad de cien millones de veces en la eficiencia de las pruebas basada en los tamaños de los que mantiene el record actual de registros para ambos tipos de pruebas.

Rápida multiplicación modular

A pesar del formato superficialmente distinto de las dos pruebas de primalidad anteriores, observamos que, tanto para la operación crucial de control de rendimiento, como en el caso de la prueba de pseudoprimo de Fermat, es la multiplicación modular quien toma el producto de un par de números de entrada y reduce el módulo de resultado en un tercio. (La operación de resta 2 adicionales de la prueba de LL es insignificante con respecto a estas escalas de trabajo asintóticas de grandes operandos). Para entradas de tamaño modesto podemos usar un algoritmo de multiplicación estándar de la "escuela de gramática" dígito por dígito, seguido de una división con resta por el módulo, pero de nuevo, para números grandes debemos ser bastante más inteligentes que esto.

La idea principal que hay detrás de los modernos algoritmos de multiplicación de grandes enteros se debe a Schönhage and Strassen (ver además http://numbers.computation.free.fr/Constants/Algorithms/fft.html para una exposición más matemática de la técnica), quienes reconocieron que la multiplicación de dos enteros equivale a la complejidad digial de las entradas. Esta idea permite que cualquiera de los conocidos algoritmos complejos de alta eficiencia del ámbito del procesamiento digital de señales se aproveche del problema. Como reflejo de la evolución del microprocesador moderno, las implementaciones más conocidas de dichos algoritmos utilizan la Transformación Rápida de Fourier (FFT) y hacen uso del hardware punto flotante del procesador, a pesar de los errores de redondeo asociados que hacen que las salidas complejas se desvíen de los valores de enteros puros, se tendrían que usar la aritmética exacta.

A pesar de esta inexactitud, el software FFT de alta calidad es notablemente agresivo en la cantidad de errores de redondeo que uno puede soportar sin corromper fatalmente las largas cadenas de multiplicación involucradas en pruebas de primalidad de enteros grandes: por ejemplo, la prueba de Lucas-Lehmer que descubrió el primo Mersenne en record de tamaño más reciente, 2^77232917−1, usaba una FFT de doble precisión que dividía cada iteración de 77232917 bits en 2²² palabras de entrada de 18 o 19 bits de longitud (hablando con exactitud 1735445 = 77232917 (mod 2²²) de palabras grandes y las 2458859 restantes de tamaño más pequeño), así se usa algo más de 18 bits por 'dígito' de entrada de la complejidad discreta. Esto proporciona salidas de complejidad de coma flotante que tiene partes fraccionales (es decir, errores de redondeo acumulados) tan grandes como casi 0,4, lo cual está notablemente cercano al fatal nivel de error 0.5 "No sé si redondear hacia arriba o hacia abajo este resultado de complejidad de punto flotante inexacta".

No obstante, las múltiples ejecuciones de verificación utilizando implementaciones FFT desarrolladas independientemente en diferentes (y mayores, por lo tanto, tienen errores de redondeo mucho más pequeños) longitudes de transformación, en varios tipos diferentes de hardware, confirmaron la exactitud de la prueba de primalidad inicial. Para un módulo de n bits, el coste de cálculo n-asintótico de tal multiplicación de FFT es O(n lg n), aunque te das cuenta de esto en la práctica, especialmente cuando los operandos son lo suficientemente grandes como para exceder los tamaños de los cachés de datos de varios niveles utilizados en los microprocesadores modernos son un desafío algorítmico y de movimiento de datos muy poco significativo. Dado que nuestras diversas pruebas de primalidad requieren O(n) de estos multiples pasos, la estimación de trabajo para una prueba de primalidad basada en FFT es O(n² lg n), que es solo un factor lg n mayor que el coste cuadrático de una sola escuela primaria multiplicar. El resultado es que para escribir un programa de prueba de primalidad de clase mundial, uno debe escribir (o hacer uso de) una complejidad basada en transformaciones de clase mundial.

Aproximadamente hace 20 años, estaba en la facultad de ingeniería en la Universidad de Case Western Reserve en Cleveland, Ohio, y estaba buscando una forma interesante de motivar la enseñanza del algoritmo de Transformación Rápida de Fourier de procesamiento de señales a los estudiantes en mis clases de sistemas computacionales. Algunas búsquedas online mostraron el uso de FFT en la multiplicación de grandes enteros y el resto era simplemente historia, como dice el refrán.

La segunda mayor aceleración algorítmica en las pruebas de primalidad modernas llegó más o menos una generación después del algoritmo Schönhage-Strassen, e involucró la longitud de FFT necesaria para multiplicar dos enteros de un determinado tamaño. Para realizar una multiplicación modular de un par de n bits, en general primero debemos calcular exactamente el producto, que es el doble de bits, y luego reducir el producto (es decir, calcular el resto) con respecto al módulo en cuestión. Para los módulos especiales "amigables con los binarios" tales como los números de Mersenne y Fermat, la reducción se puede hacer de manera eficiente utilizando aritmética bit a bit, pero calcular el producto de doble ancho todavía supuen un coste adicional, ya que nos exige aplicar cero a nuestras entradas de complejidad FFT con n 0 bits en la mitad superior y usar una longitud FFT dos veces más grande que nuestros resultados reducidos.

Todo esto cambió en 1994, cuando el fallecido Richard Crandall y Barry Fagin publiscarib un artículo que mostraba cómo, para los casos especiales de Mersenne y módulos de número de Fermat, este relleno de cero podía evitarse mediante la ponderación inteligente de las entradas de transformación con el fin de efectuar una operación de "mod implícito". Este avance produjo una mayor aceleración que el factor 2 en las pruebas de primalidad de estas dos clases de módulos, y los tipos asociados de transformaciones ponderadas discretas se han extendido posteriormente a varias clases interesantes de módulos. Para ver un ejemplo sencillo de cómo se puede usar una FFT para una modificación implícita, remitimos al lector a la post del autor en mersenneforum.org. Desgraciadamente, existen otros aspectos algorítmicos importantes que intervienen en la modulación de alta eficiencia en el hardware de los procesadores modernos: soporte de longitud de transformación no potencial de 2, FFT no recursivas en el lugar que no requieren reordenamiento de reversión de bits antiadherente de los vectores de entrada, el movimiento de datos optimizado para permitir el multihilo eficiente, etc., está más allá del alcance del artículo actual, pero aseguramos al lector/programador interesado que hay más, mucho más, de interés involucrado. Pasamos ahora a las diversas consideraciones que llevaron al esfuerzo especial de implementar soporte de código ensamblador para la plataforma ODROID y su conjunto de instrucciones de aritmética vectorial de 128 bits, y finalmente a los aspectos prácticos de la compilación y ejecución del código de prueba LL resultante en la plataforma Odroid.

¿Por qué la plataforma ARM?

Pasé la mayor parte de mi tiempo de desarrollo en los últimos 5 años paralelizando mi esquema de código FFT de Mlucas existente utilizando el entorno de trabajo pthreads de POSIX más algún código de afinidad central basado en las diversas extensiones compatibles con afinidad proporcionadas en varios sistemas Linuxy MacOS importantes. Una vez que tuve un marco de trabajo paralelo funcional que admitiese desarrollos de vectores SIMD x86 SSE2 de 128 bits y de doble escala, actualicé mis librerias SIMD en línea-ensamblaje-macro a través de sucesivas actualizaciones del conjunto de instrucciones vectoriales de Intel: primero 256 bit AVX, luego AVX2 que en primer lugar promovió las matemáticas del vector entero al mismo nivel que el punto flotante extendiendo el mencionado vector completo de 256 bits y en segundo lugar añadió soporte para instrucciones aritméticas de punto flotante (FMA3) fusionando 3 operandos . Al mismo tiempo que trabajaba en este proyecto reciente, necesitaba ser compatible con el conjunto de instrucciones AVX512 de 512 bits de Intel (que apareció por primera vez en el mercado en forma de barebones indescifrable, pero a la vez muy útil como subconjunto de instrucciones de cimentación desde principios de 2017 en estaciones de trabajo Knights Landing) el año pasado también estaba considerando una primera incursión en agregar soporte SIMD a una familia de procesadores que no sean x86. Las principales consideraciones para emprender tal esfuerzo, que normalmente son de 4 a 5 meses de codificación enfocada, depuración y puesta a punto del rendimiento, fueron las siguientes:

• Amplia Base de instalación y comunidad de desarrolladores activa;
• Compatibilidad con el compilador/depurador estándar de Linux + GCC tolestlow y la paralelización de preads de POSIX;
• Conjunto de instrucciones bien diseñado, preferiblemente al estilo RISC, con al menos tantos registros vectoriales y de propósito general como el Intel AVX’s 16-vector/16-GPR, y preferiblemente tantos registros (32 de cada tipo) como el AVX512 de Intel;
• Soporte para instrucciones FMA;
• Competitividad con las principales familias de procesadores de alta gama (por ejemplo, CPUs Intel y GPUs nVidia) en términos de rendimiento por vatio y de hardware por dólar;
• Posibilidad de mejoras continuas en las implementaciones de procesadores.

ARM (en concreto las instrucciones SIMD-vector ARMv8 de 128 bits y las diversas CPU que lo implementan) pasaron a ser rápidamente el principal candidato. Alta eficiencia energética, un generoso conjunto de registros de 32 + 32 y un excelente diseño de conjunto de instrucciones, mucho mejor que Intel SSE2 con su aspecto de Frankenstein (tuve en cuenta la compatibilidad casi cómica de instrucciones completas de vectores de 64 bits y la falta de una instrucción ROUND en la primera iteración SSE2 como dos de los muchos ejemplos incluidos aquí). Una vez que compré mi pequeño ODROID C2 y trabaje con lo problemas de crecimiento esperados de la nueva nemotecnia de instrucción y la sintaxis de ensamblaje en línea frente a SIMD x86, el desarrollo fue muy sencillo. Una primera inquietud, en la forma de la sintaxis FMA3 algo más restringida frente a la de Intel, resultó infundada, el impacto de dicha restricción se mitigó fácilmente mediante una simple reorganización del orden de los operandos, en relación con el generoso conjunto de 32 registros vectoriales, lo cual fue muy útil. ¡Esperamos que ARM tenga una actualización de 256 bits de las instrucciones del vector v8 en su hoja de ruta para un futuro no muy lejano!

Configurando el software

El software está diseñado para ser tan fácil desarrollar como sea posible en la mayor variedad distribuciones de Linux posibles. Después de haber tenido tantas experiencias malas y de pérdida de tiempo con el típico paradigma de desarrollo configure y make para freeware de Linux, decidí deliberadamente abandonarlo de una manera un tanto radical, en lugar de esforzarme al máximo en la creación de mis archivos cabezera para lo más lejos posible automatizar el proceso de identificación de los aspectos clave de la plataforma de compilación durante el preprocesamiento, tomando cualquiera de los identificadores específicos del hardware admitidos por el software (p. ej., el usuario quiere desarrollar para CPUs x86 que soporte el subconjunto de instrucción vectorial AVX2 / FMA3, o más pertinentemente para ODROIDers, las instrucciones vecto aritmética SIMD ARMv8 de 128 bits) del usuario en tiempo de compilación, además de elegir si se compila un binario de subproceso único o multiproceso.

Así pues el desarrollo se reduce a la siguiente secuencia de 3 pasos, que se presenta en la  Página de Mlucas: Insertando todos los indicadores específicos de las arquitectura aplicables, compila todos los archivos. En mi ODROID-C2, quería activar el ensamblador online que apunta a ARMv8 y quiero poder ejecutarlo en paralelo en los 4 núcleos de procesador, por lo que el comando de compilación es:

$ gcc -c -O3 -DUSE_ARM_V8_SIMD -DUSE_THREADS ../src/*.c >& build.log

Usar grep en el archivo build.log resultante del paso de compilación para los errores:

$ grep -i error build.log

Si no hay errores de compilación, enlaza un archivo binario:

$ gcc -o Mlucas *.o -lm -lpthread -lrt

Ejecuta la típica serie de autocomprabacion para el rango "mediano" de longitudes FFT que cubre todas las asignaciones de usuarios GIMPS actuales y no muy lejanas: './Mlucas -s m'. Un ODROID pasará varias horas probando las diversas formas de combinar las bases aritméticas complejas individuales que componen cada una de las distintas longitudes de FFT cubiertas por la autocomprobación, y capturando las que ofrecen el mejor rendimiento en el hardware particular del usuario a un archivo de configuración maestro llamado mlucas.cfg, que es texto sin formato al margen del sufijo no estándar.

Los tipos más comunes de errores de compilación encontrados anteriormente se pueden rastrear para usuarios que se dirigen a un conjunto de instrucciones no admitido por su CPU particular o ocasionalmente a alguna particularidad de la plataforma de compilación que necesita algunos ajustes para la lógica del preprocesador en el crucial archivo de cabecera plataforma.h. Dado que SBC ODROID viene precargado con distribuciones uniformes de sistema operativo Linux, los únicos problemas que probablemente surjan están relacionados con el conjunto de instrucciones precisas (SIMD v8 o no) compatibles con el modelo particular de ODROID que se utiliza.

En mi ODROID-C2, ejecutando la compilación SIMD de 128 bits del software en los 4 núcleos obtiene un rendimiento total equivalente al de una compilación utilizando las instrucciones vectoriales x86 SSE2 de 128 bits en un único núcleo de un ordenador portátil basada en Intel Core2 a 2 GHz . Esto se traduce en aproximadamente 1/20 del rendimiento total alcanzable en el quad Intel Haswell del autor utilizando el conjunto de instrucciones AVX2 de 256 bits, por lo que los 2 tipos de plataformas no son realmente comparables partiendo del rendimiento por núcleo. Las plataformas basadas en ARM como ODROID brillan por si solas en términos de rendimiento por vatio-hora y por rendimiento de hardware. A efectos de un usuario de ODROID el software tiene una posibilidad realista de descubrir el próximo primo Mersenne, bueno, necesitamos compensar con números lo que las plataformas de hardware de gama más alta (CPU Intel, GPU nVidia) hacen con grandes cantidades de transistores y núcleos cargados de potencia, es decir, tenemos que ser hormigas del ejército para leones. Y, por supuesto, las actualizaciones de velocidad como la prometida por el recientemente anunciado ODROID-N1 pueden ayudar. ¡Es un momento emocionante para ser micro-informático!