Factoriones

Hace tiempo que no entro a resolver problemas en Project Euler. Si te gusta resolver acertijos matemáticos, que seguido involucren habilidades de programación, es bastante recomendable. El punto es que recordé un problema que implícitamente pide encontrar factoriones.Los pide implícitamente pues sólo existen cuatro factoriones en base 10. Un factorion es un número que es igual a la suma de factoriales de sus dígitos, e.g. 145=1!+4!+5!

Para hallarlos simplemente hay que programar directamente el algoritmo para discriminar si un número es o no factorion. El verdadero problema es que el programa no pase días buscando factoriones. Para lograr esto, el procedimiento estándar es notar que k(9!) tendrá menos de k dígitos siempre que k>7, por tanto sólo hay que buscar factoriones hasta 7(9!)=2540160

Ahora bien, esto aún sigue sin ser muy eficiente. Esto es simplemente un puzzle de programación, como dije sólo existen cuatro factoriones, se conocen y punto.

El factorion más grande es 40585, así que probablemente exista una forma de recortar el límite superior a partir de un argumento lógico. Mi lenguaje ya predilecto es Python, y no sé si sea precisamente si es él el que hace lento el programa.

Acá comparto el código:

# 145 es un numero curioso, ya que 1! + 4! + 5! = 1 + 24 + 120 = 145.
# Encuentra la suma de todos los numeros iguales a la suma
# del factorial de sus digitos.
# Nota: como 1! = 1 y 2! = 2 no son sumas, no se incluyen.

import time
t0 = time.time()

def factorial(n):
    if n > 1: return (n*factorial(n-1))
    else: return 1

for i in range(1, 7*factorial(9)): # LIMITE SUPERIOR ÓPTIMO~1*(10^5)
    n, s = str(i), 0
    for j in range(len(n)):
        s += factorial(int(n[j]))
    if s == i: print(i)
print('Tiempo de ejecucion: ',time.time()-t0,'segundos.')

La espiral de Ulam y project euler 58

Llevo dos días sin poder resolver el problema 51 de Project Euler, así que decidí intentar algún otro. Encontré el que muestro a continuación, que me pareció divertido y sencillo, el único detalle es que mi solución no es tan eficiente como pensé.

El motivo de esta entrada es desafiarte a resolverlo mejorando el análisis que yo hice (es decir, el tiempo que tardó mi programa si es que programarás en Python); no importa cómo lo resuelvas, lo importante es que la solución sea elegante, sencilla y rápida.

A continuación muestro el problema, mi resultado y tanto el análisis como el programa que hice para obtenerlo. Mi recomendación es que primero lo resuelvas a tu manera. Si te interesa resolver más problemas entra a project euler y lleva tu propio récord.

Comenzando con 1 y rotando en sentido contrario a las agujas del reloj de la siguiente manera, se forma una espiral cuadrada con lados de longitud 7.

37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49

Es interesante que los cuadrados impares se sitúan en la diagonal inferior derecha, pero aún más interesante es que 8 de los 13 números de las dos diagonales son primos, es decir, una proporción de 8/13 ≈ 62%.

Si enrollamos una nueva capa en esta espiral, se formará una espiral cuadrada con lados de longitud 9. Si este proceso continúa, ¿cuál es la longitud del lado de la espiral cuadrada para la cual la relación de primos en diagonales es inferior al 10%?

Mi respuesta:

Side length of the square spiral: 26241
Time elapsed: 26.3150000572 seconds

Análisis y código fuente:
Comenzando desde 1, el número de la $i$-ésima esquina superior derecha de la espiral cuadrada, $k_i$ está dado por
$$k_i=k_{i-1}+8(i-1)+2\hspace{0.25in}\text{con}\hspace{0.25in}k_0=1$$ Para las 3 esquinas restantes (que etiqueto con 1, 2 y 3), en cada vuelta
\begin{align*}k_{(i,1)}&=k_i+2i\\
k_{(i,2)}&=k_i+4i\\
k_{(i,3)}&=k_i+6i\end{align*} Así entonces, generamos
\begin{align*}k_1&=k_0+2=1+2=3\\
k_{(1,1)}&=k_1+2(1)=3+2=5\\
k_{(1,2)}&=k_1+4(1)=3+4=7\\
k_{(1,3)}&=k_1+6(1)=3+6=9\end{align*} Falta revisar si los números son primos y contar. En el problema se menciona que en la diagonal inferior izquierda todos son cuadrados impares, así que éstos se pueden ignorar.

