Hoy voy a iniciar una serie nueva de entradas que se llamará 'python, fuente de sabiduría'. ¿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)...
In [45]:
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:
In [61]:
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')
Out[61]:
[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
2nde forma eficiente para valores denarbitrarios 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, puesdo 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}$:
En este caso sería $51 = 32 + 16 + 2 + 1$
$5^1 = 5$
$5^2 = 5 · 5 = 25$
$5^4 = 25 · 25 = 625$
$5^8 = 625 · 625 = 390625$
$5^{16} = 390625 · 390625 = 152587890625$
$5^{32} = 152587890625 · 152587890625 = 23283064365386962890625$
$5^{51} = 5^{32} · 5^{16} · 5^2 · 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.
In [1]:
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
print(n,r,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:
In [86]:
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')
Out[86]:
[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.
In [89]:
for i, t1, t2 in zip(n, t_bruta, t_bin):
print('Ganancia para n = {0:>6}, {1:>10.2f}'.format(i, t1/t2))
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.