Sumando subconjuntos – Un problema

Julio 15, 2010 at 5:51 pm (matemática)

Un breve problema para romper el largo intervalo sin posts:

En el siguiente conjunto de 70 números hay dos subconjuntos (distintos!) con igual suma.


5213588008354709077 9011219417469017946
9115813602317594993 1796745334760975384
2375118318707634486 3579709154995762145
2312952310307873817 3627590721354606439
5941472576423725122 4317696052329054505
5763055694478716846 5226146329630268999
2730952213106576953 6235647640047602135
5014371167157923471 4868653895375087301
9737387704190733086 9262565481673300485
5968266991171521006 6752113930489992229
5917882775011418152 5866436908055159779
9233099989955257688 9376338948688351159
3772194655144630358 9029836747496103755
3318990262990089104 9205546877037990475
9498032114812022495 9990105869964955299
9849598364470384044 2664344232060524242
1376783969573792128 1108556560089425769
7820574110347009988 1603454130529647419
6951628222196921724 9889945707884382295
4776591302180789869 7652184907611215542
3982809542974920136 8082643935949870308
1233783467857668500 4271233363607031147
7999940522596325715 1509145775275864006
6415171202616583365 4143260390225524497
6393762352694839041 2290598705560799669
7835010686462271800 8505865989334316749
7547888419188030222 4408107167048106771
8998470433081591390 9131522681668251869
9096632376298092495 7481066850797885510
5295758362772474604 1796212605533609345
5953431042043343946 3151838989804308537
6445353832659195698 9929081270845451034
6422250272800292361 8643312627450063997
7207546018039041794 3624820335477016277
3782849199070331143 6012991563467874119

Como podemos estar seguros de ello? Espero sus respuestas en los comentarios y daré la mía en el próximo post.

Permalink 2 comentarios

Ante la duda, aplicar Fourier

Junio 19, 2010 at 10:15 pm (física, matemática, python)

Introducción

Un clásico problema para explicar el poder del principio de superposición y de la simetría consiste en calcular la resistencia entre los puntos A y B de la siguiente red infinita de resistores:

Red infinita de resistores unitarios con dos nodos marcados.

La forma más simple de resolverlo es asumir que inyectamos una corriente de 1 por A, dejando desconectados los otros nodos. Por simetría, las corrientes que partan del nodo serán todas iguales a 1/4:

Se inyecta una corriente unitaria en el nodo A de la red anterior.

Por linealidad y simetría, si repetimos el procedimiento inyectando una corriente de -1 por el nodo B, las corrientes que partan del mismo serán iguales a -1/4. Aplicando el principio de superposición (posible por la linealidad de la Ley de Ohm), llegamos a las siguientes corrientes:

Superposición de los efectos de inyectar una corriente unitaria por el nodo A y extraer una corriente unitaria por el nodo B.

Como esto equivale a inyectar una corriente de 1 por A y extraerla por B, la diferencia de potencial entre ambos nodos será idéntica a la resistencia entre estos. La corriente que circula por el resistor unitario que los une es 1/2 y, por lo tanto, la resistencia entre A y B es 1/2.

Una variante de este problema (considerablemente más compleja) que ha aparecido en múltiples lugares (incluso XKCD) es la siguiente:

Red de resistores unitarios con los nodos A y B separados por un 'salto de caballo'.

En lo que resta del post veremos como obtener una solución numérica y, finalmente, una analítica a este problema.

Una solución numérica

Un nodo cualquiera de la red de resistores debe satisfacer la ley de Kirchhoff de conservación de corrientes. Como las corrientes a través de cada resistor unitario son iguales a la diferencia de tensión, tenemos:

(V_{i\,j} - V_{i-1\,j}) + (V_{i\,j} - V_{i+1\,j}) + (V_{i\,j} - V_{i\,j-1}) + (V_{i\,j} - V_{i\,j+1}) = 0

4V_{i\,j} - V_{i-1\,j} - V_{i+1\,j} - V_{i\,j-1} - V_{i\,j+1} = 0.

Estas ecuaciones se aplicarán a todos los nodos a excepción de los que estén conectados a fuentes externas. Si reemplazamos la red infinita de resistores por una red finita con N2 nodos, la aplicación de la ecuación de nodos nos dejará N2 – 2 ecuaciones, ya que la ecuación de nodos no se aplicará donde conectemos la fuente. Esta cantidad de ecuaciones es suficiente para determinar todas las tensiones, ya que las tensiones de dos nodos son forzadas por la fuente.

Utilizando una topología toroidal por simplicidad y aplicando Gauss Seidel + SOR con ω = 1.5 llegamos al siguiente código (ver resultados en codepad):

TOL = 1e-6
K = 1.5