Mi código en Python se ve así:

from __future__ import division
import time

t0 = time.time()

def isprime(x, i=2):
    while i <= (x/i):
        if x%i == 0: return False
        i += 1
    return True

prime, total = 0, 1
k, n = 3, 2
while prime/total >= 0.1 or prime == 0:
    for i in range(3):
        if isprime(k+i*n) == True:
            prime += 1
    k = k + 4*n + 2
    total += 4
    n += 2

print 'Side length of the square spiral:',n-1
print 'Time elapsed:',time.time()-t0,'seconds'

Apenas me he enterado que la espiral utilizada se llama Espiral de Ulam, para el que esté interesado en conocer más. La espiral en general, por ejemplo, muestra la siguiente distribución de números primos:

Números de Lychrel, Algoritmo196 y Project Euler 55

Uno de los problemas abiertos de Teoría de Números es encontrar números de Lychrel. Estos números hipotéticos desafían el llamado algoritmo196, -una secuencia que utiliza el viejo truco de sumar y agregar para obtener números palíndromos o capicúas. El algoritmo196 toma cualquier número natural, invierte sus dígitos y le suma el original; a nuestro interés, si la suma no es capicúa repite el proceso a partir de dicha suma, por ejemplo

338 + 833 = 1171
1171 + 1711 = 2882

Los números de Lychrel son entonces, los que no producen capicúas al aplicar el algoritmo196. De nuevo, esto es un problema abierto, y precisamente 196 es el primer número de Lychrel que se ha conjeturado. Consulta La conjetura del 196 en Gaussianos. En seguida muestro una parte del problema 55 de Project Euler y mi solución, que es vulgar y silvestre lógica potenciada con el lenguaje Python, en realidad no hice un análisis más allá.

[…] Aunque nadie lo ha demostrado aún, se cree que algunos números, como el 196, nunca producen un capicúa. Un número que nunca forma un capicúa mediante el proceso de sumarse a sí mismo invertido se denomina un número de Lychrel. Debido a la naturaleza teórica de estos números, y para el propósito de este problema, vamos a suponer que un número es Lychrel hasta que se demuestre lo contrario. Además, se te da la información de que cada número inferior a diez mil, (i) será un capicúa en menos de 50 iteraciones, o, (ii) nadie, con toda la potencia de cálculo que existe, ha logrado por ahora asociarlo a un capicúa.

De hecho, 10677 es el primer número que requiere más de cincuenta iteraciones antes de producir el capicúa: 4668731596684224866951378664 (53 iteraciones, 28 dígitos).

Sorprendentemente, existen números capicúa que son a su vez números de Lychrel; el primer ejemplo es 4994.

¿Cuántos números de Lychrel hay inferiores a diez mil?

Solución en Python:

import time

t0 = time.time()
L, M = 0, 10000
for i in range(1,M):
    s = str(i+int(str(i)[::-1]))
    lychrel = True
    j = 1
    while j <= 50:
        if s == s[::-1]:
            lychrel = False
            break
        j += 1
        s = str(int(s)+int(s[::-1]))
    if lychrel == True:
        L += 1
print 'Hay',L,'números Lychrel inferiores a',M
print 'Tiempo de cómputo:',time.time()-t0,'segundos'

Hay 249 números Lychrel inferiores a 10000
Tiempo de cómputo: 0.345999956131 segundos

Como he dicho, la solución no son en realidad números de Lychrel, y asimismo el encontrarlos o negarlos quizás dependa más de una astuta demostración que de potentísimos ordenadores.

