Euler #37 - "Truncatable Primes" - los algoritmos todavía importan

Retomando problemas del proyecto Euler, este en principio es bastante fácil.

https://projecteuler.net/problem=37

Es relativamente fácil, por ejemplo, con Python, convertir el número a string, ir quitando caracteres por la derecha y comprobando cada vez si el número que queda es primo. Luego, lo mismo por la izquierda. 

Mi primera implementación no tuvo en cuenta que no se considera "truncatable" si terminas con un 1. Lo mismo deberían decirlo en el enunciado, lo mismo es falta de cultura matemática por mi parte.

La función is_truncable, que se muestra a continuación, utiliza una función is_prime para trabajar. Va probando casos para probar que los números que van quedando son primos o no. Si después de analizar por la derecha y la izquierda todos han resultado ser primos, asume que es "truncable".

def is_truncable(number):

str_number = str(number)

# 2, 3, 5, 7 no son truncables

if number < 11:

return False

# si empieza o acaba en 1, no es truncable

elif str_number[0] == '1' or str_number[-1] == '1':

return False

else:

# Por la derecha

str_number = str(number)

while len(str_number) >= 1:

if not is_prime(int(str_number)):

return False

str_number = str_number[:-1]

# Por la izquierda

str_number = str(number)

while len(str_number) >= 1:

if not is_prime(int(str_number)):

return False

str_number = str_number[1:]

return True

Antes de eso, hemos hecho una pequeña función para generar una lista de números primos:

def next_prime(start):

next = start + 1

while not is_prime(next):

next += 1

return next

Y una función main() que monta todo el proceso:

def main():

truncables = []

next = 10;

# nos piden los 11 primeros truncables

while len(truncables) < 11:

next = next_prime(next, is_prime)

if (is_truncable(next, is_prime)):

truncables.append(next)

print("Truncables:", truncables)

print("Sum. truncables =", sum(truncables))

1ª versión malísima

En la primera implementación de is_prime() que hice, sin ninguna optimización, propia de un estudiante de ESO, simplemente iba dividiendo el número en cuestión por todos los números más pequeños que él mismo, si encontraba algún divisor, entonces no era primo.

def is_prime(number):

if number in (2,3):

return True

if number % 2 == 0:

return False

number_test = number - 1;

while number_test > 1:

resto = number % number_test

if resto == 0:

return False

number_test -= 1;

return True

Los problemas que plantean en proyecto Euler están muy bien pensados. Usando esta implementación ingenua, al llegar al 10º número truncable de la lista, se quedaba el problema clavado. 

4ª versión, usando librerías externas

Nos vamos al extremo contrario, usamos alguna librería existente (que seguramente esté bien optimizada y programada) para saber si un número es primo o no. El módulo sympy define una función isprime().

import sympy

def is_prime(number):

return sympy.isprime(number)

Con esta implementación, el problema se resuelve en un tiempo más que razonable, pero no vale, estamos haciendo trampa. 

2ª versión: mejor, pero sigue siendo mala

Esta implementación trata de ser un poco mejor: solo comprueba divisores hasta la mitad del número y hace una comprobación previa para los pares:

def is_prime(number):

if number in (2,3):

return True

if number % 2 == 0:

return False

for n in range(2, int(number/2)):

if number % n == 0:

return False

return True

Tampoco consigue resolver el problema en un tiempo razonable.

3ª versión: razonablemente buena

La siguente implementación de is_prime(), sin ser un prodigio, ya es bastante más óptima:

• Hace una primera prueba contra los divisores más pequeños que ya sabemos que son primos.
• Si es par, no es primo.
• Para números mayores que 5, solo hace las comprobaciones hasta la raíz cuadrada + 1 del número + 1 en cuestión, y sólo comprueba con divisores impares.

def is_prime(number):

if number <= 3 or number == 5:

return True

if number % 2 == 0 or number % 3 == 0 or number % 5 == 0:

return False;

for n in range(5, int(math.sqrt(number)) + 1, 2):

if number % n == 0:

return False

return True

Con esta implementación un poco mejorada, ya tenemos tiempos de resolución del problema razonables y del mismo orden de magnitud que usando la librería sympy.