def solve(n):
    assert n >= 4
    v = [[0.5 for y in range(n)] for x in range(n)]
    delta = 2 * TOL
    FIXED_VALUES = {
        (0, 0): 1,
        (2, 1): 0
    }
    while delta > TOL:
        delta = 0
        for x in range(n):
            for y in range(n):
                if (x, y) in FIXED_VALUES:
                    v[x][y] = FIXED_VALUES[(x, y)]
                    continue
                vo = v[x][y]
                vn = (v[x - 1][y] + v[(x + 1) % n][y] +
                      v[x][y - 1] + v[x][(y + 1) % n]) / 4.0
                v[x][y] = vn * K + vo * (1 - K)
                delta = max(abs(vn - vo), delta)
    return 1.0 / (4 * v[0][0] - v[-1][0] - v[1][0] - v[0][-1] - v[0][1])

if __name__ == '__main__':
    from time import time
    for n in (8, 16, 32, 64, 128):
        start = time()
        print 'N = %d - R = %f - Time: %f s' % (n, solve(n), time() - start)

Ejecutando el código completo (la versión de codepad excluye N = 64 y N = 128 por las limitaciones al tiempo de ejecución) obtenemos la siguiente salida:

N = 8 - R = 0.734997 - Time: 0.000000 s
N = 16 - R = 0.763460 - Time: 0.094000 s
N = 32 - R = 0.770537 - Time: 0.781000 s
N = 64 - R = 0.772408 - Time: 2.532000 s
N = 128 - R = 0.773033 - Time: 23.687000 s

Se observa una convergencia relativamente rápida a un valor cercano a 0.773 y un comportamiento bastante extraño del tiempo de ejecución.

Una solución analítica

Si llamamos Vi j a la tensión del nodo m, n e Im n a la corriente inyectada en el mismo, la conservación de la corriente nos lleva a lo siguiente:

I_{m\,n} = 4V_{m\,n} - V_{m-1\,n} - V_{m+1\,n} - V_{m\,n-1} - V_{m\,n+1}.

Las corrientes inyectadas pueden representarse como

I_{m\,n} = \delta_{m\,n} - \delta_{m-M\,n-N},

donde las “deltas” inyecta una corriente unitaria en el nodo 0, 0 y la extraen por el nodo M, N. Esta ecuación en diferencias no tiene una solución única por esencialmente el mismo motivo que el potencial de una carga puntual no tiene una única solución: el núcleo del laplaciano es no trivial. Pero hay una solución única si exigimos que el potencial se vaya a cero cuando m y n se van a infinito.

Un método muy utilizado para resolver ecuaciones diferenciales es aplicar una transformación inversible que convierta a la ecuación diferencial en una simple ecuación algebraica. Si esta última puede resolverse, aplicando la transformación inversa obtendremos una solución en el dominio original.

Si aplicamos la “versión adecuada” de la transformada de Fourier en dos dimensiones,

\displaystyle V(p,q) = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i m p} e^{-i n q},

podemos ver que la transformada de Vm-1 n es

\displaystyle \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m-1\,n} e^{i m p} e^{i n q} = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i (m + 1) p} e^{-i n q}

\displaystyle = e^{-i p} \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i m p} e^{-i n q}

\displaystyle = e^{-i p} V(p, q).

Utilizando esto, llegamos a la siguiente ecuación algebraica sobre V(p, q):

\displaystyle 1 - e^{-i M p}e^{-i N q} = V(p, q) (4 - e^{-i p} - e^{i p} - e^{-i q} - e^{i q}).

Resolviéndola:

\displaystyle V(p, q) = \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}.

Antitransformando:

\displaystyle V_{m\,n} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}e^{i m p}e^{i n q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}e^{i m p}e^{i n q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{e^{i m p}e^{i n q} - e^{i (m - M) p}e^{i (n - N) q}}{4 - 2 \cos p - 2 \cos q}.

Pero a nosotros nos interesa la resistencia entre los nodos 0, 0 y M, N, igual a la diferencia de tensión entre ambos nodos (ya que circula entre ellos una corriente unitaria). Por lo tanto,

\displaystyle R_{M\,N} = V_{0\,0} - V_{M\,N}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q} - e^{i M p}e^{i N q} + 1}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i (M p + N q)} - e^{i (M p + N q)} + 1}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{2 - 2\cos(M p + N q)}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(M p + N q)}{2 - \cos p - \cos q}

Ahora “solo” resta resolver esta integral para cuando M = 2 y N = 1:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(2 p + q)}{2 - \cos p - \cos q}.

Evaluándola numéricamente obtenemos 0.7732, un valor que coincide en 3 cifras significativas con el obtenido previamente.

Solución analítica de la integral

ADVERTENCIA: gran cantidad de fórmulas que no aportan mucho…

Podemos empezar buscando una solución via Wolfram Alpha, pero no funciona. Entonces comenzamos expandiendo el numerador:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(2 p + q)}{2 - \cos p - \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos 2p \cos q + \sin 2p \sin q}{2 - \cos p - \cos q}

Trabajando con el primer término,

\displaystyle A = \frac{1}{2 - \cos p}

\displaystyle \int_0^{2\pi} dq \frac{1}{2 - \cos p - \cos q} = \frac{1}{A^{-1}}\int_0^{2\pi} dq \frac{1}{1 - A\cos q},

llegamos a una forma en tabla (Tabla 64, Formula 12):

\displaystyle = 2A\int_0^{\pi} dq \frac{\cos 0\cdot q}{1 + (-A)\cos q}

