El algoritmo de Schönhage-Strassen es un algoritmo de multiplicación asintóticamente rápido para números enteros grandes . Fue desarrollado por Arnold Schönhage y Volker Strassen en 1971. [1] La complejidad de bits en tiempo de ejecución es, en notación Big O ,para dos números de n dígitos. El algoritmo utiliza transformadas rápidas de Fourier recursivas en anillos con 2 n +1 elementos, un tipo específico de transformada teórica de números .
El algoritmo de Schönhage-Strassen fue el método de multiplicación asintóticamente más rápido conocido desde 1971 hasta 2007, cuando se anunció un nuevo método, el algoritmo de Fürer , con menor complejidad asintótica; [2] sin embargo, el algoritmo de Fürer actualmente solo logra una ventaja para valores astronómicamente grandes y no se usa en la práctica (ver algoritmos galácticos ).
En la práctica, el algoritmo de Schönhage-Strassen comienza a superar a los métodos más antiguos como la multiplicación de Karatsuba y Toom-Cook para números más allá de 2 2 15 a 2 2 17 (10,000 a 40,000 dígitos decimales). [3] [4] [5] La biblioteca GNU Multi-Precision Library la usa para valores de al menos 1728 a 7808 palabras de 64 bits (33,000 a 150,000 dígitos decimales), dependiendo de la arquitectura. [6] Existe una implementación Java de Schönhage – Strassen que lo usa por encima de 74.000 dígitos decimales. [7]
Las aplicaciones del algoritmo de Schönhage-Strassen incluyen empirismo matemático , como Great Internet Mersenne Prime Search y aproximaciones informáticas de π , así como aplicaciones prácticas como la sustitución de Kronecker , en la que la multiplicación de polinomios con coeficientes enteros se puede reducir eficientemente a un entero grande. multiplicación; GMP-ECM lo utiliza en la práctica para la factorización de la curva elíptica de Lenstra . [8]
Detalles
Esta sección explica en detalle cómo se implementa Schönhage – Strassen. Se basa principalmente en una descripción general del método de Crandall y Pomerance en sus números primos: una perspectiva computacional . [9] Esta variante difiere algo del método original de Schönhage en que explota la transformada ponderada discreta para realizar convoluciones negacyclic de manera más eficiente. Otra fuente de información detallada es Knuth 's The Art of Computer Programming . [10] La sección comienza discutiendo los componentes básicos y concluye con una descripción paso a paso del algoritmo en sí.
Convoluciones
Suponga que estamos multiplicando dos números como 123 y 456 usando una multiplicación larga con dígitos de base B , pero sin realizar ningún acarreo. El resultado podría verse así:
1 | 2 | 3 | ||
× | 4 | 5 | 6 | |
6 | 12 | 18 | ||
5 | 10 | 15 | ||
4 | 8 | 12 | ||
4 | 13 | 28 | 27 | 18 |
Esta secuencia (4, 13, 28, 27, 18) se denomina convolución acíclica o lineal de las dos secuencias originales (1, 2, 3) y (4, 5, 6). Una vez que tenemos la convolución acíclica de dos secuencias, calcular el producto de los números originales es fácil: simplemente realizamos el transporte (por ejemplo, en la columna más a la derecha, teníamos que mantener el 8 y agregar el 1 a la columna que contiene 27). En el ejemplo, esto da como resultado el producto correcto 56088.
Hay otros dos tipos de convoluciones que serán útiles. Suponga que las secuencias de entrada tienen n elementos (aquí 3). Entonces la convolución acíclica tiene n + n −1 elementos; si tomamos los n elementos más a la derecha y agregamos los n −1 elementos más a la izquierda , esto produce la convolución cíclica :
28 | 27 | 18 | ||
+ | 4 | 13 | ||
28 | 31 | 31 |
Si realizamos la realización de la convolución cíclica, el resultado es equivalente al producto de las entradas mod B n - 1. En el ejemplo, 10 3 - 1 = 999, la realización de la realización de (28, 31, 31) produce 3141, y 3141 ≡ 56088 (mod 999).
Por el contrario, si tomamos los n elementos más a la derecha y restamos los n −1 elementos más a la izquierda , esto produce la convolución cíclica negacy :
28 | 27 | 18 | ||
- | 4 | 13 | ||
28 | 23 | 5 |
Si realizamos la realización de la convolución negacyclic, el resultado es equivalente al producto de las entradas mod B n + 1. En el ejemplo, 10 3 + 1 = 1001, la realización de la realización de (28, 23, 5) produce 3035 y 3035 ≡ 56088 (mod 1001). La convolución negacyclic puede contener números negativos, que pueden eliminarse durante el transporte utilizando préstamos, como se hace en la resta larga.
Teorema de convolución
Al igual que otros métodos de multiplicación basados en la transformada rápida de Fourier , Schönhage-Strassen depende fundamentalmente del teorema de convolución , que proporciona una forma eficiente de calcular la convolución cíclica de dos secuencias. Se afirma que:
- La convolución cíclica de dos vectores se puede encontrar tomando la transformada discreta de Fourier (DFT) de cada uno de ellos, multiplicando los vectores resultantes elemento por elemento y luego tomando la transformada discreta inversa de Fourier (IDFT).
O en símbolos:
- Convolución cíclica ( X, Y ) = IDFT (DFT ( X ) · DFT ( Y ))
Si calculamos la DFT y la IDFT usando un algoritmo de transformada rápida de Fourier e invocamos nuestro algoritmo de multiplicación de forma recursiva para multiplicar las entradas de los vectores transformados DFT ( X ) y DFT ( Y ), esto produce un algoritmo eficiente para calcular la convolución cíclica.
En este algoritmo, será más útil calcular la convolución negacyclic ; como resulta, una versión ligeramente modificada del teorema de convolución (ver transformada ponderada discreta ) también puede permitir esto. Suponga que los vectores X e Y tienen longitud n , y a es una raíz primitiva de la unidad de orden 2 n (es decir, a 2 n = 1 y a para todas las potencias menores no es 1). Entonces podemos definir un tercer vector A , llamado vector de peso , como:
- A = ( una j ), 0 ≤ j < n
- A −1 = ( a - j ), 0 ≤ j < n
Ahora, podemos afirmar:
- NegacyclicConvolution ( X , Y ) = A −1 · IDFT (DFT ( A · X ) · DFT ( A · Y ))
En otras palabras, es lo mismo que antes, excepto que las entradas se multiplican primero por A y el resultado se multiplica por A −1 .
Elección de anillo
La transformada discreta de Fourier es una operación abstracta que se puede realizar en cualquier anillo algebraico ; normalmente se realiza en números complejos, pero en realidad realizar aritmética compleja con la precisión suficiente para garantizar resultados precisos para la multiplicación es lento y propenso a errores. En su lugar, vamos a utilizar el enfoque de la cantidad teórica transformar , que es llevar a cabo la transformación de los enteros mod N para algún entero N .
Al igual que existen raíces primitivas de unidad de cada orden en el plano complejo, dado cualquier orden n podemos elegir una N adecuada tal que b sea una raíz primitiva de unidad de orden n en los enteros mod N (en otras palabras, b n ≡ 1 (mod N ), y ninguna potencia menor de b es equivalente a 1 mod N ).
El algoritmo dedicará la mayor parte de su tiempo a realizar multiplicaciones recursivas de números más pequeños; con un algoritmo ingenuo, estos ocurren en varios lugares:
- Dentro del algoritmo de transformada rápida de Fourier, donde la raíz primitiva de la unidad b se alimenta, eleva al cuadrado y multiplica repetidamente por otros valores.
- Al tomar potencias de la raíz primitiva de la unidad a para formar el vector de peso A y al multiplicar A o A −1 por otros vectores.
- Al realizar la multiplicación elemento por elemento de los vectores transformados.
La idea clave de Schönhage-Strassen es elegir N, el módulo, para que sea igual a 2 n + 1 para algún número entero n que sea un múltiplo del número de piezas P = 2 p que se están transformando. Esto tiene una serie de beneficios en los sistemas estándar que representan números enteros grandes en forma binaria:
- Cualquier valor puede reducirse rápidamente módulo 2 n + 1 usando solo cambios y adiciones, como se explica en la siguiente sección .
- Todas las raíces de unidad utilizadas por la transformada se pueden escribir en la forma 2 k ; en consecuencia, podemos multiplicar o dividir cualquier número por una raíz de unidad usando un desplazamiento, y potenciar o elevar al cuadrado una raíz de unidad operando solo en su exponente.
- Las multiplicaciones recursivas elemento por elemento de los vectores transformados se pueden realizar mediante una convolución negacyclic, que es más rápida que una convolución acíclica y ya tiene "gratis" el efecto de reducir su resultado mod 2 n + 1.
Para hacer convenientes las multiplicaciones recursivas, enmarcaremos Schönhage-Strassen como un algoritmo de multiplicación especializado para calcular no solo el producto de dos números, sino el producto de dos números mod 2 n + 1 para un n dado . Esto no es una pérdida de generalidad, ya que siempre se puede elegir n lo suficientemente grande como para que el producto mod 2 n + 1 sea simplemente el producto.
Optimizaciones de turno
En el curso del algoritmo, hay muchos casos en los que la multiplicación o la división por una potencia de dos (incluidas todas las raíces de la unidad) pueden ser reemplazadas de manera rentable por una pequeña cantidad de cambios y sumas. Esto hace uso de la observación de que:
- (2 norte ) k ≡ (−1) k mod (2 norte + 1)
Tenga en cuenta que un número de k dígitos en base 2 n escrito en notación posicional se puede expresar como ( d k −1 , ..., d 1 , d 0 ). Representa el numero. También tenga en cuenta que para cada d i tenemos 0 ≤ d i <2 n .
Esto simplifica la reducción de un número representado en el modo binario 2 n + 1 : tome los n bits más a la derecha (menos significativos) , reste los n bits siguientes , agregue los n bits siguientes , y así sucesivamente hasta que se agoten los bits. Si el valor resultante aún no está entre 0 y 2 n , normalícelo sumando o restando un múltiplo del módulo 2 n + 1 . Por ejemplo, si n = 3 (y entonces el módulo es 2 3 +1 = 9) y el número que se reduce es 656, tenemos:
- 656 = 1010010000 2 ≡ 000 2 - 010 2 + 010 2 - 1 2 = 0 - 2 + 2 - 1 = −1 ≡ 8 (mod 2 3 + 1).
Además, es posible efectuar cambios muy grandes sin construir nunca el resultado cambiado. Supongamos que tenemos un número A entre 0 y 2 n , y deseamos multiplicarlo por 2 k . Dividiendo k entre n encontramos k = qn + r con r < n . Resulta que:
- A (2 k ) = A (2 qn + r ) = A [(2 n ) q (2 r )] ≡ (−1) q Shift_Left (A, r ) (mod 2 n + 1).
Dado que A es ≤ 2 n y r < n , un desplazamiento a la izquierda r tiene como máximo 2 n −1 bits, por lo que solo se necesita un desplazamiento y una resta (seguido de normalización).
Finalmente, para dividir por 2 k , observe que al elevar al cuadrado la primera equivalencia anterior se obtiene:
- 2 2 norte ≡ 1 (mod 2 norte + 1)
Por eso,
- A / 2 k = A (2 - k ) ≡ A (2 2 n - k ) = Desplazamiento_Izquierda (A, 2 n - k ) (mod 2 n + 1).
Algoritmo
El algoritmo sigue una división, evaluación (FFT directa), multiplicación puntual, interpolación (FFT inversa) y combina fases similares a los métodos Karatsuba y Toom-Cook.
Números dados de entrada x y y , y un número entero N , el siguiente algoritmo calcula el producto xy mod 2 N + 1. Siempre N es suficientemente grande esto es simplemente el producto.
- Dividir cada número de entrada en vectores de X e Y de 2 k piezas de cada uno, donde 2 k divisiones N . (por ejemplo, 12345678 → (12, 34, 56, 78)).
- Para progresar, es necesario utilizar una N más pequeña para las multiplicaciones recursivas. Para este propósito, elija n como el número entero más pequeño al menos 2 N / 2 k + k y divisible por 2 k .
- Calcule el producto de X e Y mod 2 n + 1 usando la convolución negacyclic:
- Multiplica X e Y cada uno por el vector de peso A usando cambios (desplaza la j- ésima entrada a la izquierda por jn / 2 k ).
- Calcule la DFT de X e Y usando la FFT de la teoría numérica (realice todas las multiplicaciones usando turnos; para la 2 k -ésima raíz de la unidad, use 2 2 n / 2 k ).
- Aplique de forma recursiva este algoritmo para multiplicar los elementos correspondientes de las X e Y transformadas.
- Calcule la IDFT del vector resultante para obtener el vector de resultado C (realice todas las multiplicaciones usando cambios). Esto corresponde a la fase de interpolación.
- Multiplique el vector de resultado C por A −1 usando cambios.
- Ajustar signos: algunos elementos del resultado pueden ser negativos. Calculamos el valor positivo más grande posible para el j- ésimo elemento de C, ( j + 1) 2 2 N / 2 k , y si lo excede restamos el módulo 2 n + 1.
- Por último, realice el transporte mod 2 N + 1 para obtener el resultado final.
El número óptimo de piezas para dividir la entrada es proporcional a , donde N es el número de bits de entrada, y esta configuración alcanza el tiempo de ejecución de O ( N log N log log N ), [1] [9] por lo que el parámetro k debe configurarse en consecuencia. En la práctica, se establece empíricamente en función de los tamaños de entrada y la arquitectura, por lo general a un valor entre 4 y 16. [8]
En el paso 2, se utiliza la observación de que:
- Cada elemento de los vectores de entrada tiene como máximo n / 2 k bits;
- El producto de cualesquiera dos elementos del vector de entrada tiene como máximo 2 n / 2 k bits;
- Cada elemento de la convolución es la suma de como máximo 2 k de tales productos, por lo que no puede exceder de 2 n / 2 k + k bits.
- n debe ser divisible por 2 k para garantizar que en las llamadas recursivas se cumpla la condición "2 k divide N " en el paso 1.
Optimizaciones
Esta sección explica una serie de optimizaciones prácticas importantes que se han considerado al implementar Schönhage – Strassen en sistemas reales. Se basa principalmente en un trabajo de 2007 de Gaudry, Kruppa y Zimmermann que describe las mejoras de la biblioteca de precisión múltiple de GNU . [8]
Por debajo de cierto punto de corte, es más eficiente realizar las multiplicaciones recursivas utilizando otros algoritmos, como la multiplicación de Toom-Cook . Los resultados deben reducirse mod 2 n + 1, lo que se puede hacer de manera eficiente como se explicó anteriormente en Optimizaciones de turnos con turnos y sumas / restas.
Calcular la IDFT implica dividir cada entrada por la raíz primitiva de la unidad 2 2 n / 2 k , una operación que frecuentemente se combina con la multiplicación del vector por A −1 posteriormente, ya que ambas implican la división por una potencia de dos.
En un sistema donde un gran número se representa como una matriz de palabras de 2 w bits, es útil asegurarse de que el tamaño del vector 2 k sea también un múltiplo de los bits por palabra eligiendo k ≥ w (por ejemplo, elija k ≥ 5 en una computadora de 32 bits yk ≥ 6 en una computadora de 64 bits); esto permite que las entradas se dividan en partes sin cambios de bits, y proporciona una representación uniforme para los valores mod 2 n + 1 donde la palabra alta solo puede ser cero o uno.
La normalización implica sumar o restar el módulo 2 n + 1; este valor tiene solo dos bits establecidos, lo que significa que esto se puede hacer en tiempo constante en promedio con una operación especializada.
Los algoritmos iterativos de FFT como el algoritmo de Cooley-Tukey FFT , aunque se utilizan con frecuencia para FFT en vectores de números complejos, tienden a exhibir una localidad de caché muy pobre con las entradas de vectores grandes utilizadas en Schönhage-Strassen. La implementación recursiva directa, no en el lugar de FFT es más exitosa, con todas las operaciones encajando en la caché más allá de un cierto punto en la profundidad de la llamada, pero aún hace un uso subóptimo de la caché en profundidades de llamada más altas. Gaudry, Kruppa y Zimmerman utilizaron una técnica que combina el algoritmo de 4 pasos de Bailey con transformaciones de radix superiores que combinan múltiples pasos recursivos. También mezclan fases, yendo lo más lejos posible del algoritmo en cada elemento del vector antes de pasar al siguiente.
El "truco de la raíz cuadrada de 2", descrito por primera vez por Schönhage, es notar que, siempre que k ≥ 2, 2 3 n / 4 −2 n / 4 sea una raíz cuadrada de 2 mod 2 n +1, y por lo tanto un 4 n -ésima raíz de la unidad (ya que 2 2 n ≡ 1). Esto permite extender la longitud de la transformada de 2 k a 2 k + 1 .
Finalmente, los autores tienen cuidado de elegir el valor correcto de k para diferentes rangos de números de entrada, notando que el valor óptimo de k puede ir y venir entre los mismos valores varias veces a medida que aumenta el tamaño de entrada.
Referencias
- ↑ a b A. Schönhage y V. Strassen, " Schnelle Multiplikation großer Zahlen ", Computación 7 (1971), págs. 281-292.
- ^ Martin Fürer, " Multiplicación de enteros más rápida Archivado el25 de abril de 2013en la Wayback Machine ", Actas de STOC 2007, págs. 57–66.
- ^ Van Meter, Rodney; Itoh, Kohei M. (2005). "Exponenciación modular cuántica rápida". Revisión física . A. 71 (5): 052320. arXiv : quant-ph / 0408006 . doi : 10.1103 / PhysRevA.71.052320 . S2CID 14983569 .
- ^ Descripción general de las características de Magma V2.9, sección aritmética. Archivado el 20 de agosto de 2006 en la Wayback Machine : analiza los puntos de cruce prácticos entre varios algoritmos.
- ^ Luis Carlos Coronado García, " ¿Puede la multiplicación de Schönhage acelerar el cifrado o descifrado RSA? Archivado ", Universidad de Tecnología, Darmstadt (2005)
- ^ "MUL_FFT_THRESHOLD" . El rincón de los desarrolladores de GMP . Archivado desde el original el 24 de noviembre de 2010 . Consultado el 3 de noviembre de 2011 .
- ^ "Una clase BigInteger mejorada que utiliza algoritmos eficientes, incluido Schönhage-Strassen" . Oracle . Consultado el 10 de enero de 2014 .
- ^ a b c Pierrick Gaudry, Alexander Kruppa y Paul Zimmermann. Una implementación basada en GMP del algoritmo de multiplicación de enteros grandes de Schönhage-Strassen Archivado . Actas del Simposio Internacional de 2007 sobre Computación Simbólica y Algebraica, pp.167-174.
- ^ a b R. Crandall y C. Pomerance. Números primos: una perspectiva computacional . Segunda edición, Springer, 2005. Sección 9.5.6: Método de Schönhage, p. 502. ISBN 0-387-94777-9
- ^ Donald E. Knuth, El arte de la programación informática , Volumen 2: Algoritmos seminuméricos (3ª edición), 1997. Addison-Wesley Professional, ISBN 0-201-89684-2 . Sección 4.3.3.C: Transformadas discretas de Fourier, pág.305.