Para poder hacer todas las diferentes pruebas, he metido todo el código común en un módulo y desde el script principal, pasamos la función is_prime correspondiente.

# Fichero euler37lib.py

def next_prime(start, is_prime):

next = start + 1

while not is_prime(next):

next += 1

return next

def is_truncable(number, is_prime):

str_number = str(number)

# 2, 3, 5, 7 no son truncables

if number < 11:

return False

# Si empieza o acaba en 1, no es truncable

elif str_number[0] == '1' or str_number[-1] == '1':

return False

else:

# Por la derecha

str_number = str(number)

while len(str_number) >= 1:

if not is_prime(int(str_number)):

return False

str_number = str_number[:-1]

# Por la izquierda

str_number = str(number)

while len(str_number) >= 1:

if not is_prime(int(str_number)):

return False

str_number = str_number[1:]

return True

def main(is_prime):

truncables = []

next = 10;

# Nos piden los 11 primeros truncables

while len(truncables) < 11:

next = next_prime(next, is_prime)

if (is_truncable(next, is_prime)):

truncables.append(next)

print("Truncables:", truncables)

print("Sum. truncables =", sum(truncables))

Y aquí están las cuatro diferentes implementaciones. Primero, haciendo trampa.

# Fichero 037-sympy.py

import sympy

import euler37lib

def is_prime(number):

return sympy.isprime(number)

if __name__ == '__main__':

euler37lib.main(is_prime)

Refinando un poco el algoritmo:

# Fichero 037-refinado.py

import math

import euler37lib

def is_prime(number):

if number <= 3 or number == 5:

return True

if number % 2 == 0 or number % 3 == 0 or number % 5 == 0:

return False;

for n in range(5, int(math.sqrt(number)) + 1, 2):

if number % n == 0:

return False

return True

if __name__ == '__main__':

euler37lib.main(is_prime)

Implementación ligeramente mejor a la ingenua, pero aún así impracticable:

# Fichero 037-pasable.py

import euler37lib

def is_prime(number):

if number in (2,3):

return True

if number % 2 == 0:

return False

for n in range(2, int(number/2)):

if number % n == 0:

return False

return True

if __name__ == '__main__':

euler37lib.main(is_prime)

Implementación ingenua:

# Fichero 037-ingenuo.py

import euler37lib

def is_prime(number):

if number in (2,3):

return True

if number % 2 == 0:

return False

number_test = number - 1;

while number_test > 1:

resto = number % number_test

if resto == 0:

return False

number_test -= 1;

return True

if __name__ == '__main__':

euler37lib.main(is_prime)

Comparativa de rendimiento

Usando librería: 

Implementación "refinada"

Nuestra implementación "refinada" tiene tiempos comparables a la implementación usando el módulo sympy.

Implementación mala:

Esta otra pasa de apenas un segundo de las anteriores a 9 minutos. Es 540 veces más lenta.

Implementación inservible:

5.100 segundos. Las dos primeras implementaciones tardaban como 1 segundo. Un algoritmo 5.000 veces más lento.

Conclusiones

1. Los algoritmos sí que importan. Si este programa, por ejemplo, fuese código para producción, ninguna de las dos implementaciones sencillas sería aceptable.
1. Con matices: si en el problema se nos garanzizase que, por ejemplo, siempre tendremos datos de entrada con números por debajo de 100, no hay ninguna mejora apreciable en el rendimiento por usar una u otra versión de is_prime(). 
   
   
2. Si existe una librería o módulo que hace lo que necesitas y es de pago, cómpralo o convence a tu jefe/cliente para que lo compre [*]:
1. Funcionará mejor. 
2. Tendrá soporte.
3. A la larga será más económico.

[*] Es una generalización un poco burda, tampoco vamos a estar metiendo librerías o código de terceros para todo, pero es un poco la idea: los "vendor" suelen hacer las cosas mejor que lo que tú puedas hacer.

Euler #35

https://projecteuler.net/problem=35