Aplicando la fórmula y simplificando:

\displaystyle = 2A\frac{\pi}{\sqrt{1-A^2}}\left(\frac{\sqrt{1-A^2}-1}{-A}\right)^0

\displaystyle = \frac{2}{2 - \cos p}\frac{\pi}{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}}

\displaystyle = \frac{2}{2 - \cos p}\frac{\pi(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}}

\displaystyle = \frac{2\pi}{\sqrt{(2 - \cos p)^2-1}}

Del mismo modo podemos proceder con el segundo término:

\displaystyle \int_0^{2\pi} dq \frac{\cos 2p \cos q}{2 - \cos p - \cos q} = \frac{\cos 2p}{A^{-1}}\int_0^{2\pi} dq \frac{\cos q}{1 - A\cos q}

\displaystyle = A \cos 2p\int_0^{2\pi} dq \frac{\cos q}{1 - A\cos q}

\displaystyle = 2A\cos 2p\int_0^{\pi} dq \frac{\cos 1\cdot q}{1 + (-A)\cos q}

\displaystyle = 2A\cos 2p\frac{\pi}{\sqrt{1-A^2}}\left(\frac{\sqrt{1-A^2}-1}{-A}\right)^1

\displaystyle = -2A\cos 2p\frac{\pi}{\sqrt{1-A^2}}\frac{\sqrt{1-A^2}-1}{A}

\displaystyle = -2\frac{1}{2 - \cos p}\cos 2p\frac{\pi}{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}}\frac{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}-1}{\frac{1}{2 - \cos p}}

\displaystyle = -\frac{2\pi\cos 2p}{(2 - \cos p)\frac{\sqrt{(2 - \cos p)^2-1}}{2 - \cos p}}\frac{\frac{\sqrt{(2 - \cos p)^2-1}}{2 - \cos p}-1}{\frac{1}{2 - \cos p}}

\displaystyle = -\frac{2\pi\cos 2p}{\sqrt{(2 - \cos p)^2-1}}\left(\sqrt{(2 - \cos p)^2-1} - (2 - \cos p)\right)

\displaystyle = - 2\pi\cos 2p + \frac{2\pi\cos 2p(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}}

La integral del último término se anula, como puede verse a continuación:

\displaystyle \int_0^{2\pi} dq \frac{\sin 2p \sin q}{2 - \cos p - \cos q} = \frac{1}{A^{-1}}\int_0^{2\pi} dq \frac{\sin 2p \sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^{2\pi} dq \frac{\sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} - A\sin 2p \int_{2\pi}^\pi dq \frac{\sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} + A\sin 2p \int_0^\pi dr \frac{\sin (2\pi - r)}{1 - A\cos (2\pi - r)}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} + A\sin 2p \int_0^\pi dr \frac{-\sin r}{1 - A\cos r}

\displaystyle = 0

Mediante estas integraciones, la integral doble original

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos 2p \cos q + \sin 2p \sin q}{2 - \cos p - \cos q}

se reduce a una integral simple:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \frac{2\pi}{\sqrt{(2 - \cos p)^2-1}} - \frac{2\pi\cos 2p(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}} + 2\pi\cos 2p

\displaystyle = \frac{1}{2\pi} \int_0^{2\pi} dp \frac{1- 2\cos 2p + \cos p\cos 2p}{\sqrt{(2 - \cos p)^2-1}}.

Utilizando Wolfram Alpha:

\displaystyle R_{2\,1} = \frac{1}{2\pi} \int_0^{2\pi} dp \frac{1- 2\cos 2p + \cos p\cos 2p}{\sqrt{(2 - \cos p)^2-1}}

\displaystyle = \frac{1}{2\pi} (8-\pi)

\displaystyle = \frac{4}{\pi} - \frac{1}{2}

\displaystyle \approx 0.7732,

coincidiendo con el resultado aproximado anterior.

Conclusiones y reconocimientos

Esto vindica, como Fer bien conoce, la importancia de aplicar Fourier cuando uno se encuentra frente a un problema que no sabe como resolver. :-P

Soluciones esencialmente idénticas a esta fueron desarrolladas mucho antes por las fuentes citadas y por “frooha”, que desarrolló una solución muy similar a esta, entre otros. También existe al menos una solución por la via experimental, preferida por los electrónicos.

Permalink Dejar un comentario

Dinámica de fluidos en Javascript

Mayo 30, 2010 at 8:51 pm (CFD, física, javascript)

Inspirado por un ejemplo portado por Oliver Hunt de C a JS (y frustrado porque no se me ocurría como resolver el nivel Ophanim en Manufactoria :-D ), decidí portar mi código de simulación de fluidos mediante Lattice Boltzmann (LB) a Javascript y agregarle interactividad. La visualización es algo primitiva, ya que me concentré en optimizar los cálculos sin prestarle mucha atención a la optimización del manejo de canvas; posiblemente mejore ese aspecto en una próxima versión.