La millonésima permutación lexicográfica

De nuevo Project Euler, tal vez con el programa más difícil que me he encontrado hasta ahora. Desde que empecé a programar, la recursión (recursividad o recurrencia) me ha parecido de las formas más difíciles y sin embargo más elegantes de programar (programación funcional). Para hacer permutaciones me parece inevitable utilizar recursión. Cito el problema 24 de Project Euler:

Una permutación es una disposición ordenada de objetos. Por ejemplo, 3124 es una  posible permutación de los dígitos 1, 2, 3, 4. Si todas las permutaciones se listan numérica o alfabéticamente, a esta disposición le llamamos orden lexicográfico. Las permutaciones lexicográficas de 0, 1, 2 son: 012   021   102   120   201   210. ¿Cuál es la millonésima permutación lexicográfica de los dígitos 0, 1, 2, 3, 4, 5, 6, 7, 8, 9?

El problema parece sencillo; después de programar los problemas 19 y 28 de forma casi trivial, me topé con este problema que simplemente no podía lograr. Prácticamente siempre empiezo un programa de la forma más sencilla, lo pruebo así y lo llevo después a la forma complicada. El ejemplo que da Project Euler puede resultar engañoso, y como yo, cualquiera puede caer en un pozo sin fondo.

Sabemos que el número de permutaciones posibles de un conjunto de $n$ objetos es igual a $n!=n\cdot(n-1)\cdot3\cdot2\cdot1$. El problema de comenzar con casos sencillos, es que, por ejemplo para $n=3$, un simple bucle anidado de 2 niveles es suficiente para repetir 3 veces un bucle que se repite 2 veces, obteniendo las $3!=6$ permutaciones buscadas, en cambio, si tomamos el caso $n=4$, necesitaríamos un bucle de 3 niveles que hiciera 4 veces un bucle que repite 3 veces otro bucle que se repite 2 veces, obteniendo $4!=24$ permutaciones... y así sucesivamente para cualquier $n$ usando este razonamiento, lo cual no es práctico ni mucho menos, recordando que el problema que queremos resolver nos pide $10!$ permutaciones.

Bueno, pues ya que sabemos que se obtendrán $10!$ permutaciones, debemos recordar cómo utilizar la técnica de recursividad. Regularmente cuando se introduce este concepto se utiliza el ejemplo del factorial. En Python se ve algo así:

def f(n):
        if n > 1 :
            return (n*f(n-1))
        else:
            return 1

Bien pues se necesita usar este tipo de razonamiento pero buscando formar permutaciones. Una buena forma de lograrlo es intentarlo primero con un conjunto de cuatro objetos, lo más sencillo a empezar con 0123. Un razonamiento que funciona es intercambiar un par externo e ir avanzando hacia el siguiente término, como se observa el patrón que a continuación muestro. Las primeras permutaciones se verían así:

0123
0132
0213
0231
0312
0321
1023
1032
1203
1230
1302
1320
2013
$\vdots$

Notamos fácilmente que para la primer posición de izquierda a derecha, el 0 se mantiene fijo y los términos adyacentes permutan $6=3!$ veces, y de igual modo, si nos fijamos en la posición siguiente, en todas las permutaciones que comienzan con 0, cada objeto tiene $2=2!$ permutaciones, y esto fácilmente lo podemos predecir para cualquier combinación siguiente. Si observáramos una posición aún anterior, veríamos que 0, 1, 2, y 3 forman exactamente $4!$ permutaciones, que es lo que buscamos. También es fácil notar que no necesitamos ordenar permutaciones para formar un orden lexicográfico si comenzamos con una cadena ordenada, i.e. 0123456789.

Ahora el trabajo arduo (al menos para mí lo fue) es describir esto en un lenguaje estricto y obtener lo que se busca. Recomiendo primero encontrar el patrón que describen las permutaciones y luego intentar codificarlo. Casi siempre que resuelvo un nuevo problema de Project Euler parece que ha sido el más retador, sin embargo éste se ha llevado definitivamente el título, además de que la recursividad como dije, no me parece algo tan sencillo.