La (pequeña) dificultad de este problema consiste en calcular todas las posibles rotaciones de un número. La sintaxis de Python hace muy cómoda esta tarea.

import sympy

def check_rotations(number):
 all_primes = True
 number_str = str(number)
 number_str = number_str[1:] + number_str[:1]
 while int(number_str) <> number:
  if not sympy.isprime(int(number_str)):
   all_primes = False
   break
  number_str = number_str[1:] + number_str[:1]
 return all_primes

circular_primes = []
for n in range(1, 1000000):
 if sympy.isprime(n) and check_rotations(n):
  circular_primes.append(n)

print circular_primes
print len(circular_primes)

Puntos a destacar:

• sympy es una librería de Python especializada en matemáticas con símbolos (polinomios, resolución de ecuaciones, matrices,...) En nuestro script simplemente utilizamos la función isprime() para comprobar si el número es primo. Podríamos haber reutilizado alguna de las funciones que ya escribimos en problemas anteriores.
• La expresión "cadena"[1:] devuelve todos los caracteres exceptuando el primero: "adena"
• La expresión "cadena"[:1] devuelve el primer caracter, podríamos haber escrito también "cadena"[0]
• La sintaxis de subíndices de Python garantiza siempre que s == s[:i] + s[i:] 

Una posible optimización es recorrer solo los números impares, fijando 2 como único par a incluir:

circular_primes = [2]
for n in range(3, 1000000, 2): # El tercer parámetro de range es "step"

Euler #34

Tras un largo parón, otro problema del proyecto Euler, de nuevo lo solucionamos con Python. Tratando los números como cadenas de texto y convirtiendo los dígitos de nuevo a enteros el problema es muy sencillo de resolver.

https://projecteuler.net/problem=34

def fact(n):
 if n <= 1:
  return 1
 else:
  return n * fact(n - 1)

candidates = []
n = 3
while True:
# print "Analizyng", n
 n_str = str(n)
 partial_sum = 0
 for digit in n_str:
  partial_sum += fact(int(digit))
 if partial_sum == n:
  candidates.append(n)
 if n > 1000000:
  break
 n += 1

#print candidates
print sum(candidates)

Euler #33

https://projecteuler.net/problem=33

De nuevo, la facilidad que ofrece Python para pasar de unos tipos a otros y los módulos estándar que trae hacen que este problema sea trivial.

Se puede compactar más la solución, pero así se ve claro cómo vamos filtrando poco a poco los números candidatos hasta llegar a la solución. También ayuda mucho el enunciado del problema, estableciendo cuántos resultados son esperables.

#!/usr/bin/env python
from fractions import Fraction
from operator import mul

fracts = []
for d in range(10, 100):
 for n in range(10, 100):
  if n / d < 1 and n % 10 != 0 and d % 10 != 0:
   num = str(n)
   den = str(d)
   for x in range(1, 10):
    xstr = str(x)
    if (num[0] == xstr and den[1] == xstr) or (num[1] == xstr and den[0] == xstr):
     new_n = int(num.replace(xstr, '', 1))
     new_d = int(den.replace(xstr, '', 1))
     if (new_n * d == new_d * n):
      fracts.append(Fraction(new_n, new_d))

print reduce(mul, fracts, 1)

La función "reduce" es muy útil para operar con listas en bloque. Es un atajo para no tener que escribir bucles:

tot = 1
for f in fracts:
    tot *= f
print tot

Así mismo, la clase Fraction es lo suficientemente inteligente para lidiar con enteros, operar con fracciones, simplificarlas,...

Euler #32

Llevaba bastante tiempo parado con la resolución de problemas del proyecto Euler.
Vamos con el número 32. La resolución con Python es muy fácil convirtiendo enteros en "strings":

#!/usr/bin/env python

numbers = []

for a in range(9999):
  for b in range(9999):
    p = a*b
    asString = str(a) + str(b) + str(p)
    if len(asString) > 9:
      break
    if len(asString) == 9 and \
     '1' in asString and \
     '2' in asString and \
     '3' in asString and \
     '4' in asString and \
     '5' in asString and \
     '6' in asString and \
     '7' in asString and \
     '8' in asString and \
     '9' in asString:
     if not p in numbers:
       numbers.append(p)
       print a, "*", b, " = ", p, " => ", asString