El código presenta algunas limitaciones propias de utilizar una grilla de 64 X 64 como base (el máximo tamaño con performance “realtime“) y, muy probablemente de mi desconocimiento sobre como optimizar los parámetros de LB. Es particularmente clara la compresibilidad del fluido y no es muy difícil hacer “explotar” la simulación (se recupera mediante un reset).

Como WordPress no permite incluir código Javascript en los posts, para visualizar la simulación pueden hacer click en la imagen de referencia de los controles.

Referencia de los controles de la simulación. Hacer click para arrancar la simulación.

El uso de los controles es muy simple:

  • Play/pause: arranca o detiene la simulación.
  • Reset: reinicia la simulación desde cero.
  • Visualization reset: reinicia la visualización, creando una nueva serie de partículas para marcar el movimiento del fluido (sin alterar el estado del mismo).

Para darle un impulso a una parte del fluido basta con arrastrar sobre el área de visualización.

En un posterior post incluiré detalles acerca de como funciona LB y de las optimizaciones realizadas (tengo que leer más, sobre todo respecto a la conexión con los valores físicos :-D ).

Permalink 4 comentarios

Porqué no se puede ganar siempre

Mayo 2, 2010 at 8:15 pm (informática, javascript, matemática, probabilidad)

Respuesta al post anterior

La técnica mostrada en el post anterior se denomina martingala y es muy antigua. Obviamente, a pesar de lo que puedan decir algunos, no funciona. Pero porqué la simulación anterior parecía mostrar lo contrario?

El motivo principal del fallo de la simulación es no limitar la cantidad de dinero del jugador disponible para apostar. Como la cantidad que se requiere apostar crece exponencialmente con la cantidad de jugadas, en una cantidad relativamente pequeña de estas se excederán las reservas del jugador y no podrá continuar apostando.

Otro problema de la simulación es que, incluso una cantidad grande de simulaciones, puede obviar acontecimientos muy improbables. Esto no es problemático cuando los eventos en cuestión no tienen consecuencias excepcionales pero, cuando conducen a un impacto proporcional a su improbabilidad, puede distorsionarse mucho el valor esperado (probar la nueva simulación Javascript con M = 100000 para N = 1000 y N = 1000000; observar las diferencias en los valores esperados estimados).

Pero puede objetarse (con razón!) que solo demostramos el fallo de una estrategia de apuesta, no de todas las posibles. Por lo tanto, a continuación, demostraremos que esto se aplica a cualquier estrategia.

Porqué los sistemas no funcionan

Como en toda demostración matemática, necesitaremos hacer ciertas suposiciones sobre el problema en cuestión:

  1. El juego consiste en una cierta serie de jugadas, cada una asociada con una variable aleatoria X_i. Asumimos que todas las variables X_i tienen la misma distribución.
  2. En la jugada i el jugador apuesta una cantidad m_i a un evento aleatorio e_i (donde tomamos e_i como un subconjunto de lso valores que podría tomar X_i).
  3. Si el evento ocurre al que el jugador apostó ocurre, este recibe el monto apostado multiplicado por una constante k(e_i) dependiente del evento al que apostó. En caso contrario, pierde lo apostado.
  4. k(e) P(e) \le 1.
  5. Los valores de las variables aleatorias X_1, X_2, ..., X_{i-1} no aportan información sobre el valor de X_i.

Las primeras 3 suposiciones son simplemente una descripción general de un juego de azar, mientras que 5 nos indica que las anteriores jugadas no nos dan información sobre la jugada actual. La suposición 4 parece más compleja, pero esencialmente nos indica que los factores k han sido correctamente elegidos por el casino.

Para ver esto, supongamos que existe un evento e tal que k(e) P(e) > 1. Entonces, si elegimos una estrategia de apostar N veces al evento e un monto unitario, nuestro balance esperado será:

\displaystyle E[b] = E[\sum_{i=1}^N b_i]

\displaystyle = \sum_{i=1}^N E[b_i]

\displaystyle = \sum_{i=1}^N ((k(e)-1)P(e) - (1-P(e)))

\displaystyle = N (k(e)P(e) - 1)

\displaystyle > 0,

donde b es el balance general del jugador y b_i el balance de la jugada i.

Eso implicaría que podríamos ganar simplemente repitiendo una apuesta una y otra vez (puede demostrarse sin mucha dificultad que ni siquiera necesitaríamos reservas muy grandes, solo de orden \sqrt(N)). Como eso claramente no ocurre en los casinos, es obvio que los multiplicadores están elegidos correctamente (de hecho, por ejemplo en la ruleta, están elegidos “a favor” del casino por la existencia del cero).

Busquemos ahora el valor esperado de una estrategia de apuestas m_i(i, X_1, X_2, ..., X_{i-1}), e_i(i, X_1, X_2, ..., X_{i-1}):

\displaystyle E[b] = E[\sum_{i=1}^N b_i]

\displaystyle = E[\sum_{i=1}^N b_i]

\displaystyle = \sum_{i=1}^N E[b_i]

\displaystyle = \sum_{i=1}^N ((k(e_i)-1)P(e_i) - (1-P(e_i)))m_i

\displaystyle = \sum_{i=1}^N (k(e_i)P(e_i) - 1)m_i

