lunes, 5 de junio de 2017

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,...

sábado, 27 de mayo de 2017

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)

martes, 11 de abril de 2017

Meld, tremenda herramienta

Hasta ahora había utilizado siempre Meld para apañar merges y conflictos cuando trabajo con Git, pero tiene muchas más utilidades. Acabo de descubrir, tras unos cuantos años, que también compara directorios, y lo hace de una forma muy fácil y vistosa.

Puede mostrar los ficheros que están un directorio y los que no, los que son más recientes, los borrados... y por supuesto, el "side to side" de cualquier fichero.


Por ponerlo en contexto, lo que se muestra es un "theme" de Wordpress que ha sido muy modificado y se le quiere parchear con los nuevos cambios del desarrollador. En verde, a la derecha, nos muestra los ficheros nuevos; tachados, a la izquierda, los que no están presentes; en negrita azul, los que tienen cambios,...

Por supuesto, se puede hacer un diff de los ficheros y editarlos directamente.


sábado, 18 de marzo de 2017

Más conjuntos de Mandelbrot


En la anterior entrada vimos por encima qué era el conjunto de Mandelbrot y cómo se podía generar el gráfico. Recordamos que todo se basaba en la función generadora z2+c. Se puede cambiar la función y repetir el proceso con otras. Estudiemos otras variantes:

  • (z+c)(z-c)
  • z3+c
  • z4+c

Para implementar estas nuevas funciones, dado que hay que hacer algunas operaciones un tanto tediosas, pensé en utilizar alguna librería Javascript que facilitase el trabajo con números complejos. Probé Mathjs, pero es demasiado lenta, así que se ha quedado todo como Javascript plano.


De nuevo, obtenemos unos bonitos fractales.



(z+c)(z-c)


z3+ c



z4+c






En esta página tenemos el código modificado:

jueves, 16 de marzo de 2017

El conjunto de Mandelbrot para 'dummies'




¿Quién no se ha sentido atraído por la belleza de esta figura? La representación gráfica del conjunto de Mandelbrot es el típico ejemplo de figura fractal y es, probablemente, la más popular.


Pero, ¿qué es exactamente, y cómo se llega a esta figura?

Definiendo el conjunto de Mandelbrot


Pertenecen al conjunto de Mandelbrot todos aquellos números complejos que cumplen una propiedad: que el valor de una sucesión determinada, que ahora veremos, para un determinado número, es convergente, es decir, no tiende a infinito.

La sucesión de la que hablamos es recursiva, y  es esta:

zn+1 = z2 + C

donde C es el número que estamos tratando de comprobar si pertenece al conjunto o no y z se va calculando recursivamente, esto es, tomando como entrada de cada iteración el valor obtenido anteriormente.

Recordemos que la aritmética de números complejos involucra operaciones con binomios del tipo a+bi, donde i es la unidad compleja. Recordemos que las operaciones en las que interviene i son un poco especiales, en este caso debemos recordar que i2 es -1.

Vamos a hacer unas cuentas con la ayuda del intérprete de Python, que tiene soporte de serie para números complejos (la unidad imaginaria se denota con j en Python).

Primero, definimos la función que aplicaremos de forma recursiva:
>>> def fm(z, c):
...     return z**2 + c

Ahora, vamos a probar qué pasa con algunos números al aplicarle la función. Empezaremos con el número 1 + i, en notación Python habría que escribirlo como 1 + 1j. Aplicaremos la función anterior cinco veces seguidas, en la primera iteración tanto z como c serán 1 + i y en las siguientes, el resultado de la función se usa como entrada para la siguiente. Es más difícil decirlo que hacerlo:
>>> z = c = 1+1j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(1+3j)
(-7+7j)
(1-97j)
(-9407-193j)
(88454401+3631103j)

Vemos que tras apenas 5 iteraciones, los valores se van disparando, así que, claramente, la sucesión de valores es divergente, se va a infinito.

Probemos con 0.5 + 0.3j:
>>> z = c = 0.5 + 0.3j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(0.66+0.6j)
(0.5756+1.092j)
(-0.36114864+1.5571104j)
(-1.79416445761-0.82469660658j)
(3.03890160806+3.25928267968j)

También parece que diverge, pero no tan rápido como el caso anterior.