Python (y supongo que otros lenguajes más) tienen funciones predefinidas que dan permutaciones, sin embargo ¿qué de divertido habría en sólo buscar la millonésima de éstas, o si acaso tener que ordenar las permutaciones?.

Este tipo de problemas es digno de discusión con profesores o colegas/compañeros si es que te estancas; en mi caso pasé 2 días y 2 noches pensando y lidiando con la recursividad, aunque el patrón sea muy sencillo de descifrar.

Mi código es un tanto feo aunque no es tan largo y tuve algunos problemas con las posiciones del arreglo que contiene cada permutación. Dentro del programa cree un archivo en el cual escribí las permutaciones y al final la millonésima permutación lexicográfica, lo que como vimos arriba sería innecesario considerando que $10!>1,000,000$ pudiéndonos detener entonces antes de que el programa se ejecute totalmente.

Aquí dejo las permutaciones del conjunto 012345678 para que te puedas guiar con más datos acerca del patrón formado y el algoritmo que debes codificar.

HAPPY CODING!

:)

Narcisistas, ¿con Python o C++?

Bien, pues he comenzado a descubrir ventajas de Python o que tengo debilidades en C++. Los siguientes programas encuentran números narcisistas, esto es, un número de n dígitos que es la suma de la enésima potencia de sus dígitos (por ejemplo ${407=4^3+0^3+7^3}$ es un 3-narcisista).

En los programas P es el número de dígitos de los narcisistas y ambos programas tardan un tiempo razonable (máx. 1 min) con ${P\leq{7}}$.

C++:

long int potencia (int,int);
long int look (int);

#define P 7

int main(void)
{
       cout << "Numeros " << P << "-narcisistas:" << endl;
       for(int i=potencia(10,P-1); i<potencia(10,P); i++)
       {
             if(look(i)==i){cout << i << endl;}
       }
       cout << endl << "Clic para finalizar" << endl;
       getch();
       return 0;
}

long int potencia(int a, int b)
{
       if(b>0) return a*potencia(a, b-1);
       else return 1;
}

long int look (int i)
{
       if (i>9) { return look(i/10)+look(i%10);}
       else { return potencia(i%10,P);}
}

Implementé el programa de este modo porque no pude hacerlo de modo más sencillo, como lo hice con Python. En Python el programa simplemente aumenta número por número y digamos que “desmiembra” dígito por dígito, a partir de un entero. Aquí simplemente no pude hacerlo; intenté convertir todo el entero a cadena y convertir caracter por caracter a entero. No funcionó. Así limpiara el arreglo original, usara punteros o memoria dinámica. Y bueno, el código que muestro sólo tiene de más la función "potencia", pero me quise ahorrar el problema de meter una variable double y hacer de más malabares para no tener problemas. También es notable el que haya usado recursión en la función "look" porque aún me cuesta bastante usar esta forma de programación. Aunque el resultado me parece elegante, fue bastante difícil lograrlo; el algoritmo no es inmediato y la recursividad, aunque ahorra fácilmente 5 o más líneas de declaraciones e intercambios, me resulta muy complicada para lo que busco que haga el programa.

En cambio véase con qué facilidad se codifica en Python:

t,P = 0,7

print "Numeros  ",P,"- narcisistas:"

for i in range(10**(P-1), 10**P):
        for j in str(i):
                t += int(j)**P
        if t == i:
                print t
        t=0

Ya todo lo he descrito y el programa habla por sí solo con su sencillez. Lo único que me pareció inconveniente, y aún no me lo explico del todo, es que tardó mucho más la ejecución que con C++.

En fin, ¡a seguir descubriendo Python y a afinarnos en C++ se ha dicho!

El periodo más grande en una parte fraccionaria

