PFS: Multiplicación de Karatsuba

Vamos a seguir con nuestra serie sobre 'python, fuente de sabiduria' hablando sobre la multiplicación.

¿Te has preguntado alguna vez cómo se multiplican dos números en un ordenador o calculadora? Normalmente, los circuitos de los procesadores incluyen operaciones básicas como la multiplicación de enteros de tamaño pequeño. Sobre estas operaciones básicas se construye todo lo demás.

El caso que nos ocupa hoy, multiplicar números enteros muy grandes, se consigue mediante el uso de algoritmos vía software. Una forma sencilla de calcular enteros es usando la misma metodología que usábamos en la escuela, con la tabla de multiplicar del uno al diez y con sumas:

    23958233
        5830 ×
------------
    00000000 ( =      23,958,233 ×     0)
   71874699  ( =      23,958,233 ×    30)
 191665864   ( =      23,958,233 ×   800)
119791165    ( =      23,958,233 × 5,000)
------------
139676498390 ( = 139,676,498,390        )

(ejemplo extraído de wikipedia)

El algoritmo que hace uso de esa metodología se suele llamar long, grade-school, naïve,..., multiplication. En este caso vamos a hacer una implementación sucia de este algoritmo descomponiendo el 'multiplicador' en dígitos que multiplicaremos al 'multiplicando':

def naive(x, y):
    total = 0
    yy = [int(i) for i in str(y)]
    for i, yyy in enumerate(yy[::-1]):
        total += x * yyy * (10 ** i)
    return total

Para el ejemplo anterior (multiplicando = 23958233 y multiplicador = 5830) lo que estariamos haciendo sería lo siguiente:

  • el multiplicador lo descomponemos en una lista de los dígitos que lo componen, [5, 8, 3, 0]
  • el 0 lo multiplicamos por el multiplicando
  • el 3 lo multiplicamos por 10 y por el multiplicando y lo sumamos al resultado del paso anterior
  • el 8 lo multiplicamos por 100 y por el multiplicando y lo sumamos al resultado del paso anterior
  • el 5 lo multiplicamos por 1000 y por el multiplicando y lo sumamos al resultado del paso anterior
  • devolvemos el resultado final

Esta es una forma de hacer multiplicaciones pero cuando los números involucrados en la multiplicación son muy grandes podemos encontrarnos con problemas de rendimiento y perder mucho tiempo haciendo un cálculo.

Como en el desarrollo de Python hay gente muy inteligente, si buceas por el código fuente de Python verás que para multiplicación de long's se usa el algoritmo de la multiplicación rápida de Karatsuba. ¿Y quién demonios es este japonés? Pues este japonés es un ruso al que, con 23 añitos (en 1960), se le ocurrió llevar la contraria al mismísimo Kolmogorov.

Lo mejor es que le echéis un ojo a la siguiente entrada (o en la Wikipedia en español) para entender más claramente el funcionamiento (DRY). Básicamente, de lo que se trata es de minimizar el número de multiplicaciones (aunque no siempre se consigue) usando una estrategia de 'divide y vencerás'.

A continuación os dejo una implementación muy compacta que he encontrado en project nayuki (si no entendéis algo del código preguntad en los comentarios):

## http://nayuki.eigenstate.org/res/karatsuba-multiplication/karatsuba.py
## http://nayuki.eigenstate.org/page/karatsuba-multiplication
_CUTOFF = 1024
def kara(x, y):
    if x.bit_length() <= _CUTOFF or y.bit_length() <= _CUTOFF:  # Base case
        #print('naive')
        return naive(x, y)
    else:
        #print('karatsuba')
        n = max(x.bit_length(), y.bit_length())
        half = (n + 32) // 64 * 32
        mask = (1 << half) - 1
        xlow = x & mask
        ylow = y & mask
        xhigh = x >> half
        yhigh = y >> half
        a = kara(xhigh, yhigh)
        b = kara(xlow + xhigh, ylow + yhigh)
        c = kara(xlow, ylow)
        d = b - a - c
        return (((a << half) + d) << half) + c

Bien, pues vamos a ver si esto es capar de optimizar mi función naive. Para ello vamos a multiplicar el siguiente número (tened cuidado si vuestro ordenador no es muy potente ya que el cálculo se puede llegar a alargar):

