Un ejercicio de programación básica consiste en hallar los números naturales que se pueden descomponer en sumas de dos cubos de dos (o más) modos distintos. Se sabe, por ejemplo, que 1729 es el primer número que tiene esta propiedad, a saber:
1729 = 10^3 + 9^3 = 12^3 + 1^3
Tomando esto como idea base, podríamos preguntarnos por los siguientes números con tal característica. El siguiente programa permite encontrar estos números. De hecho, 1729 es el primer resultado:
#include <iostream>
using namespace std;
void prueba(long n);
int main()
{
long z;
for(z=0L;;z++)
prueba(z);
return 0;
}
void prueba(long n)
{
long x,y,x1,y1;
int maneras=0;
for(x=1;;x++)
{
for(y=1;y<=x;y++)
{
if(x*x*x+y*y*y==n)
{
if(maneras==0)
{
x1=x;y1=y;
maneras=1;
}
else
{
cout << "numero: " << n << endl;
cout << "x=" << x1 << " y=" << y1 << endl;
cout << "x=" << x << " y=" << y << endl;
return;
}
}
}
if(x*x*x>n)
return;
}
}Los siguientes son algunos valores que arroja la ejecución de “cubos0.cpp”:
numero: 1729 x=10 y=9 x=12 y=1 numero: 4104 x=15 y=9 x=16 y=2 numero: 13832 x=20 y=18 x=24 y=2 ...
Un aspecto que debe ser perfectamente comprendido es el camino de búsqueda de los pares (x,y) que elevados al cubo suman el valor buscado. Este recorrido se hace incrementando a “x” de uno en uno, y para cada uno de estos, se hace recorrer a “y” en el rango [1,x]. Esto garantiza que se recorran todas las posibles parejas (x,y) pero sin repetirse el caso invertido (y,x). En otras palabras, estamos evitando soluciones erróneas de tipo
N = a^3+b^3 = b^3 + a^3.
La ejecución genera números con rapidez, pero pronto se percibe una desaceleración conforme se busca entre números mayores.
A fin de enfocar nuestro problema con más claridad, restringiremos el ejercicio a encontrar “el siguiente número que cumple con la propiedad mencionada, pero que es mayor que cierto valor”. Por ejemplo, empezaremos buscando nuestros números a condición que sean mayores a 7'000,000. La elección de este valor ha sido totalmente arbitraria, y tan solo obedece a que mi computador se tarda unos pocos segundos en obtener la solución. El lector podría usar otro valor sin ningún inconveniente.
La siguiente versión incorpora los cambios mencionados:
/* cubos1.cpp: buscar el siguiente numero que cumple
la propiedad es mayor que 7 millones */
#include <iostream>
#include <cstdlib> // PRIMER CAMBIO
using namespace std;
void prueba(long n);
int main()
{
long z;
for(z=7000000L;;z++) // SEGUNDO CAMBIO
prueba(z);
return 0;
}
void prueba(long n)
{
.... TODO SIGUE IGUAL ....
cout << "x=" << x << " y=" << y << endl;
exit(0); // TERCER CAMBIO
}
}
}
if(x*x*x>n) return;
}
}Como se aprecia, hemos evitado repetir todo el programa, señalando los tres únicos cambios necesarios. A continuación se muestra la compilación y ejecución del programa con ayuda del comando “time” que permite obtener estadisticas temporales:
$ g++ -o cubos1 cubos1.cpp $ time ./cubos1 numero: 7081984 x=160 y=144 x=192 y=16 real 0m18.884s user 0m18.880s sys 0m0.000s $ time ./cubos1 numero: 7081984 x=160 y=144 x=192 y=16 real 0m18.885s user 0m18.890s sys 0m0.000s
Es sencillo verificar que la solución encontrada es correcta. Por otro lado, el tiempo de ejecución es de aproximadamente 18,8 segundos (en mi computador) y este valor es consistente a lo largo de varias ejecuciones. También apreciamos que el valor corresponde exclusivamente a tiempo “user”, lo cual es de esperarse pues no hay I/O involucrado; de otro lado, el tiempo “user” es prácticamente el mismo que el tiempo “real” debido a que mi computador no está haciendo absolutamente nada más en este momento.
La manera más rápida de mejorar los tiempos de ejecución es con ayuda de las opciones de optimización incluidas en el compilador. Estas opciones y su efectividad varían de computador en computador, de sistema operativo en sistema operativo y de compilador en compilador. En el caso de mi compilador g++ (gcc) emplearemos las opciones -O1, -O2 y -O3 para apreciar las mejoras:
$ g++ -O1 -o cubos1 cubos1.cpp $ time ./cubos1 $ time ./cubos1 numero: 7081984 x=160 y=144 x=192 y=16 real 0m9.723s $ g++ -O2 -o cubos1 cubos1.cpp $ time ./cubos1 numero: 7081984 x=160 y=144 x=192 y=16 real 0m7.742s $ g++ -O3 -o cubos1 cubos1.cpp $ time ./cubos1 numero: 7081984 x=160 y=144 x=192 y=16 real 0m7.740s
Por simplicidad hemos omitido los valores “user” y “sys”. En la optimización -O1 se aprecia una reducción de aproximadamente 49%, mientras que con -O2 y -O3 esta reducción alcanza aproximadamente 59%.
Los compiladores suelen tener muchas opciones de optimización, y el lector está invitado a probarlas. En particular “gcc” proporciona opciones de “arquitectura” que generan código específico para el procesador que se está empleando. Esto en ciertos casos mejora los tiempos dramáticamente, aunque en el presente programa no fue así por lo que no seguiré con este tema1. En lo que sigue, sólo asumiré compilaciones con -O3.
Ahora empezaremos a optimizar el algoritmo y veremos que los tiempos se pueden mejorar dramáticamente.
El programa anterior ciertamente fue bastante ingenuo: probaba con cada número para ver si se cumplía o no la propiedad. El problema es que con cada número se tiene que hacer una serie de cálculos que esencialmente son los mismos. Nos referimos a hallar el resultado de la operación
f=x^3+y^3
para cada pareja (x,y). Una manera distinta de proceder consiste en recorrer todas las parejas (x,y) por una única vez, y almacenar todos los resultados “f”, buscando repeticiones de estos últimos. Evidentemente, sólo tomaremos en consideración los “f” que sean mayores a nuestro valor de inicio de 7'000,000.
Sin embargo, este modo de proceder ocasiona un pequeño inconveniente. Si durante los cálculos hallamos dos valores coincidentes “f1=f2”, nada nos garantiza que éste sea exactamente el siguiente número mayor a 7'000,000 con la propiedad deseada. Ciertamente, será un número que cumple la propiedad, pero no necesariamente el inmediato consecutivo.
Una manera de salvar el asunto consiste en encontrar la menor de estas coincidencias “fa=fb”. Esta búsqueda se puede limitar con relativa facilidad, pues si se verifica
x^3 >= fu
(donde fu es la menor de las coincidencias halladas hasta el momento) entonces no hay manera de que los sucesivos “x^3+y\^3” obtengan un valor menor a “fu”.
En el siguiente listado, la memoria de valores hallados se guarda en el array “res[]”. La variable “st” (al inicio igual a -1) es un indice que apuntará al menor valor que cumple la propiedad hasta el momento. “res_offset” mantiene la última posición asignada de “res”. Nótese que hemos asumido que no requeriremos almacenar más de “MAX=1'000,000” valores para “res[]” (y “res_x[]”, “res_y[]” que son las parejas originantes.) Esto lo corregiremos más adelante mediante asignación dinámica de memoria. Finalmente, “st_x” y “st_y” contienen la “otra pareja” que genera la coincidencia “f1=f2”.
#include <iostream>
#include <cstdlib>
using namespace std;
#define MAX 1000000
long res[MAX];
long res_x[MAX];
long res_y[MAX];
int res_offset=0;
int main()
{
long n=1000000000L;
long x,y,st_x,st_y;
long r,j;
int st=-1; // Todavia no hay numero
for(x=1;;x++)
{
// Ya hay coincidencias, y estamos mas
// alla del rango
if(st!=-1 && x*x*x>=res[st])
{
cout << "numero: " << res[st] << endl;
cout << "x=" << res_x[st] << " y=" <<
res_y[st] << endl;
cout << "x=" << st_x << " y=" << st_y << endl;
cout << "offset=" << res_offset << endl;
exit(0);
}
for(y=1;y<=x;y++)
{
r=x*x*x+y*y*y;
if(r<=n)
continue;
// Buscar si ya esta en la bolsa
for(j=0;j<res_offset;j++)
{
if(res[j]==r)
{
// Este es el primer numero hallado con la
// propiedad, pero no necesariamente el <
if(st==-1)
{
st=j;
continue;
}
// Tenemos un nuevo numero que cumple
// la propiedad. Sera menor que el actual?
if(r<res[st])
{
st=j;
st_x=x;
st_y=y;
}
}
}
// No esta, asi que lo introducimos
res[res_offset]=r;
res_x[res_offset]=x;
res_y[res_offset]=y;
res_offset++;
}
}
}La ejecución no puede ser más alentadora:
$ time ./cubos2 numero: 7081984 x=160 y=144 x=192 y=16 offset=2380 real 0m0.012s user 0m0.010s sys 0m0.000s
Esto significa una reducción en tiempo de más de 99.9%! Este éxito rotundo no suele ser muy frecuente, pero es sorprendente la cantidad de lugares donde se podría obtener y se prefiere optar por una solución más simple: upgrade de hardware.
Es conveniente en este punto hacer un pequeño análisis acerca de el por qué de esta reducción. El lector que desea seguir con las optimizaciones puede saltarse esta sección sin perjuicio alguno.
La manera más usual de analizar y contrastar algoritmos corresponde al análisis del número de operaciones matemáticas involucradas
durante su ejecución. Considerando el primer algoritmo (cubos0.cpp, cubos1.cpp) tenemos que para cada número “n” se realiza un
conjunto de pruebas consistente en recorrer las parejas (x,y) hasta que x^3 > n, es decir, hasta que x > n^(1/3). Pero hasta
que “x” alcanza este valor, se ha hecho el cálculo de f = x^3 + y^3 en x(x+1)/2 oportunidades, debido a que “y” varía en [1,x] con
cada “x”. Si consideramos como operación fundamental dicho cálculo de “f”, entonces tenemos que para cierto “n” el cálculo se hace aproximadamente:
[n^(1/3)][n^(1/3)+1]/2 operaciones = K . n^(2/3) aprox.
donde K es una constante.
Sin embargo, hemos visto que esto se hace (ingenuamente) con todos los valores 1, 2, 3… hasta llegar a “n”, así que debemos sumar:
K [ 1^(2/3) + 2^(2/3) + 3^(2/3) + ... + n^(2/3) ]
Esta suma se puede aproximar mediante una integral, reduciendose a:
T total = K . n^(5/2)
Para el segundo algoritmo la cosa es diferente. A fin de llegar a “n” efectivamente tenemos que barrer con las parejas (x,y). Si x^3 >= n entonces con seguridad habremos alcanzado “n”, aunque el algoritmo normalmente se detendrá antes. Dado que el recorrido sigue el mismo camino que el anterior, también tendremos x(x+1)/2 operaciones, pero además tendremos por cada una de ellas una búsqueda en la “bolsa”. Esta búsqueda evidentemente tiene una duración proporcional al tamaño de la bolsa, el cual es directamente proporcional a la cantidad de operaciones “f” realizadas hasta el momento.
Por tanto, las operaciones “f” tomarán un tiempo directamente proporcional a:
x(x+1)/2 = K . n^(2/3) aprox.
Mientras que las búsquedas toman
M.1 + M.2 + M.3 +... + M[x(x+1)/2] = M . k^4 = M . n^(4/3) aprox. donde M es otra constante.
El tiempo total será por lo tanto:
T total = K . n^(2/3) + M. n^(4/3)
Es claro que conforme “n” crece, el segundo término sobrepasará largamente al primero. Por tanto, para “n” grande:
T total = M . n^(4/3)
Si se compara con el algoritmo anterior, es claro que el segundo es muy superior.
Antes de continuar, permítaseme cambiar el límite inferior del problema (7 millones) a un valor más elevado. El motivo es poder apreciar mejor las futuras optimizaciones, cosa que sería difícil de apreciar si nuestro tiempo ya ha sido reducido a 0,012 segundos. En ese sentido, elevaré el valor inicial de búsqueda a 1,000'000,000 (1E9.)
Siguiendo con las optimizaciones, un simple análisis demuestra que es un desperdicio iniciar nuestro recorrido de parejas (x,y) con x=1. En efecto, si “y” siempre está en [1,x], entonces:
y <=x ; luego: x^3 + y ^3 <= 2 x^3
Por tanto, x^3 >= (x^3 + y^3)/2 = n/2 ;
Finalmente, x >= (n/2)^(1/3).
Esta modificación se aprecia en el listado cubos3.cpp:
/* cubos3.cpp: x no se inicia en 1 */
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
#define MAX 1000000
long res[MAX];
long res_x[MAX];
long res_y[MAX];
int res_offset=0;
int main()
{
long n=1000000000L;
long x,y,st_x,st_y;
long r,j;
long start=(long)pow(n/2.0,1.0/3);
int st=-1; // Todavia no hay numero
cout << "inicio: " << start << endl;
for(x=start;;x++)
{
.... EL RESTO PERMANECE IGUAL ....Como se aprecia, hemos introducido un “n” inicial de 1E9 como anunciamos arriba, tanto en cubos2.cpp como en cubos3.cpp, con lo que obtenemos los siguientes resultados:
$ time ./cubos2 numero: 1003625000 x=900 y=650 x=975 y=425 offset=59980 real 0m4.990s $ time ./cubos3 inicio: 793 numero: 1003625000 x=900 y=650 x=975 y=425 offset=59980 real 0m4.973s
El valor de “n” incial ha elevado el tiempo a alrededor de 5 segundos. Sin embargo, no hay mayor diferencia en lo que respecta a nuestra última optimización; es decir, cubos2 y cubos3 son muy similares.
Hay dos explicaciones para ésto. En primer lugar, si bien hemos evitado probar para “x” en el intervalo [1,793>, los valores posteriores generan muchos más recorridos para “y” por lo que este ahorro se hace pequeño. En segundo lugar, la obtención de la raíz cúbica parece tomar un tiempo significativo. Esto lo comprobé introduciendo “x=793” en vez del cálculo respectivo.
Sin embargo, en general esta optimización no es muy trascendente.
En cubos3.cpp, y más precisamente en el recorrido de (x,y), la variación de “x” ha sido limitada tanto en el inicio (variable “start”) como en su final mediante la condición:
if(st != -1 && x*x*x >= res[st])Sin embargo nada hemos hecho con respecto al lazo interno. Veamos el inicio del lazo: Dado que tenemos un valor inicial de “n”, entonces:
x^3 + y^3 >= n (inicial), si se conoce “x” (pues estamos en el lazo interno), entonces tenemos:
y >= ( n - x^3)^(1/3)
Conforme crece “x”, este valor se podría hacer negativo, en cuyo caso forzaríamos “y=1”.
Esta optimización mejora ligeramente la performance, pero no drásticamente (sugiero que lo pruebe el lector.) Por otro lado, el extremo final del lazo puede cerrarse anticipadamente cuando ya tenemos una coincidencia (un valor de “f”) y estamos obteniendo valores superiores a éste (como sabemos, sólo nos interesan los posibles valores menores.) En tal sentido, podríamos plantear:
if(st!=-1 && r>res[st]) break;Este último cambio ciertamente es dramático:
$ time ./cubos4 inicio: 793 numero: 1003625000 x=900 y=650 x=975 y=425 offset=21343 real 0m0.596s
Es decir, hemos logrado una reducción de aproximadamente 88%! Otro aspecto que hemos dejado de lado hasta este momento corresponde al consumo real de memoria que se puede apreciar tras el texto “offset=” en la salida. La última optimización redujo este consumo casi a la tercera parte. A continuación los cambios necesarios en cubos3.cpp para convertirlo en cubos4.cpp. Obsérvese que sólo corresponden a “start2” y a la condición señalada arriba.
.... SIGUE IGUAL ....
int main()
{
long n=1000000000L;
long x,y,st_x,st_y;
long r,j;
long start=(long)pow(n/2.0,1.0/3);
int st=-1;
long start2; // OJO!
cout << "inicio: " << start << endl;
for(x=start;;x++) {
start2=(long)pow(n-x*x*x ,1.0/3.0);
if(start2<1)
start2=1;
if(st!=-1 && x*x*x>=res[st]) {
cout << "numero: " << res[st] << endl;
cout << "x=" << res_x[st] << " y=" <<
res_y[st] << endl;
cout << "x=" << st_x << " y=" << st_y << endl;
cout << "offset=" << res_offset << endl;
exit(0);
}
for(y=start2;y<=x;y++) {
r=x*x*x+y*y*y;
if(r<=n)
continue;
if(st!=-1 && r>res[st]) // OJO!
break;
for(j=0;j<res_offset;j++)
.... SIGUE IGUAL .....Hasta el momento hemos empleado un array de tamaño fijo para nuestro almacenamiento. Esto puede generar un desperdicio de memoria física y un desbordamiento si la “bolsa” se llena. La siguiente versión emplea asignación dinámica. Como se ve, no hay mayores diferencias, y el tiempo de ejecución es casi idéntico:
$ time ./cubos5 inicio: 793 numero: 1003625000 x=900 y=650 x=975 y=425 offset=21343 real 0m0.522s
Esta similitud de tiempos me sorprende, y de hecho se debe al excelente trabajo de los desarrolladores de la GNU Libc que estoy empleando. En general, es de esperarse que esta versión sea ligeramente menos veloz que la que emplea arrays.
/* cubos5.cpp: asignacion dinamica */
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
long *res=NULL;
long *res_x=NULL;
long *res_y=NULL;
int res_offset=0;
int main()
{
.... SIGUE IGUAL ....
res=(long *)realloc(res,(res_offset+1)*sizeof(long));
res_x=(long *)realloc(res_x,(res_offset+1)*sizeof(long));
res_y=(long *)realloc(res_y,(res_offset+1)*sizeof(long));
res[res_offset]=r;
res_x[res_offset]=x;
res_y[res_offset]=y;
res_offset++;
}
}
}En la sección de análisis de los algoritmos concluimos que el tiempo venía asociado principalmente a la duración de las búsquedas. En ese sentido se propone emplear árboles binarios como optimización final.
Como se sabe, los árboles binarios reducen el tiempo de búsqueda en una muestra de tamaño “n” a un valor proporcional a “log2(n)”, mientras que una búsqueda lineal como la que hemos venido efectuando toma un tiempo directamente proporcional a “n”2.
Lamentablemente, la implementación requiere de varios añadidos al programa. En particular, los arrays que implementaban la “bolsa” son reemplazados por un conjunto de estructuras de tipo “struct nodo” (rebautizadas como “Nodo”) que tienen la siguiente definición:
struct nodo
{
struct nodo *mayor;
struct nodo *menor;
long valor;
long x;
long y;
};Estas estructuras se enlazan entre sí mediante los punteros “mayor” y “menor”, conteniendo el valor de la operación “f”, y la pareja (x,y). En el esquema que mostraremos, cada nodo tiene (posiblemente) asociados dos nodos adicionales, correspondientes a aquellos que tienen un valor mayor o menor:
NODO 1
/ \
NODO NODO
MENOR MAYOR
Suponiendo que durante el recorrido (x,y) se hubiese obtenido los valores 80, 12, 34, 89, 13, 35 (en ese orden); entonces el árbol sería como sigue:
80
/ \
12 89
/ \ /\
34
/ \
13 35
/\ /\
La ventaja que proporciona este esquema radica en la facilidad de búsqueda, dado que con cada “salto” en promedio se elimina la mitad de la lista (en vez de un número a la vez en una búsqueda secuencial.)
El programa que se presenta a continuación implementa estas ideas. Nótese que se ha tratado de mantener cierta similitud con las variables de cubos5.cpp:
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
long *res=NULL;
long *res_x=NULL;
long *res_y=NULL;
int res_offset=0;
int main()
{
long n=1000000000L;
long x,y,st_x,st_y;
long r,j;
long start=(long)pow(n/2.0,1.0/3);
int st=-1; // Todavia no hay numero
long start2;
cout << "inicio: " << start << endl;
for(x=start;;x++)
{
start2=(long)pow(n-x*x*x ,1.0/3.0);
if(start2<1)
start2=1;
if(st!=-1 && x*x*x>=res[st])
{
cout << "numero: " << res[st] << endl;
cout << "x=" << res_x[st] << " y=" <<
res_y[st] << endl;
cout << "x=" << st_x << " y=" << st_y << endl;
cout << "offset=" << res_offset << endl;
exit(0);
}
for(y=start2;y<=x;y++)
{
r=x*x*x+y*y*y;
if(r<=n)
continue;
if(st!=-1 && r>res[st])
break;
for(j=0;j<res_offset;j++)
{
if(res[j]==r)
{
if(st==-1)
{
st=j;
continue;
}
if(r<res[st])
{
st=j;
st_x=x;
st_y=y;
}
}
}
res=(long *)realloc(res,(res_offset+1)*sizeof(long));
res_x=(long *)realloc(res_x,(res_offset+1)*sizeof(long));
res_y=(long *)realloc(res_y,(res_offset+1)*sizeof(long));
res[res_offset]=r;
res_x[res_offset]=x;
res_y[res_offset]=y;
res_offset++;
}
}
}Y a continuación, los tiempos:
$ time ./cubos6 inicio: 793 numero: 1003625000 x=900 y=650 x=975 y=425 real 0m0.031s $ time ./cubos6 inicio: 793 numero: 1003625000 x=900 y=650 x=975 y=425 real 0m0.051s
Es decir, nuevamente tenemos una reducción de más de 90 %!Lamentablemente, los tiempos nuevamente son muy pequeños y por tanto no son muy estables como se puede apreciar. Podríamos elevar nuevamente nuestro límite inferior y continuar intentando, pero tenemos que tener cuidado: los enteros tipo “long” típicamente tienen un límite de aproximadamente 2E9 (o 4E9 si no tienen signo) así que deberíamos analizar si disponemos de enteros de 64 bits3.
-
El compilador puede hacer mucho por nosotros sin pedir mucho a cambio. Hay que utilizarlo!
-
El análisis cuidadoso de los algoritmos puede proporcionar mejoras dramáticas o inverosímiles en el tiempo de ejecución de los programas. En muchos casos esto ocurrirá sin mucho esfuerzo mientras que en otros se requiere una buena dosis de complejidad
-
Los programas con malos algoritmos suelen desencadenar un innecesaria adquisición de hardware. Por ejemplo, en los análisis de “capacity planning”, se suele considerar la carga del CPU como factor decisivo. Sin embargo, los algoritmos mal implementados del inicio de este artículo, si bien son muy lentos, de todas maneras cargarán al máximo al CPU por lo que exteriormente parecerán faltar recursos de hardware.
Un agradecimiento especial a Carlos Salvatierra que revisó este texto e hizo ciertas sugerencias que espero aparezcan en una próxima versión.