Me tomó como una semana resolver este problema, y aun así lo he hecho de una forma bastante grotesca pues no he tenido tiempo para estudiarlo o para entender formas más elegantes de resolverlo. He estado bastante ocupado con las actividades curriculares, por lo que no he trabajado más en ello, aunque al menos he aprendido un poco más el uso de Matlab e intentaré aprender también SciPy (de Python) para resolver problemas de análisis numérico.

El problema es el 26 de Project Euler y se puede resumir a: Encuentra el valor de d<1000 para el cual 1/d contiene el periodo más grande en su parte fraccionaria. Mi sorpresa ha sido que la solución se obtiene elegantemente con teoría de números y álgebra abstracta al resolver el algoritmo discreto ${10^k\equiv1\mod{n}}$. Sin embargo son herramientas de las cuales aún no tengo pleno conocimiento; como sea, quien esté interesado puede consultar esta liga.

Estuve bastante ansioso por resolver el problema, y finalmente lo hice, aunque fuera de forma rudimentaria, pero encontrando algunas cosillas interesantes.

Lo difícil fue determinar cuándo la parte fraccionaria decimal era periódica. Digamos que tenemos 0.3333… bueno, resulta muy sencillo, ahora digamos que tenemos 0.323232… aún es sencillo, ahora bien, para el número
0.61029367467310296743597463610293674673102967435974636102936743…
bueno, tal vez ya no sea tan sencillo.

Debemos notar que para 1/d con parte fraccionaria periódica, el periodo máximo que se puede obtener es d-1; ésta es una propiedad demostrada y fácilmente se puede intuir. Bien, de esto sabemos que un algoritmo se volverá más eficiente si comenzamos a buscar a partir del d más grande posible, con ello al momento de encontrar un número con periodo fraccionario decimal d-1, nos detenemos y damos el resultado.

Bien, mi solución guarda el doble del máximo de dígitos que pueda tener cada periodo fraccionario en un arreglo; si tenemos 1/7=0.142857..., el arreglo contendrá 142857142857, y seguido podemos dividir éste arreglo en dos partes iguales y revisamos si éstas son idénticas (d-1 será par). Esto no garantiza que el periodo sea máximo, pues podrían haber, digamos, "subperiodos", desde simples como 1/9=0.111... hasta por ejemplo, 1/977 que tiene periodo 166, pero 976/166=6, es decir, ¡en cada mitad del arreglo habrán 6 periodos de 166 que serán idénticos!, lo que nos podría llevar a pensar incorrectamente que el periodo es 976 (digo incorrectamente pues aunque este último sí es un periodo, no es el que nos interesa).

De lo anterior se sigue inmediatamente un criterio para discernir cuándo el periodo buscado es d-1; i.e. cuando las mitades de una sola mitad del arreglo sean distintas, como en 0.142857 podemos ver que 142$\neq$857, por ejemplo. Tendríamos que expandir o reconsiderar este criterio si es que pudieran haber cantidades impares de subperiodos (¿puedes decir si esto ocurre o no, y por qué?), además de considerar los más simples (mi programa no los considera todos).

El siguiente es un código en Python que resuelve en problema. Yendo un poco más allá del problema de project euler, para d<10^5, funciona en 14.2419729233 segundos bajo Linux (Ubuntu) dando el resultado d=99989.

from decimal import *
import time
t0 = time.time()

for i in range((10**5)-1, 1, -2):
    getcontext().prec = 2*(i-1)
    cad = str(Decimal(1)/Decimal(i))
    if len(cad) > i+1:
        a, b = [], []
        for x in range(i+1, 2*i):
            a.append(cad[x-i+1])
            b.append(cad[x])
        c = 0
        for x in range(i-2):
            if a[x] == a[x+1]:
                c += 1
        if c >= i-2:
            continue
        else:
            if a != b:
                continue
            else:
                c, d = [], []
                for n in range(0,len(a)/2):
                    c.append(a[n])
                    d.append(a[(len(a)/2)+n])
                if c == d:
                    continue
                else:
                    print "Respuesta:",i
                    break
print "Duracion:",time.time()-t0, "segundos"