num = 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890*(10**14000)
print(num)

1234567890123456789012345678901234567890123456789012345678901234567890
1234567890123456789012345678900000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000000000
000000000000000000000000000000

( Qué post más largo... :-) )

Usando la función naive.

a = naive(num, num)
%timeit naive(num, num)

1 loops, best of 3: 1.31 s per loop

Usando la función kara

b = kara(num, num)
%timeit kara(num, num)

10 loops, best of 3: 172 ms per loop

print(a == b)
print(1.31 / 0.173)

True
7.572254335260117

Vemos que hemos conseguido una ganancia de 7.5x!!!

Aunque si comprobamos como lo hace Python vemos que no hemos hecho una gran optimización (aunque no era lo importante de esta entrada):

%timeit num * num

1000 loops, best of 3: 455 µs per loop

Parece que solo es alrededor de 375 veces más rápido ;-)

Este código lo he probado en Python 3.3, si veis alguna errata o alguna barbaridad que haya podido decir, por favor, avísadme.

Esta entrada la puedes descargar en formato notebook desde nuestro repositorio en github.

P.D.: Pretendía que esta entrada fuera más rigurosa y completa pero mi tiempo se ha visto disminuido de forma drástica últimamente. La familia manda :-)

PFS: Exponenciación binaria

Hoy voy a iniciar una serie nueva de entradas que se llamará 'python, fuente de sabiduría' (PFS para abreviar). ¿Qué pretende esta serie? La idea básica es mostrar la influencia y realimentación entre las matemáticas y la programación. La idea e inspiración para estas entradas salió después de leer esta entrada.

Inauguraremos esta entrada con una forma de calcular la potencia entera de un número entero. Si tuviera que crear una función que calculase esto, la primera idea que se me ocurre es usar la fuerza bruta (soy así de primario)...

def pow_fuerzabruta(x,n):
    y = x
    for i in range(1,n):
        y = y * x
    return y

Con la anterior función, la cual no he pensado mucho, lo que estamos haciendo son n-1 multiplicaciones. A medida que aumente el exponente aumentará el tiempo de cálculo. Veámoslo en una gráfica:

from matplotlib import pyplot as plt
n = [5, 25, 100, 1000, 10000, 100000]       # exponente
t = [1.13, 3.93, 17.8, 330, 15000, 1260000] # t en microsegundos
plt.xscale('log')
plt.yscale('log')
plt.plot(n,t)
plt.xlabel('valor del exponente n')
plt.ylabel('t en microsegundos')

pow_bruto

[Nota: la escala de los ejes en la anterior figura es logarítmica]

Esto, cuando usamos números muy grandes, no es muy deseable. Por ello, vamos a mostrar el algortimo de hoy, la exponenciación binaria, también llamado algoritmo de izquierda a derecha binario, que nos permitirá reducir de forma drástica el número de operaciones a realizar. Su uso en Python lo podéis ver, por ejemplo, en la siguiente implementación de long_pow (ver a partir de línea 3376).

Podemos pensar que estas cosas son muy modernas y sofisticadas pero, en este caso, no es así. Los primeros orígenes de esta implementación para calcular potencias parece que sería anterior a 2200 años atrás. Según el libro de Donald Knuth titulado Seminumerical Algorithms, volume 2 of The Art of Computer Programming, página 441, podemos leer:

El método es bastante antiguo; apareció antes del año 200 A.C. en Pingala's Hindu classic Chandah-sutra [ver B. Datta y A.N. Singh, History of Hindu Mathematics 1, 1935]; sin embargo, parece que existen otras referencias a este método fuera de la India a lo largo de los 1000 años posteriores. Una discusión muy clara sobre como calcular 2n de forma eficiente para valores de n arbitrarios sería debida a al-Uqlidisi de Damscus en 952 D.C.; podéis ver The Arithmetic of al-Uglidisi por A.S. Saidan (1975), p. 341-342, donde las ideas generales están ilustradas para el caso de n=51. Ver también al-Biruni's Chronology of Ancient Nations (1879), p. 132-136; durante el siglo XI, los avances del mundo árabe tuvieron una gran influencia.