siendo e_i y m_i las funciones anteriormente descritas. Pero la suposición 5 nos dice que

P(e_i|i,X_1=x_1, X_2=x_2, ..., X_{i-1} = x_{i-1}) = P(e_i),

es decir que esa información no nos modifica la probabilidad de ningún evento y, en particular, de e_i. Combinando con la suposición 4 llegamos a

\displaystyle E[b] \le 0 ,

siempre que m_i sea positivo, como es razonable.

Caminos de escape

Aunque una demostración fuera perfecta, sus conclusiones nunca será más fuertes que sus hipótesis. En esta sección analizaremos algunas formas en las que las hipótesis pueden ser falsas:

mi < 0: Este es esencialmente el método utilizado por los casinos. Suele ser bastante complejo de seguir, pero muy redituable :-)

Fallo en la hipótesis 4: Los multiplicadores son fáciles de determinar conociendo las probabilidades y generalmente se utilizan probabilidades “nominales” (por ejemplo, suponer equiprobables a los 37 o 38 números de la ruleta). Pero si algún observador encuentra evidencia que lo mueva a probabilidades distintas de las nominales, puede que existan apuestas con valor esperado positivo.

Fallo en la hipótesis 5: Pueden existir correlaciones no esperadas entre las variables aleatorias. Esto es más común en los juegos por computadora, que suelen usar números pseudoaleatorios y, en algunos casos, tienen errores de implementación que vuelven predecibles las secuencias generadas.

Permalink Dejar un comentario

Como ganar siempre :-)

Abril 29, 2010 at 12:57 am (C++, informática, javascript, probabilidad, python)

Un sistema de apuestas, muy utilizado por lo que puede verse en Google, sigue el siguiente esquema (donde lo aplicamos a una apuesta doble o nada con una moneda “justa”):

  1. Empezar apostando una unidad.
  2. Apostar a “cara”.
  3. En caso de ganar, retirarse.
  4. En caso de perder, doblar la apuesta y volver al paso 2.

Analicemos este sistema mediante un millón de simulaciones de la técnica, utilizando 3 lenguajes de programación: C++ (codepad), Python (codepad, N = 1000) y Javascript.

Como puede observarse, los resultados son esencialmente idénticos:

Number of simulations: 1000000
Minimum player balance: 1
Maximum player balance: 1
Average number of rounds before winning: 2.000831

En todas las ejecuciones observamos una ganancia fija de una unidad y que se obtiene en promedio en una pequeña cantidad de tiradas.

Donde está el problema con este sistema? En el próximo post la respuesta… :-D

Permalink 2 comentarios

Detectando… paralelepípedos

Abril 26, 2010 at 3:46 am (informática, matemática, python)

(Continuación de Detectando los cubos.)

En esta primera parte veremos como detectar configuraciones de puntos que pueden formar la proyección paralela de los vértices de un paralelelepípedo, resolviendo parte del problema planteado anteriormente. Analizaremos como determinar si corresponden a un cubo en un posterior post.

La idea básica detrás del algoritmo de detección es observar que las caras de un paralelepípedo son paralelogramos y que estos preservan el paralelismo de sus lados al ser proyectados en forma paralela.

Puede dividirse a este algoritmo en cinco pasos:

  1. Buscar el punto que se encuentre más a la izquierda (menor valor de X) y, dentro de ellos más arriba (menor valor de Y, en nuestro caso). Ese punto será denominado p0 y será tomado como referencia.
  2. Imaginando un semieje en dirección +X partiendo desde p0, se eligen los dos puntos que posean el mayor y el menor ángulo (teniendo en cuenta el signo!) respecto a ese eje. Estos puntos son denominados p1 y p2 y formarán la base para identificar la primera cara del paralelepípedo.
  3. Se busca un punto p12 que cierre el paralelogramo formado por p0, p1 y p2. Estos puntos formarán la que podemos denominar “cara -X” del paralelepípedo.
  4. Opuesta a la “cara -X” debe existir una “cara +X”, que además debe poseer vértices “análogos” a cada uno de los hallados anteriormente, solo que trasladados en forma paralela. Podemos buscar el análogo de p0 repitiendo el procedimiento utilizado para hallar p0 pero omitiendo los puntos p0, p1, p2 y p12.
  5. Una vez hallado p3, conocemos el desplazamiento relativo entre ambas caras (-X y +X). Por lo tanto podemos buscar los “análogos” a p1, p2 y p12: p13, p23 y p123. Si todos estos puntos pueden ser hallados, entonces la serie de puntos puede ser la proyección paralela de un paralelepípedo.

En la siguiente figura podemos ver el funcionamiento del algoritmo:

Funcionamiento normal del algoritmo (detección exitosa).

Los únicos pasos donde el algoritmo puede fallar, detectando que no estamos trabajando con la proyección de un paralelepípedo, son los pasos 3 y 5:

Funcionamiento con fallo en la detección del algoritmo (paso 3).

Funcionamiento con fallo en la detección del algoritmo (paso 5).