print sum(numbers)

Euler #30

https://projecteuler.net/problem=30 => "Digit fifth powers"

Un poco de culturilla matemática, se trata de "números narcisistas"

def checkSumDigits(number, power):
 _sum = 0
 for digit in str(number):
  _sum += (int(digit))**power
 return _sum == number

final_sum = 0
for a in range(2, 10**6):
 if checkSumDigits(a, 5):
  print a
  final_sum += a
print final_sum

Euler #29

https://projecteuler.net/problem=29 => Distinct Powers

numbers = []
for a in range(2, 101):
 for b in range(2, 101):
  n = a**b
  if not n in numbers:
   numbers.append(n)

numbers.sort() 
print len(numbers)

Euler #28

https://projecteuler.net/problem=28 => Number spiral diagonals

Los problemas en "2D" son siempre entretenidos. El problema es bastante sencillo, pero hay que tener cuidado y llevar un control de en qué dirección vamos y a cuándo "hay que torcer".

import sys

if len(sys.argv) != 2:
 print "Usage: 028.py <SIZE>"
 sys.exit()

table = []
size = int(sys.argv[1])
if size % 2 == 0:
 print "Only odd numbers"
 sys.exit()

for _ in range(0, size):
 row = []
 for _ in range(0, size):
  row.append(0)
 table.append(row)

directions = ((0,1), (1,0),(0,-1),(-1,0))

def get_next_direction(current): 
 i = directions.index(current)
 if i == len(directions) - 1:
  i = 0
 else:
  i += 1
 return directions[i]

max_num = size**2
x = size/2
y = size/2
diagonal_sum = 0
current_direction = (-1,0)
for a in range(1, max_num+1):
 if x == y or x + y == size - 1:
  diagonal_sum += a
 table[x][y] = a
 new_direction = get_next_direction(current_direction)
 x1 = x + new_direction[0]
 y1 = y + new_direction[1]
 if table[x1][y1] == 0:
  x, y = x1, y1
  current_direction = new_direction
 else:
  x, y = x + current_direction[0], y + current_direction[1]

#for row in table:
# print row
print "diagonal_sum", diagonal_sum

Euler #27

https://projecteuler.net/problem=27 => Cuadratic primes

Este problema es muy sencillo, simplemente hay que iterar entre todos los valores posibles para los dos coeficientes (-999 <= a <= 999 y -1000 <= b <= 1000) y comprobar cuántos números primos consecutivos permite calcular el polinomio n² + an + b.

La primera implementación es sencilla y bastante rápida, apenas unos segundos. Pese a que es un problema con un orden de complejidad polinómico por definición (cuadrático para las posibles combinaciones de a y b más otros ciclo para comprobar los números primos), dado que los valores del programa son pequeños, se resuelve en un tiempo razonable.

Pero claro, siempre se quiere optimizar. Así que pensé "si al comprobar si un número es primo o no, puedo tenerlos ya guardados en una lista y así no tengo que calcular lo mismo una y otra vez".

Veamos la implementación, sólo cambia lo siguiente: una lista de números primos inicialmente vacía y una comprobación previa. El resto del código es igual. Sin embargo, al ejecutarlo, nos encontramos con que tarda bastante más tiempo.

¿Qué ha pasado? Fácil, que la supuesta optimización ha resultado ser peor. Con pensar 5 segundos se da uno cuenta por qué: cada vez que queremos saber si un número es primo debemos buscarlo en una lista que cada vez es más grande. Además, debemos asignar memoria y hacer diversas operaciones internas en la lista para redimensionarla al añadir un nuevo elemento. Desconozco la implementación exacta de la búsqueda en listas de Python, pero no puede ser mejor que del orden de O(log n).





Otra forma de crear esta caché de números primos es utilizar un diccionario Python, que requiere menos manejo de memoria, pero la verdad, la mejora, para las dimensiones de este problema, es mínima.

Conclusión: si optimizas, ten cuidado, no sea que lo estropees.

Euler #26