Más información histórica sobre los orígenes la podéis encontrar en este artículo.

Vamos a ver como funciona, más o menos, este algoritmo. La idea básica detrás del mismo es la siguiente (ver wikipedia para una explicación más completa):

Las potencias más rápidas de calcular son las que tienen la siguiente forma \({x^2}^n\). Estos es debido a que se pueden encontrar las potencias calculando el cuadrado del número de forma repetida.

Vamos a explicar el mismo ejemplo que al-Uqlidisi de Damscus (unos diez siglos más tarde...) considerando un número elevado a 51. Por tanto, puedo calcular \(5^2\) con solo una multiplicación, \(5^4={5^2}^2\) con solo dos multiplicaciones, \(5^8={{5^2}^2}^2\) en tres multiplicaciones,...

Creo que vais viendo por donde van los tiros. El algoritmo explicado de forma muy básica sería para el ejemplo \(5^{51}\):

  • Obtener el exponente como suma de potencias de 2

En este caso sería 51 = 32 + 16 + 2 + 1

  • Encontrar la base de potencias de dos que aparecen en el exponente

\(5^1=5\)

\(5^2 = 5 cdot 5 = 25\)

\(5^4= 25 cdot 25 = 625\)

\(5^8= 625 cdot 625 = 390625\)

\(5^{16} = 390625 cdot 390625 = 152587890625\)

\(5^{32} = 152587890625 cdot 152587890625 = 23283064365386962890625\)

  • Y, por último, se multiplicarían de forma conjunta las potencias ya encontradas

\(5^{51} = 5^{32} cdot 5^{16} cdot 5^2 cdot 5 = 444089209850062616169452667236328125\)

De esta forma, en lugar de multiplicar 5 por si mismo 50 veces solo hemos de realizar 8 multiplicaciones.

El que vamos a implementar es un algoritmo de exponenciación binaria pero de derecha a izquierda. Podéis ver los detalles en la wikipedia.

def my_pow(x,n):
    r = 1
    y = x
    while n > 1:
        if n%2:
            r = r * y
        n = int(n/2)
        y = y * y
    r = r*y
    return r

Si calculamos los tiempos que hemos obtenido con la nueva función (my_pow) y los comparamos con la anterior implementación (pow_fuerzabruta) veremos la ganancia que hemos obtenido:

from matplotlib import pyplot as plt
n = [5, 25, 100, 1000, 10000, 100000]       # exponente
t_bruta = [1.13, 3.93, 17.8, 330, 15000, 1260000] # t  para la primera función
t_bin = [1.34, 2.52, 3.68, 19.8, 765, 32800] # t para la segunda función
plt.plot(n,t_bruta, label='bruta')
plt.plot(n,t_bin, label='bin')
plt.legend()
plt.xlabel('valor del exponente n')
plt.ylabel('t en microsegundos')

pow_bruto_bin

[Nota: Ahora, las escalas de los ejes son lineales]

Vemos que la ganancia, dependiendo del valor del exponente, es cada vez más importante a medida que aumenta el valor del mismo.

for i, t1, t2 in zip(n, t_bruta, t_bin):
    print('Ganancia para n = {0:>6}, {1:>10.2f}'.format(i, t1/t2))


Ganancia para n = 5, 0.84
Ganancia para n = 25, 1.56
Ganancia para n = 100, 4.84
Ganancia para n = 1000, 16.67
Ganancia para n = 10000, 19.00
Ganancia para n = 100000, 38.00

Para números pequeños vemos que incluso puede llegar a ser peor pero para números grandes la ganancia es significativa llegando a valores en torno a 40 para un número n de orden \(10^5\)

Por último, solo comentar que las funciones `pow` y `math.pow` son más eficientes que lo que he hecho (programadas en C). La idea de la entrada solo es mostrar que, gracias al trabajo y saber colectivo, la inteligencia que se esconde en las tripas de Python (y el resto de lenguajes de programación) es mucha y en el día a día ni nos damos cuenta de ello.

P.D.: Las funciones implementadas no contemplan exponentes negativos ni 0. He intentado simplificar el código al máximo para que se vea el bosque en lugar de los árboles.

Esta entrada la puedes descargar en formato notebook desde nuestro repositorio en github.