El código correspondiente en Python, que es una traducción casi directa de los pasos mencionados anteriormente, puede encontrarse en el repositorio SVN de este blog. La única diferencia digna de mención es la adición de un “paso 0″ que comprueba la cantidad de vértices.

Para probarlo basta con ejecutar el generador de casos de prueba (que creará dos archivos: input.txt y key.txt), luego el detector de paralelepípedos (que creará un archivo de salida output.txt) y, por último, el comparador (que comparará el archivo output.txt y el archivo key.txt, indicando la cantidad de falsos negativos y falsos positivos).

Permalink Dejar un comentario

Otra forma de pensar Taylor

Marzo 11, 2010 at 5:17 pm (matemática)

En primer lugar, en el post anterior había dejado planteado el siguiente problema:

Dado un conjunto de 51 números dentro del intervalo 1-100 (incluyendo los extremos), demostrar que existen al menos dos números en ese conjunto que no comparten ningún divisor (a excepción de 1, obviamente).

Así que, antes de empezar con el tema en sí de este post, no puedo hacer menos que presentar mi demostración ;-)

El problema de los 51 números

La prueba que pensé se divide en dos pasos:

  • Paso 1: Determinar que al menos dos de los 51 números deben ser adyacentes.
  • Paso 2: Determinar que dos números adyacentes no comparten divisores no triviales.
Paso 1

Es intuitivamente razonable que, dados 51 números en un intervalo de 100, al menos dos deben ser adyacentes. Pero como la demostración que encontré no es muy elegante, la separé en una página auxiliar. (Dejo la primera demostración como ejemplo de cuanto más complejo de lo necesario puede hacerse algo…)

Definamos 100 variables 0-1 x1, x2, …, x100 para representar el subconjunto en cuestión. Claramente la condición sobre la cantidad de elementos puede representarse como

\displaystyle \sum_{i=1}^{100} x_i = 51

Si agrupamos la suma en 50 pares obtenemos

\displaystyle \sum_{i=0}^{49} x_{2i+1} + x_{2i+2} = \sum_{i=0}^{49} y_i = 51

Esto implica que el valor medio de los yi es 51/50 = 1.02 > 1. Pero como los yi solo pueden tomar los valores enteros 0, 1 y 2, eso implica que existe al menos un yi que vale 2. En consecuencia, existen necesariamente en el subconjunto dos números adyacentes.

Paso 2

Es bien conocido, y la base del algoritmo de Euclides, que

\mathrm{gcd}(a, b) = \mathrm{gcd}(a, b \% a),

donde \% se utiliza para indicar reducción modular al estilo C.

Ahora, si dos números a y b son adyacentes, entonces podemos decir sin pérdida de generalidad que b = a + 1. Por lo tanto:

\mathrm{gcd}(a, b) = \mathrm{gcd}(a, a+1) = \mathrm{gcd}(a, (a+1) \% a) = \mathrm{gcd}(a, 1) = 1.

Por ende, dos números adyacentes no pueden tener factores en común.

Conclusión

Como los 51 números en {1, 2, …, 100} deben incluir al menos un par de números adyacentes y los números adyacentes no tienen factores en común, entonces los 51 números en {1, 2, …, 100} deben incluir al menos un par de números sin factores comunes.

Taylor

Las demostraciones normales del Teorema de Taylor siempre me parecieron poco motivadas: son indudablemente correctas, pero no queda claro porqué uno iría por ese camino. Hace unos años encontré una demostración que, aunque no me cabe duda que debe haber sido conocida por los egipcios :-P , nunca vi en otro lado.

Para simplificar realicemos la expansión alrededor de cero (serie de Maclaurin). No cambia los resultados pero disminuye el “ruido” en las expresiones.

La idea básica es utilizar el hecho de que toda función derivable f(x) puede ser representada en función de su valor en 0 y su derivada del siguiente modo:

\displaystyle f(x) = f(0) + \int_0^x dx_1 f'(x_1)

Si la función puede ser derivada múltiples veces, el proceso puede ser aplicado en forma recursiva:

\displaystyle f'(x) = f'(0) + \int_0^x dx_1 f''(x_1)