https://projecteuler.net/problem=26 => Reciprocal cycles

Este ha costado un poco más sacarlo, y ha sido gracias al número 1/17 cuando me he dado cuenta que un número sí que se podía repetir en el periodo (asumí que no). Sabiendo eso, la comprobación para ver si estamos empezando un nuevo periodo solo debe tener en cuenta si estamos repitiendo un resto. Veamos una implementación en Python:

from decimal import Decimal, getcontext

getcontext().prec = 40

def Div(D, d):
 return (D/d, D%d)

def GetPeriod(d):
 period = []
 remainders = []
 index = 0
 c, r = Div(1, d)
 if r == 0:
  return []
 if c != 0:
  period.append(c)
 while True:
  D = 10*r
  c, r = Div(D, d)
  if r == 0:
   return []
  period.append(c)
  if r in remainders:  
   return period[period.index(c):index]
  index += 1
  remainders.append(r)

 return period


longestList = []
longestInt = 1
for a in range(1, 1000):
 period = GetPeriod(a)
 if len(period) > len(longestList):
  longestList = period
  longestInt = a

print "Longest period for ", longestInt, " => ", longestList, " (len: ", len(longestList), ")"

Euler #25

https://projecteuler.net/problem=25 => 1000-digit Fibonacci number

Volvemos a intentar cumplir los nuevos propósitos: vamos a probar con un lenguaje que no conocemos: Go.

Primera aproximación, por fuerza bruta, tras varias horas, no obtenemos resultado alguno:

package main

import ("math"
        "fmt")

func fib(n int) int {
    if n <= 2 {
        return 1
    }
    return fib(n-1) + fib(n-2)
}

func main() {
    nDigits := 1000
    limit := int(math.Pow(10, float64(nDigits) - 1)) - 1    
    index := 1
    for {
        num := fib(index)
        if num > limit {
            fmt.Printf("Index => %d\n", index)
            return
        }
        index++
    } 
}

El problema es que cada vez que iteramos, se empieza todo de nuevo. La mejora que vamos a intentar es ir guardando los números de Fibonacci que vamos generando en algún tipo de estructura de datos. Al igual que con C cuando manejamos enteros muy grandes, nos volvemos a encontrar con problemas de desbordamiento, hay que utilizar números enteros de precisión arbitraria (definidos en el package math/big de Go). Al implementarlo, de nuevo, me siento desbordado con los detalles del lenguaje y me doy por vencido, como en el anterior intento, que dejé la implementación en Perl sin hacer.

Implementamos en Python:

fib_numbers_cache = {}

def fib(num):
 if num in fib_numbers_cache:
  return fib_numbers_cache[num]
 if num <= 2:
  fibNum = 1
 else:
  fibNum = fib(num - 1) + fib(num - 2)
 fib_numbers_cache[num] = fibNum
 return fibNum

n = 1
n_digits = 1000
fib_num_max = 10**(n_digits - 1) - 1
while True:
 fn = fib(n)
 if fn > fib_num_max:
  print "END: ", n, "=>", fn
  break
 n += 1

Tras varios intentos de implementar las cosas en distintos lenguajes, siempre vuelvo al que más me gusta, que es Python ;)

Euler #24

https://projecteuler.net/thread=24 => Lexicographic permutations

En el anterior problema me propuse utilizar lenguajes con los que no estuviese muy familiarizado, y este problema lo empecé a hacer con Perl... pero es que no me gusta nada, nada, nada.

Al final, a fuerza bruta y en Python:

def get_permutations (arr):
 if len(arr) < 2:
  return [[arr]]
 elif len(arr) == 2:
  return [[arr[0], arr[1]], [arr[1], arr[0]]]
 else:
  permuts = []
  for a in arr:
   arr2 = arr[:]
   arr2.remove(a)
   for p in get_permutations(arr2):
    p.insert(0, a)
    permuts.append(p)
  return permuts

res = get_permutations([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
print res[999999]

Con el módulo itertools se hace en prácticamente una línea, pero no vale así, es trampa :D

En el foro de Euler comentan varias formas alternativas de afrontar el problema.