Otro caso más, con el número 0.2 - 0.1j:


>>> z = c = 0.2 - 0.1j
>>> for x in range(5):
...     z = fm(z, c)
...     print z
... 
(0.23-0.14j)
(0.2333-0.1644j)
(0.22740153-0.17670904j)
(0.220485371029-0.180367812122j)
(0.216081251188-0.179536927955j)

Tras cinco iteraciones, no parece que se vaya a infinito. Probemos otras cinco iteraciones:


>>> for x in range(10):
...     z = fm(z, c)
...     print z
... 
(0.214457598616-0.177589128054j)
(0.214454163201-0.176170675885j)
(0.214954481072-0.175561069755j)
(0.21538373972-0.175475277291j)
(0.215598582395-0.175589042903j)
(0.215651236743-0.175713497468j)
(0.215630222717-0.175785666083j)
(0.215595792549-0.175809404656j)
(0.215572598999-0.175807535868j)
(0.215563255771-0.175798574862j)

 Parece claro que este número, 0.2 - 0.1j, al aplicarle la función del principio, la sucesión que se obtiene es convergente.

Se dice que un número pertenece al conjunto de Mandelbrot cuando la sucesión es convergente, es decir, no se va a infinito. Se puede probar que en cuanto que uno de los resultado de la sucesión se sale del rango de ± 2 ± 2i, la sucesión será divergente, esto es, el número no pertenecería al conjunto. La forma más rápida de calcular esto es comprobando, para cada uno de los resultados intermedios, si el módulo del número complejo es > 2, o lo que es lo mismo, si la suma de los cuadrados de sus partes real e imaginaria es > 4. De esta última forma llegamos al mismo resultado ahorrándonos hacer una raíz cuadrada, que siempre es más costosa computacionalmente.

Sabiendo esto, vamos a empezar a comprobar todos los posibles números imaginarios que hay entre -2 - 2i y 2 + 2i, y viendo si la sucesión que se deriva de cada uno de ellos diverge. Obviamente, hay infinitos números en esta región, así que tenemos que ir calculando el mayor número de ellos con pequeños incrementos.

De la misma forma, tendremos que hacer unas cuantas iteraciones para saber si la sucesión diverge o no. Por ejemplo, para el número 0.5 + 0.3j que vimos antes tuvimos que hacer 5 iteraciones o calcular 5 términos de la sucesión antes de darnos cuenta de que realmente divergía.

Empecemos de forma sencilla: moviéndonos en pasitos de 0.1 y haciendo 3 iteraciones. Tendríamos esta representación gráfica:



Vemos que tenemos muy poco nivel de detalle, incrementemos paulatinamente el nivel de detalle (haciendo los incrementos más pequeños y aumentando el número de iteraciones para "cazar" más números que están fuera del conjunto), ahí van unas cuantas muestras:





En este punto, ya tenemos un detalle bastante fino, pero es probable que, como hemos hecho pocas iteraciones, haya números que pensemos que pertenecen al conjunto, pero en realidad, no pertenecen. Aumentemos las iteraciones.





Hasta ahora nos hemos limitado a colorear en rojo los números que, al aplicar la función, la sucesión diverge; en negro, los que converge.

El siguiente paso en detalle, y aquí es donde aparece la "magia" de este conjunto, es asignar diferentes colores a la velocidad con la que los números que no pertenecen al conjunto, al evaluar la sucesión, divergen. Esto es, si asignamos colores diferentes al número de iteraciones necesarias para determinar si el número no pertenece al conjunto, empezamos a obtener unos bonitos gráficos.

La gracia en este caso es encontrar una gradación de colores y asignación al número de iteraciones que quede bonita. Por ejemplo, una sencilla, basada en amarillos:



Un último ejemplo, con mayor resolución, utilizando una gama más o menos aleatoria de colores:


Implementación


Es bastante sencillo implementar una representación gráfica del fractal. Esta que se muestra está hecha en HTML Canvas y Javascript, usando "plain Javascript", sin plugins ni librerías.

Lo primero que debemos hacer es definir un canvas en nuestro HTML. También hemos agregado unos input para poder jugar con el incremento y el número de iteraciones a probar.


