Lo que mas chocaría a Gauss si volviera a la vida sería saber que los estudios que hizo para determinar la excentricidad del asteroide Juno (0.254236) se usan hoy día para multiplicar números naturales.
Esa sorpresa la vivió en primera persona Kolmogorov. En los años 50, cuando tenía unos 50 años, influido por los trabajos de Wiener y Claude Shannon, se interesó por los algoritmos de computación, fue de los primeros en interesarse por la complejidad de los algoritmos.
Kolmogorov se planteó el problema de cuantificar el trabajo de efectuar un algoritmo. ¿Cuánto mas cuesta multiplicar dos números de n cifras que sumarlos? La respuesta es contar cuantas operaciones elementales tenemos que realizar. Podríamos dar una definición precisa de qué es una operación elemental, pero creo que en lo que sigue será suficiente las ideas intuitivas que tenemos. Para sumar los dos números tenemos que leer cada cifra (pues al modificar una sola el resultado cambiaría), son 2n operaciones y tenemos que escribir el resultado, otras n operaciones. El trabajo con cada par de cifras es mas o menos constante. De esté modo vemos que el número S(n) de operaciones necesarios para sumar dos números de n dígitos verifica las acotaciones 2n≤S(n)≤Cn donde C es cierta constante. Diremos que el número de operaciones necesarias para la suma es S(n)=O(n). Esto es, del orden de n operaciones.
Andrei N. Kolmogorov
Qué ocurre con el producto. Denotamos por M(n) el orden de operaciones elementales necesarias para multiplicar números de ncifras con el mejor algoritmo posible. Tenemos también que M(n)≥2n porque los números, al menos, hay que leerlos. En el algoritmo que aprendemos en la escuela debemos hacer una operación elemental con cada pareja de dígitos uno del multiplicando y otro del multiplicador. Son n2 parejas. El trabajo restante es de orden menor. De manera que llegamos a que el algoritmo usual prueba que M(n)≤cn2.
El algoritmo usual es muy antiguo. Tanto que en realidad no sabemos cuando se originó. Un algoritmo semejante existía en tiempos de los sumerios o egipcios, hace unos 4000 años, y posiblemente viene de más atrás en el tiempo. Kolmogorov pensaba que si existiese otro algoritmo mejor ya habría sido descubierto. En sus seminarios planteaba la hipótesis de que M(n)=O(n2). Esto es, que el mejor método para multiplicar números de longitud n necesitaría del orden de n2 operaciones elementales.
Anatoly A. Karatsuba
Anatoly Karatsuba, un alumno de Kolmogorov, con 23 años en aquel momento, cuenta en [2] lo que sucedió:
En el otoño de 1960 tuvo lugar un seminario sobre problemas en cibernética en la Facultad de Mecánica y Matemáticas de la Universidad de Moscú bajo la dirección de Kolmogorov. En ella Kolmogorov enunció la conjetura n2 y propuso algunos problemas relativos a la complejidad de la solución de sistemas lineales de ecuaciones y algunos cálculos similares. Comencé a pensar activamente sobre la conjetura n2, y en menos de una semana encontré que el algoritmo con el que pensaba deducir una cota inferior para M(n) proporcionaba una cota de la forma
M(n)=O(nlog23),log23=1.5849…
Al finalizar la siguiente reunión del seminario le conté a Kolmogorov el nuevo algoritmo y la refutación de la conjetura n2. Kolmogorov estaba muy agitado porque esto contradecía su muy plausible conjetura. En la siguiente reunión del seminario, Kolmogorov mismo contó a los participantes sobre mi método y en este punto el seminario acabó.
Son 4000 años vencidos en una semana por una mente joven, abierta, frente a un problema bien propuesto. Es interesante el análisis. Kolmogorov se planteó el problema de determinar el orden de M(n). En los 4000 años anteriores se había buscado un buen método y se había encontrado. Gauss, por ejemplo, no se planteó mejorarlo, era un buen método y rara vez se tenían que multiplicar números de muchos dígitos. Pero Kolmogorov planteó el problema de ver que no hay otro método mejor, probar que M(n)≥cn2. Es decir, propuso probar que no importa el método que usemos necesitaremos del orden de n2 operaciones elementales.
Entonces llegó Karatsuba con su mente abierta y trató de probar la conjetura de Kolmogorov. ¿Cómo hacerlo? Su estrategia fue suponer conocido M(n) y trató de encontrar una relación entre M(n) y M(n/2). Entonces hay otra cuestión. Si sé multiplicar números de n/2 cifras, ¿en que me ayuda esto para multiplicar números de 2n cifras? Como veremos esto lleva directamente al nuevo algoritmo.
Buenas preguntas y mente abierta son dos de los ingredientes de las matemáticas.
El algoritmo de Karatsuba.
El algoritmo de Karatsuba es un ejemplo de la estrategia general divide y conquista. La idea es simple: sean Ay B los números a multiplicar con 2n cifras; los separamos en grupos de n cifras
A=A110n+A0,B=B110n+B0.
entonces
AB=(A110n+A0)(B110n+B0)=A1B1102n+(A1B0+A0B1)10n+A0B0.
Notemos que
A1B0+A0B1=(A1−A0)(B0−B1)+A1B1+A0B0.
De manera que obtenemos
AB=(102n+10n)A1B1+10n(A1−A0)(B0−B1)+(10n+1)A0B0
Hemos reducido así el producto de dos números de 2n dígitos a tres productos de números de n dígitos A1B1, A0B0 y (A1−A0)(B0−B1) más unas cuantas traslaciones y sumas. Luego hemos probado que
M(2n)≤3M(n)+8S(n)≤3M(n)+an,
para cierta constante a. Por inducción probamos que para una potencia de dos n=2m se tiene
M(2m)≤M(1)3m+a(2m+2m−13+2m−232+⋯+2⋅3m−1)
y siendo 2α=3, esto es, α=log23, obtenemos
M(2m)≤3m(M(1)+2a)=C2αm=C(2m)α.
Ahora para cualquier n, existe m tal que 2m−1≤n<2m y cualquier número de n dígitos puede considerarse un número de 2m dígitos, luego
M(n)≤M(2m)≤C(2m)α=C2α(2m−1)α≤3Cnα.
Esto es una versión del método de Karatsuba, simplificando para la exposición el método contenido en [2].
Nuevos progresos, el algoritmo de Schönhage-Strassen.
Después de la sorpresa de Karatsuba quedó planteado el problema del verdadero orden de M(n). Toom y Cook (1966) extendieron las ideas de Karatsuba y consiguieron probar que M(n)=O(n(logn)eclogn√). Nos centraremos en lo que sigue en explicar el método de Schönhage-Strassen (1971) que es el que se usa en la práctica hoy día cuando se quieren multiplicar números realmente grandes, por ejemplo con cientos de miles o hasta billones de dígitos.
Schönhage y Strassen jugando al ajedrez
Todos los métodos conocidos se basan en una idea:
Multiplicar números es equivalente a multiplicar polinomios.
Los factores están escritos en una base, realmente los métodos de multiplicación usan la base 2, pero lo ilustramos con nuestra base decimal que nos es más familiar. La notación 325 por ejemplo significa
3×102+2×10+5=p(10),dondep(x)=3x2+2x+5.
Y multiplicar dos números es casi equivalente a multiplicar dos polinomios. Por ejemplo
2682583297072775253556x46x421x3+4x3+25x33x22x29x2+14x2+10x2+33x2+2x+7x+6x+35x+41x+5+3+1515
El producto de polinomios es mas ordenado, no hay transferencias de una columna a otra. Además es prácticamente equivalente pues el polinomio producto evaluado en x=10 nos da el valor del producto de los números. De hecho conociendo los coeficientes del polinomio producto, se puede obtener el producto, como se ve en el siguiente diagrama.
Los coeficientes del polinomio producto serán menores que nb2 si estamos trabajando en base b, y el número de operaciones elementales para obtener el producto a partir del polinomio será aproximadamente nlog(nb2).
En el método de Schönhage-Strassen si queremos multiplicar números de n dígitos binarios, se buscan números k y ℓ tales que 2n≤2kℓ<4n y entonces trabajamos con números de 2k dígitos en base b=2ℓ. De este modo, obtenido el producto de los polinomios, necesitamos 2klog(2k22ℓ)≈2k(k+2ℓ). Tomando k y ℓ del mismo orden conseguimos obtener el producto con O(n) operaciones elementales, suponiendo conocido el producto de polinomios.
Convolución y Transformada de Fourier.
Recordemos que para multiplicar números de $n$ dígitos escogemos 2n≤2kℓ<4n y trabajamos en base b=2ℓ, es decir agrupando las cifras de los factores en grupos de ℓ dígitos binarios.
Los factores tendrán aproximadamente n/ℓ<2k−1 dígitos en base 2ℓ. Definiendo K=2k podremos escribir los factores en la base b en la forma
(UK−1,UK−2,…,U0)b,(VK−1,VK−2,…,V0)b.
Como eran números de menos de K/2 cifras, tendremos Uj=Vj para j≥K/2. Eso lo hemos hecho intencionadamente por lo que sigue.
Debemos calcular los coeficientes del polinomio producto
P(x)=U(x)V(x)=(UK−1xK−1+UK−2xK−2+⋯+U1x+U0)(VK−1xK−1+VK−2xK−2+⋯+V1x+V0)
que son los números
Pr=UrV0+Ur−1V1+⋯+U0Vr
gracias a que hemos añadido los Uj=Vj=0 para los j grandes se tiene que esto es igual a la convolución en el grupo Z/(KZ) de los números modulo K es decir que
Pr=∑a+b=rUaVb.
La operación de convolución está estrechamente relacionada con la transformada de Fourier. En el caso del grupo finito la trasformada de (Uj)Kj=0 es la sucesión
Uˆa=∑j=0K−1Ujωaj,ω=e−2πi/K.(1)
Es decir
Uˆa=U(ωa),Vˆa=V(ωa)
son los valores de los polinomios U(x) y V(x) en las raíces de la unidad.
De esto se sigue la propiedad fundamental de la transformada y la convolución. Debido a que el polinomio P(x)=U(x)V(x), tendremos que la transformada de la convolución es el producto de las transformadas
P(ωa)=U(ωa)V(ωa),esto esPˆa=UˆaVˆa.
Así que, no los coeficientes Pa, pero sí sus transformados de Fourier se obtienen fácilmente mediante tan solo K productos.
Al rescate viene entonces el hecho de que la inversa de la transformada es esencialmente otra transformada de Fourier
Pa=1K∑j=0K−1Pˆaω−aj.
Tenemos entonces el siguiente esquema para multiplicar números. Encontrar la transformada de Fourier de los factores Uˆa y Vˆa. Formar K productos Pˆa=UˆaVˆa y realizar una nueva transformada de Fourier (inversa) para obtener los coeficientes del producto.
El problema es que tal como hemos definido la transformada en (1), necesitamos K2 productos, y esto es demasiado grande. En este punto nos ayuda una técnica que ya era conocida por Gauss pero que fue redescubierta y popularizada en Cooley y Tukey en 1965.
Transformada de Fourier rápida.
La técnica de la FFT (Fast Fourier Transform), es complicada. Su utilidad práctica se extiende a muchos campos de la técnica. Cualquier procedimiento en que se aplique la transformada de Fourier se aprovecha de ella. Gauss [1] la estudio en relación a sus cálculos astronómicos. Gauss quería determinar la órbita de un asteroide a partir de una serie de mediciones. Expresaba la órbita como una serie trigonométrica y deseaba obtener los coeficientes a partir de las mediciones. Esto es en realidad una transformada de Fourier.
Diagrama de mariposa para la FFT
FFT se basa en la identidad de Danielson-Lanczos. Suponiendo K par, se tiene
Esta igualdad expresa una transformada de Fourier, que requiere K2 productos, mediante dos transformadas de longitud mitad que requieren 2(K/2)2=K2/2 multiplicaciones. Si K es una potencia de dos, podemos iterar el procedimiento y al final obtenemos un método de calcular la transformada que tan solo requiere O(KlogK) operaciones.
Repito que hay libros enteros dedicados a la técnica FFT y que por tanto a nuestra explicación, siendo correcta, le faltan algunos detalles.
El Álgebra y el juego de los anillos.
Llegamos al último ingrediente del método de Schönhage-Strassen. Tal como lo hemos descrito ha sido implementada con éxito. En ella se ha usado el anillo de los números complejos. La raíz de la unidad ω es un número complejo, los cálculos han de hacerse en C y las transformadas Uˆj son números complejos. Esto hace necesario trabajar con números aproximados y requiere un estudio riguroso de los errores, aparte de que hay un poco de suerte en que los errores no se propagan demasiado en los cálculos a realizar.
Pero es posible sustituir el cuerpo C por un anillo en que tengamos una raíz primitiva K-ésima de la unidad y en el que K sea invertible. Con ello podremos definir la transformada sustituyendo el número complejo ωpor la raíz de la unidad. Las fórmulas anteriores siguen siendo válidas.
Por ejemplo uno puede usar el anillo Z((2m+1)Z) tomando m>k+2ℓ. De este modo los números Pamod(2m+1) permiten recuperar Pa ya que Pa<2m.
Tomamos m que sea multiplo de K, de este modo ω=22m/K es una raíz primitiva de la unidad de orden K. La elección dado n de k, ℓ K, y m es uno de los problemas prácticos del método de Schönhage-Strassen. Pero unas buenas elecciones garantizan que el método se realiza con O(nlognloglogn)operaciones elementales.
Avances posteriores sobre el valor de M(n).
La conjetura de Kolmogorov no era válida, y el resultado de Karatsuba lo demuestra. Pero entonces ¿cuál es el verdadero orden de M(n)? No lo sabemos, desde el resultado de Karatsuba se han ido construyendo mejores algoritmos que rebajan la cota superior de M(n). Los avances posteriores al método de Schönhage-Strassen no han tenido repercusión práctica. En primer lugar son solo modificaciones, podemos decir que ligeras, del método de Schönhage-Strassen del que se separan principalmente en la elección de anillo o anillos en los que se trabaja.
Segundo, los avances se refieren al último factor en nlognloglogn. Es sustituido por otro factor que crece mas lentamente, pero dado que loglogn crece tan lento, la ventaja no tendría lugar más que para números que realmente no están al alcance de los ordenadores actuales.
Para explicar hasta donde llegan estas mejoras, debemos definir una función muy lenta, extremadamente lenta log∗x, definida para x≥1. Dado x>1, pongamos x0=x, x1=logx, x2=logx1, …, xn=logxn−1, parando en el momento que encontremos el primer xn≤1. Entonces definimos log∗x=n. La función log∗x vale 1 en el intervalo (1,e], log∗x=2 para x∈(e,ee], log∗x=3 para x∈(ee,eee], … vale 4 cuando x tiene menos de 1656520 dígitos. La función log∗x tiende a infinito pero de la manera mas perezosa imaginable (si no fuéramos matemáticos, que tenemos una imaginación inimaginable).
Ahora estamos en condiciones de explicar la evolución de nuestro conocimiento sobre M(n).
Hemos escrito esta entrada con motivo de la aparición del trabajo de Harvey, y van der Hoeven [11] en arXiv en febrero de este año. En trabajos anteriores se había conseguido la multiplicación en O(nlogn4log∗n)operaciones elementales, pero suponiendo la existencia de ciertos primos que eran necesarios para construir los anillos en que se trabajaba. Esto hacía el resultado condicional, porque aunque los primos adecuados se supone que existen, eso no estaba probado. En este trabajo se obtienen anillos adecuados sin necesidad de hacer hipótesis sobre la existencia de ciertos números primos.
Hay una pregunta que no quiero dejar de formular antes de terminar la entrada.
Problema 1. ¿Es realmente M(n)≥Cnlogn?
Para la suma conocemos el orden exacto ya que es claro que S(n)>Cn, porque al menos tenemos que leer los datos, y el método usual de suma da una cota superior del mismo orden. En nuestra cuestión es necesario considerar cualquier método posible para multiplicar y ver que todos necesitan hacer del orden de nlogn operaciones elementales. Este tipo de problemas son actualmente muy difíciles.
Para saber más.
El aspecto histórico de como Gauss fue el primero en definir la técnica de la FFT se encuentra bien desarrollado en
[1] Heideman, M.T., Johnson, D.H., y Burrus, C. S., Gauss and the History of the Fast Fourier Transform, IEEE ASSP Magazine 1 (4) (1984) 14-21.
Copia en pdf.
La historia de Karatsuba la he sacado de
[2] Karatsuba, A. A., The Complexity of computations, Proc. Steklov Inst. Math., 211 (1995), 169-183.
Pdf del artículo.
Esta historia tiene una segunda parte, un tanto desgraciada. Karatsuba lo cuenta, pero podemos ver algunos detalles y/o comentarios en:
[3] Pregunta sobre la relación Kolmogorov-Karatsuba en Stackexchange sobre teoría de la computación.
El método de Schönhage-Strassen está expuesto en diversas fuentes. Recomendamos especialmente
[4] Knuth, D. E. The art of computer programming. Vol. 2. Seminumerical algorithms, Third edition. Addison-Wesley, Reading, MA, 1998.
o también
[5] Crandall, R., Pomerance, C., Prime numbers. A computational perspective, Second edition. Springer, New York, 2005.
El trabajo original está escrito en Alemán:
[6] Schönhage, A.; Strassen, V., Schnelle Multiplikation grosser Zahlen, Computing (Arch. Elektron. Rechnen) 7(1971), 281-292.
Para el que desee información en detalle puede ser útil
[7] Gaudry, P., Kruppa, A., Zimmermann, P., A GMP-based implementation of Schönhage-Strassen’s large integer multiplication algorithm, ISSAC 2007, 167–174, ACM, New York, 2007. Pdf del artículo.
Otros trabajos que he consultado para esta entrada:
[8] Pollard, J. M., The Fast Fourier Transform in a Finite Field, Math. of Comp. 25, (1971) 365-374. Pdf del artículo.
[9] Toom, A. L., The complexity of a scheme of functional elements simulating the multiplication of integers, Soviet Math. Dokl., 3 (1963) 714–716.
El pdf del artículo puede encontrase en su página web con algún otro material interesante.
[10] Fürer, M., Faster integer multiplication, STOC’07—Proceedings of the 39th Annual ACM Symposium on Theory of Computing, 57–66, ACM, New York, 2007.
Pdf del artículo.
[11] Harvey, D. y J. van der Hoeven, Faster Integer multiplication using short lattice vectors. febrero 2018, arXiv:1802.07932.