\displaystyle f''(x) = f''(0) + \int_0^x dx_1 f'''(x_1)

\displaystyle \ldots

Combinando estas representaciones podemos representar a una función en término de su valor, el de las n-1 primeras derivadas en cero y su derivada n-ésima:

\displaystyle f(x) = f(0) + \int_0^x dx_1 \left( f'(0) + \int_0^{x_1} dx_2 \left( f''(0) + \ldots \right) \right)

Para poder simplificar estas expresiones es necesario resolver integrales de la forma

\displaystyle \int_0^x dx_1 \int_0^{x_1} dx_2 \ldots \int_0^{x_{n-1}} dx_n.

Esto podría hacerse por inducción pero, como estoy un poco quemado de utilizar esa técnica, podemos hacerlo por inspección simplemente:

\displaystyle \int_0^x dx_1 = x

\displaystyle \int_0^x dx_1 \int_0^{x_1} dx_2 = \int_0^x dx_1 x_1 = \frac{x^2}{2}

\displaystyle \int_0^x dx_1 \int_0^{x_1} dx_2 \int_0^{x_3} dx_3 = \int_0^x dx_1 \frac{x_1^2}{2} = \frac{x^3}{3 \cdot 2}

\displaystyle \int_0^x dx_1 \int_0^{x_1} dx_2 \int_0^{x_3} dx_3 \int_0^{x_4} dx_4 = \int_0^x dx_1 \frac{x_1^3}{3 \cdot 2} = \frac{x^4}{4 \cdot 3 \cdot 2}

Observando como opera esta progresión, no es difícil observar (o demostrar por inducción) que

\displaystyle \int_0^x dx_1 \int_0^{x_1} dx_2 \ldots \int_0^{x_{n-1}} dx_n = \frac{x^n}{n!}.

Si aplicamos esto a la expresión anterior de f(x) obtendremos

\displaystyle f(x) = f(0) + f'(0) x + f''(0) \frac{x^2}{2} + \ldots + \int_0^x dx_1 \ldots \int_0^{x_{n-1}} dx_n f^{(n)}(x_n),

donde f^{(n)} representa la derivada n-ésima de f.

Acotando la última integral, podemos ver que

\displaystyle \left| \int_0^{x_{n-1}} dx_n f^{(n)}(x_n) \right| < \int_0^{x_{n-1}} dx_n M = M x_{n-1},

siendo M una cota del valor absoluto de f^{(n)}(x) en el intervalo [0, x_{n-1}] (supongo x > 0 sin pérdida de generalidad). Eso implica que

\displaystyle \left| \int_0^x dx_1 \ldots \int_0^{x_{n-1}} dx_n f^{(n)}(x_n) \right| < M \int_0^x dx_1 \ldots \int_0^{x_{n-1}} dx_n,

donde M es ahora una cota del valor absoluto de f^{(n)}(x) en el intervalo [0, x]. Esta última integral podemos evaluarla explíctamente, dándonos una cota explícita para este último término en función de x y una cota para f^{(n)} en el intervalo [0, x]:

\displaystyle \left| \int_0^x dx_1 \ldots \int_0^{x_{n-1}} dx_n f^{(n)}(x_n) \right| < M \frac{x^n}{n!}.

Combinando todos estos resultados, y usando la nomenclatura f^{(0)} = f por consistencia, podemos decir que, si f es n veces derivable en un entorno de 0,

\displaystyle f(x) = \sum_{i=0}^{n-1} \frac{f^{(i)}(0)}{i!} + Err(x),

estando el término de error Err(x) acotado del siguiente modo:

\displaystyle |Err(x)| = \left| \int_0^x dx_1 \ldots \int_0^{x_{n-1}} dx_n f^{(n)}(x_n) \right| < M \frac{x^n}{n!}.

Esto significa que, si asumimos que la derivada n-ésima toma valores O(1), truncar la serie después del término con la derivada de orden n-1 nos dejará una aproximación O(x^n) a la función.

Próximo post

Con algo de suerte, será sobre la detección de cubos…

Permalink Dejar un comentario

Una respuesta y otro problema

Marzo 7, 2010 at 9:02 pm (C++, informática, matemática)

[Continuacion de Excusas et cetera]

Además de las soluciones de Demian y Guille, armé 5 soluciones propias con distintas características. Pueden verse todos los tests realizados a estas soluciones en el repositorio SVN “asociado” a este blog.

La solución de Guille fallaba a causa del problema identificado por Demian y puede verse a mi solución 4 como una adaptación del algoritmo recursivo de Guille para realizar el bubble sort completo. Aunque bubble sort suele ser duramente criticado, es de hecho óptimo (en las circunstancias apropiadas, ver Improved Bounds on Sorting by Length-Weighted Reversals, Lemma 59 :-D ).

Mis soluciones

A continuación muestro las distintas soluciones que planteé, cada una demostrando algo en particular.

Solución 1 – Arrays

Esta es una solución simple, esencialmente idéntica a la primera de Demian.

void f1(int* a, size_t n)
{
    size_t i = 0, j;
    for (j = 0; j < n; j++)
        if (a[j])
            a[i++] = a[j];
    for (; i < n; i++)
        a[i] = 0;
}
Solución 2 – Punteros

Esta es una adaptación simple de la Solución 1 para utilizar punteros.

void f2(int* p, size_t n)
{
    int *q, *r;
    for (q = p, r = p + n; q < r; q++)
        if (*q)
            *p++ = *q;
    for (; p < r; p++)
        *p = 0;
}
Solución 3 – Obfuscated

Esta es una versión optimizada para ser difícil de entender, empleando solo una definición y un “statement”.

void f3(int* p, size_t n)
{
    int *q = p, *r = p + n;
    while ((r - p) *
           (r - q ? *q ? (*p++ = *q++) : ++q - p : (*p++ = 0) + 1));
}
Solución 4 – Recursiva

Esta es la solución de Guille, adaptada para realizar una pasada completa de “burbujeo” en cada llamada recursiva. De las soluciones que encontré, es la única que no requiere variables adicionales (aunque obviamente utiliza múltiples instancias del stack frame).

void f4(int* p, size_t n)
{
    if (n == 0)
        return;
    f4(p + 1, n - 1);
    while (n >= 2)
    {
        if (!p[0] && p[1])
        {
            /* xor swap */
            /* Proof:
               p0' = p0 ^ p1
               p1f = p1 ^ p0'
                   = p1 ^ p0 ^ p1
                   = p0
               p0f = p0' ^ p1f
                   = p0 ^ p1 ^ p0
                   = p1
            */
            p[0] ^= p[1];
            p[1] ^= p[0];
            p[0] ^= p[1];
        }
        p++;
        n--;
    }
}
Solución 5 – Utilizando solo una variable extra

Esta solución es una “mezcla” de las soluciones 1 y 2 modificada para utilizar solo una variable temporal. Dudo que sea posible resolver el problema sin variables temporales extras (aparte de los parámetros), pero no parece simple probarlo…

void f8(int* a, size_t n)
{
    size_t j;
    for (j = 0; j < n; j++)
    {
        if (a[j])
            a[0] = a[j], a++, j--, n--;
    }
    for (j = 0; j < n; j++)
        a[j] = 0;
}

Un problema matemático

Hace unos días vi en el blog de Michael Nielsen un interesante problema:

Dado un conjunto de 51 números dentro del intervalo 1-100 (incluyendo los extremos), demostrar que existen al menos dos números en ese conjunto que no comparten ningún divisor (a excepción de 1, obviamente).

(No, no vi inmediatamente la respuesta :-D )

Permalink 2 comentarios

Excusas et cetera

Marzo 2, 2010 at 4:50 pm (C++, informática)

Originalmente tenía pensado publicar la primera parte de la respuesta al post sobre detección de cubos este fin de semana pasado, pero los videos de Susskind y Coleman sobre QFT se interpusieron :-D

Siempre creí que las clases expositivas eran bastante redundantes pero, cuando el material es dificultoso como en este caso, el inevitable “unpacking” que ocurre al explicar los temas “en vivo” se vuelve una ventaja importante. Los videos de Susskind son bastante fáciles de seguir para los que fueron siguiendo la serie en general, pero a veces pueden resultar un poco “handwavy”. Los de Coleman son bastante más rigurosos pero tienen mala calidad de filmación (incluso para 1975) y presuponen más conocimientos matemáticos y de mecánica cuántica de los que tengo, lo que me dificulta seguirlos.

Para no dejar todo en la nada, ofrezco un problema de programación C inspirado en uno visto (si no recuerdo mal) en algún comentario en el blog de Eric Raymond:

Enunciado:

Dado un vector de enteros, implementar una función que mueva “en el lugar” los ceros al final del vector dejando los otros enteros en el orden original. El prototipo de la función deberá ser:

void foo(int* a, size_t n).

Ejemplos de efectos sobre el vector pasado por parámetro:

{5, 0, 4, 0, 0, 3, 0, 0, 0, 2} --> {5, 4, 3, 2, 0, 0, 0, 0, 0, 0}
{0, 0, 0, 1, 2, 3, 0, 4} --> {1, 2, 3, 4, 0, 0, 0, 0}

Algunas respuestas a esto, salvando fuerza mayor, este fin de semana. Ustedes pueden dejar otras respuestas en los comentarios que serán debidamente analizadas… :-D

Permalink 10 comentarios

Detectando los cubos

Febrero 15, 2010 at 1:11 am (informática, matemática)

Para romper el largo lapso sin posts, presento el siguiente problema (robado de inspirado en uno visto en el blog de Petr Mitrichev):

Desarrollar un algoritmo que, dados 8 puntos en el plano (x_1, y_1), (x_2, y_2), \ldots, (x_8, y_8), determine si existen valores z_1, z_2, \ldots, z_8 tales que los puntos (x_1, y_1, z_1), (x_2, y_2, z_2), \ldots, (x_8, y_8, z_8) forme un cubo en el espacio. O, en otras palabras, que determine si los 8 puntos en el plano pueden formar la proyección paralela de un cubo.

Figura 1. 8 puntos en el plano que pueden verse como la proyección paralela de un cubo.

Figura 2. 8 puntos en el plano que no pueden verse como la proyección paralela de un cubo.

El algoritmo debe poder determinar que 8 puntos tales como (-0.343731, 0.404338), (-0.383761, 0.233634), (0.582329, -0.512121), (0.474738, 0.331386), (-0.236140, -0.439169), (0.622359, -0.341416), (-0.276170, -0.609873) y (0.514768, 0.502091) pueden verse como la proyección de un cubo (ver figura 1); por otro lado, tiene que clasificar al conjunto de puntos (0.085135, -0.025253), (-0.196956, 0.110847), (0.098387, -0.104055), (-0.091524, 0.268373), (0.051493, -0.111503), (-0.183703, 0.032046), (0.190567, 0.132272) y (-0.289135, -0.125480) como incapaz de ser la proyección de un cubo (ver figura 2).

En no menos de una semana, propondré una respuesta (no robada :-D ) al problema.

Permalink 2 comentarios

Next page »