<canvas id="canvas"></canvas>
<form>
<p>Increment: <input name="increment" value="0.01"></p>
<p>Iterations: <input name="iterations" value="1000"></p>
<p><input type="button" value="Run!" id="run"></p>
<p id="message"></p>
</form>

A continuación, escribimos el código Javascript que hará los cálculos y representará los puntos en pantalla. El primer paso, es acceder e inicializar el elemento canvas:


var width = 600;
var height = 600;
var c = document.getElementById("canvas");
c.setAttribute('width', width);
c.setAttribute('height', height);
var ctx = c.getContext("2d");

Ahora definimos algunas constantes y funciones:
  • adjust_scale(i, j) => centra los puntos correctamente en el gráfico. Nuestro sistema de coordenadas es el mátemático (0, 0) en el centro, en un canvas el (0, 0) es la esquina superior izquierda.
  • get_color(nIter) => devuelve un color en función del número de iteraciones. En esta función se puede ser creativo y experimentar para probar diferentes gamas de colores y gradaciones. En nuestro ejemplo nos limitamos a una gama de rojos.
  • mandelbrot(z, c) => calcular un término de la sucesión aplicando la fórmula zn+1 = z2 + C, ya explicada.
  • in_mandelbrot_set(c, maxIter) => determina si un número C está en el conjunto de Mandelbrot utilizando un máximo de iteraciones. Devuelve un objeto con dos claves: si la sucesión es convergente y el número de iteraciones necesarias para llegar a dicha conclusión.
  • run() => lanza el gráfico.

El código completo es este:

<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<canvas id="canvas"></canvas>
<form>
<p>Increment: <input name="increment" value="0.01"></p>
<p>Iterations: <input name="iterations" value="1000"></p>
<p><input type="button" value="Run!" id="run"></p>
<p id="message"></p>
</form>

<script>
var width = 600;
var height = 600;
var c = document.getElementById("canvas");
c.setAttribute('width', width);
c.setAttribute('height', height);
var ctx = c.getContext("2d");

var mandelbrot_limit = 1.8;
var mandelbrot_max_module = 4;



function adjust_scale(i, j) {
    return [(width/(2*mandelbrot_limit))*i + (width/2), (height/(2*mandelbrot_limit))*j + (height/2)];
}

function get_color(nIter) {
    var factor = 50 + nIter;
    factor = parseInt(factor);
    if (factor > 255) factor = 255;
    return "rgb("+factor+",0,0)";
}

function mandelbrot(z, c) {
    var x = z[0]; 
    var y = z[1];
    var real = x*x - y*y + c[0];
    var imag = 2*x*y + c[1];
    return [real, imag];
}

function in_mandelbrot_set(c, maxIter) {
    var x = 0;
    var res = c
    while (x < maxIter) {
        res = mandelbrot(res, c);
        if ((res[0]*res[0] + res[1]*res[1]) > mandelbrot_max_module) {
            return {'status': false, 'nIter': x}; // No está en el conjunto, hemos tardado x interaciones en saberlo
        }
        x++;
    }
    return {'status': true, 'nIter': x};
}

function run() {
    ctx.fillStyle = "#FFFFFF";
    ctx.fillRect(0, 1, width, width);

    var mandelbrot_inc = parseFloat(document.forms[0].elements.increment.value);
    var mandelbrot_test_iterations = parseInt(document.forms[0].elements.iterations.value);
    var point_width = 1;
    if (mandelbrot_inc >= 0.0009) point_width = 2;
    if (mandelbrot_inc >= 0.009) point_width = 8;       
    if (mandelbrot_inc >= 0.09) point_width = 16;       

    document.getElementById("message").innerHTML = "Increment: "+mandelbrot_inc+", Iterations: "+mandelbrot_test_iterations;
    for(var i=-mandelbrot_limit; i<=mandelbrot_limit; i += mandelbrot_inc) {
        for (var j=-mandelbrot_limit; j<=mandelbrot_limit; j += mandelbrot_inc) {
            check = in_mandelbrot_set([i, j], mandelbrot_test_iterations);
            if (check.status) {
                ctx.fillStyle = "#000000";
            } else {
                ctx.fillStyle = get_color(check.nIter);
            }
            var point = adjust_scale(i, j);
            ctx.fillRect(point[0], point[1], point_width, point_width);
        }
    }
}

document.getElementById("run").addEventListener('click', run);

</script>