# Expected value – Some “solutions”

(Follow-up to Expected value.)

This problem is known as the St. Petersburg paradox, and it was studied in the 18th century by Daniel Bernoulli and other mathematicians. In this post we will do a simple expected value analysis, but any kind of realistic analysis must consider the phenomenon of risk aversion.

#### Expected value

Let’s call Hi the event where a head appears at the i-th toss, and Hi the event where a tail appears at the same toss number. If we call Pn the probability that n tosses are required to get a head, we have:

$\displaystyle P_1 = P(H_1) = \frac{1}{2}$

$\displaystyle P_2 = P(T_1) P(H_2) = \frac{1}{2}\frac{1}{2} = \frac{1}{4}$

$\displaystyle P_3 = P(T_1) P(T_2) P(H_3) = \frac{1}{2}\frac{1}{2} \frac{1}{2} = \frac{1}{8}$

and, generalizing,

$\displaystyle P_n = P(T_1) P(T_2) \hdots P(T_{n-1}) P(H_n) = \left(\frac{1}{2}\right)^n = \frac{1}{2^n}$

As the bettor receives $2n from the house if n tosses are required, the expected balance is $\displaystyle E[H] = \sum_{n=1}^\infty P_n 2^n - 100 = \sum_{n=1}^\infty \frac{2^n}{2^n} - 100 = \sum_{n=1}^\infty 1 - 100 = \infty$. So it would seem a perfect bet to make: what can be better than a bet with “infinite expected value”? 😀 #### Simulation Let’s make a simulation to check if this result is reasonable: import random N = 10000 def get_tosses_until_head(): tosses = 0 while True: coin = random.randint(0, 1) tosses += 1 if coin == 0: break return tosses def get_balance(): balance = -100 tosses = get_tosses_until_head() balance += 2 ** tosses return balance def main(): random.seed(0) min_balance = max_balance = accum_balance = get_balance() for i in xrange(1, N): balance = get_balance() if balance < min_balance: min_balance = balance if balance > max_balance: max_balance = balance accum_balance += balance print 'Statistics after %d tries:' % N print ' Minimum balance: %d' % min_balance print ' Average balance: %d' % (float(accum_balance) / N) print ' Maximum balance: %d' % max_balance if __name__ == '__main__': main()  Running the simulation we get: Statistics after 10000 tries: Minimum balance: -98 Average balance: -75 Maximum balance: 32668  The average seems quite negative, how can we reconcile this with an “infinite expected value”? The answer is the same as in the case of the martingale: low probability events with huge returns skew the expected balance. Let’s see how the normal expected value estimator works with a random variable $X$, with finite expected value $\bar{x}$ and finite variance $\sigma^2_X$: $\displaystyle \widehat{x} = \frac{1}{N}\sum_{i=0}^N X_i$ $\displaystyle E[\widehat{x}] = E[\frac{1}{N}\sum_{i=0}^N X_i]$ $\displaystyle E[\widehat{x}] = \frac{1}{N}\sum_{i=0}^N E[X_i]$ (linearity of expectation) $\displaystyle E[\widehat{x}] = \frac{1}{N}\sum_{i=0}^N \bar{x}$ ($X_i$ are identically distributed) $\displaystyle E[\widehat{x}] = \bar{x}$ (the estimator is unbiased) $\displaystyle \sigma^2(\widehat{x}) = E[\widehat{x}^2 - E[\widehat{x}]^2]$ $\displaystyle \sigma^2(\widehat{x}) = E[\widehat{x}^2 - \bar{x}^2]$ ($\bar{x}$ is an unbiased estimator) $\displaystyle \sigma^2(\widehat{x}) = E[\widehat{x}^2] - \bar{x}^2$ (linearity of expectation; $\bar{x}$ is a number) Replacing $\widehat{x}$ by its definition and focusing in $E[\widehat{x}^2]$: $\displaystyle E[\widehat{x}^2] = E\left[\left(\frac{1}{N}\sum_{i=0}^N X_i\right)^2\right]$ $\displaystyle E[\widehat{x}^2] = \frac{1}{N^2}E\left[\left(\sum_{i=0}^N X_i\right)^2\right]$ (linearity of expectation) $\displaystyle E[\widehat{x}^2] = \frac{1}{N^2}E\left[\sum_{i=0}^N\sum_{j=0}^N X_i X_j\right]$ (square of a multinomial) $\displaystyle E[\widehat{x}^2] = \frac{1}{N^2}\sum_{i=0}^N\sum_{j=0}^N E[X_i X_j]$ (linearity of expectation) $\displaystyle N^2 E[\widehat{x}^2] = \sum_i E[X_i^2] + \sum_{\substack{i, j \\ j \ne i}} E[X_i] E[X_j]$ ($X_i$ is independent of $X_j$ if $i \ne j$) $\displaystyle N^2 E[\widehat{x}^2] = N E[X_i^2] + (N^2 - N) \bar{x}^2$ ($X_i$ are identically distributed) Integrating this result in the $\sigma^2(\widehat{x})$ expression: $\displaystyle \sigma^2(\widehat{x}) = \frac{N}{N^2} E[X_i^2] + \frac{N^2 - N}{N^2} \bar{x}^2 - \bar{x}^2$ $\displaystyle \sigma^2(\widehat{x}) = \frac{N}{N^2} E[X_i^2] - \frac{N}{N^2} \bar{x}^2$ $\displaystyle \sigma^2(\widehat{x}) = \frac{1}{N} \sigma^2_X$ This guarantees the convergence of the estimator to the expected value. But in our case there is no expected value (“it’s infinite”) and, consequently, we cannot even define the variance. Then the simulation is not very useful to get the expected value, as the estimator doesn’t need to converge to the expected value. #### House without unbounded wealth One (very realistic) situation where even the expected balance is negative is when the house has reasonably bounded wealth. For example, if we are betting against an agent whose possessions are merely equal to the entire world (!), the amounts we can be paid are bounded by ~$ 200·1012, giving a very different expected value:

$\displaystyle E[H] = \sum_{n=1}^\infty P_n \min(2^n, 200\cdot10^{12}) - 100$

$\displaystyle E[H] = \sum_{n=1}^{47} \frac{2^n}{2^n} + \sum_{n=48}^\infty \frac{200\cdot10^{12}}{2^n} - 100$

$\displaystyle E[H] = 47 + \frac{200\cdot10^{12}}{2^{47}}\sum_{n=1}^\infty \frac{1}{2^n} - 100$

$\displaystyle E[H] = 47 + \frac{200\cdot10^{12}}{2^{47}} - 100 < 47 + 2 - 100 = -51$

# Expected value

(In memory of Toxie 😀 )

Let’s suppose the following bet were proposed:

• The bettor pays $100 to the house. • An honest coin is tossed until a head appears. • The bettor receives$ 2n from the house, being n the number of coin tosses.

Is this bet worth making? Why?

An analysis will be given in the next post.

# Blue-eyed islanders: a solution

Follow up to Product & Sum.

As in the “three wise men” puzzle, it’s useful to begin by analyzing simplified variants of the problem. Lets start analyzing the variant where only one islander has blue eyes and the other 999 have brown eyes.

When the foreigner gives his message, as the blue-eyed islander can see that the other islanders are all brown-eyed, he now knows his own eye color and must commit suicide the following day at noon. But, when the blue-eyed islander kills himself, all the other islanders know that they are brown-eyed and they must commit suicide at noon the next day (we can assume that each islander knows that his eye color is either blue or brown).

Now we can examine the variant where two of the 1000 islanders are blue-eyed. In this case the foreigner’s message does not give any information immediately, as any inhabitant of the island knows by observation that there is at least one blue-eyed islander. But the absence of a suicide at noon the following day shows to each of the two blue-eyed islanders that the other one is not the only blue-eyed inhabitant of the island. Then, in the second noon after the message, the blue-eyed islanders commit suicide, followed by the rest of the tribe the next day.

This reasoning can be extended to the full case, where the 100 blue-eyed islanders will commit suicide at the 100th noon since the foreigner’s message and the 900 brown-eyed ones will do the same at the following noon (the 101st).

# Solving the product & sum riddle

Followup to Product & Sum.

Contradicting what was said in the previous post, the blue eyed islanders puzzle will be solved in a future post to avoid making this post excessively long. 😀

#### Product & Sum

One way to visualize the structure of this problem is to take each assertion by P and S as a filter to be applied over the initial set of integer pairs. This solution will then represent each of these filters as a block of Python code, with additional text explaining how each filter is connected with the associated assertion.

Given are X and Y, two integers, greater than 1, are not equal, and their sum is less than 100. S and P are two mathematicians; S is given the sum X+Y, and P is given the product X*Y of these numbers.

It is clear that we can restrict the first number to be strictly less than the second number, as the order of the numbers cannot be determined from the given data. Then we can get all the allowed pairs with the following Python code:

all_pairs = [(x, y)
for y in range(2, 100) for x in range(2, y)
if x + y < 100]


– P says “I cannot find these numbers.”

This implies that there are multiple allowed pairs whose products match the value that was given to P. We can start counting the number of pairs with each possible product:

num_pairs_by_prod = {}
for x, y in all_pairs:
if x * y not in num_pairs_by_prod:
num_pairs_by_prod[x * y] = 0
num_pairs_by_prod[x * y] += 1


The pairs allowed by P’s assertion are those whose product is shared with other pairs:

pairs_1 = [(x, y) for x, y in all_pairs if num_pairs_by_prod[x * y] > 1]


– S says “I was sure that you could not find them.”

Then we know that the value of the sum is enough to know that the product cannot identify the pair of integers. The set of sums of integer pairs that can be identified by their products is:

identif_by_prod_pairs_sums = set(x + y for x, y in all_pairs
if num_pairs_by_prod[x * y] == 1)


As the sum is enough to know that the integer pair cannot be identified by its product, the sum of the pair cannot be in the above set:

pairs_2 = [(x, y) for x, y in pairs_1
if x + y not in identif_by_prod_pairs_sums]


– P says “Then, I found these numbers.”

This indicates that in pairs_2, the list of pairs restricted by the first two assertions, the correct pair can be identified by its product. Then we can do essentially the same we did in the first step but now keeping the pairs identifiable by their products:

num_pairs_by_prod_2 = {}
for x, y in pairs_2:
if x * y not in num_pairs_by_prod_2:
num_pairs_by_prod_2[x * y] = 0
num_pairs_by_prod_2[x * y] += 1
pairs_3 = [(x, y) for x, y in pairs_2 if num_pairs_by_prod_2[x * y] == 1]


– S says “If you could find them, then I also found them.”

This implies that the pair can be identified by its sum from the list restricted by the first three assertions. We can get the final pairs list repeating the previous step, but now searching for a unique sum:

num_pairs_by_sum_3 = {}
for x, y in pairs_3:
if x + y not in num_pairs_by_sum_3:
num_pairs_by_sum_3[x + y] = 0
num_pairs_by_sum_3[x + y] += 1
pairs_4 = [(x, y) for x, y in pairs_3 if num_pairs_by_sum_3[x + y] == 1]


Putting the code together, adding a print statement to get the final list of pairs and running the code we get

[(4, 13)]


matching the results in the literature.

# Product & Sum

Follow up to Three wise men.

#### Solution to “three wise men”

A certain king wishes to determine which of his three wise men is the wisest. He arranges them in a circle so that they can see and hear each other and tells them that he will put a white or black spot on each of their foreheads but that at least one spot will be white. In fact all three spots are white. He then offers his favor to the one who will first tell him the color of his spot. After a while, the wisest announces that his spot is white. How does he know?

One way of visualizing the solution is to think hypothetical scenarios where some of the men have black spots in their foreheads. Lets start supposing that two of the men have black spots: then it will be obvious to the third that he has a white spot, as at least one of the three must have one.

Now suppose that only one of the men has a black spot and take the point of view of the “wisest” (quicker at least) man with a white spot. He cannot directly infer the color of his spot, as he sees one spot of each color. But, as we have seen in the last paragraph, the man in front of him with a white spot would have known the color of his spot if he (the observer) would have had a black spot in his forehead. As he didn’t speak, then he knows that the color of his spot must be white.

Finally, lets attack the full problem taking the place of the wisest man. He knows that, if he would have had a black spot in his forehead, one of the other men would have spoken, by the previously mentioned reasons. As they haven’t done that, he can conclude that the spot in his forehead must be white.

#### Blue eyed islanders puzzle

A similar, but more difficult, problem is the blue eyed islanders puzzle (this version was stolen from Terry Tao):

There is an island upon which a tribe resides. The tribe consists of 1000 people, 100 of which are blue-eyed and 900 of which are brown-eyed. Yet, their religion forbids them to know their own eye color, or even to discuss the topic; thus, one resident can see the eye colors of all other residents but has no way of discovering his own (there are no reflective surfaces). If a tribesperson does discover his or her own eye color, then their religion compels them to commit ritual suicide at noon the following day in the village square for all to witness. All the tribespeople are highly logical, highly devout, and they all know that each other is also highly logical and highly devout. One day, a blue-eyed foreigner visits to the island and wins the complete trust of the tribe. One evening, he addresses the entire tribe to thank them for their hospitality. However, not knowing the customs, the foreigner makes the mistake of mentioning eye color in his address, mentioning in his address “how unusual it is to see another blue-eyed person like myself in this region of the world”. What effect, if anything, does this faux pas have on the tribe?

#### Product & sum

This is another classical problem, first proposed (apparently) by Hans Freudenthal en 1959, and also called the “Impossible” Puzzle:

Given are X and Y, two integers, greater than 1, are not equal, and their sum is less than 100. S and P are two mathematicians; S is given the sum X+Y, and P is given the product X*Y of these numbers.
P says “I cannot find these numbers.”
S says “I was sure that you could not find them.”
P says “Then, I found these numbers.”
S says “If you could find them, then I also found them.”
What are these numbers?

Hint: a computer is very useful to solve this problem.

Solution for both puzzles in the next post.

# Three wise men

An interesting traditional puzzle (this formulation was stolen from John McCarthy):

A certain king wishes to determine which of his three wise men is the wisest. He arranges them in a circle so that they can see and hear each other and tells them that he will put a white or black spot on each of their foreheads but that at least one spot will be white. In fact all three spots are white. He then offers his favor to the one who will first tell him the color of his spot. After a while, the wisest announces that his spot is white. How does he know?

# 400 teracycles, 200 gigabytes, 7 collisions

(Follow-up to Sumando subconjuntos – La respuesta.)

#### Introduction

In a previous post (and in a comment from Guille [ɡiˈʒe] 😀 ) we have seen how the pigeonhole principle implies that a set of 70 numbers in the range [1018, 1019) must have two subsets with equal sum. But this is a non-constructive proof, as it doesn’t give us the two subsets with the same sum. To rectify this omission, in this post we will see how this “sum collision” can be obtained.

#### Finding collisions: simple cases

We can start by defining the problem in a more general way: given a sequence of elements xi and a function f(), find two elements of the sequence, xj and xk, such that f(xj) = f(xk). In this way the sum collision problem can be reduced to finding a duplicate in the associated sequence of sums, f(xi).

Two common ways to get duplicates in a sequence are the following:

def find_duplicate_sort(g):
sl = sorted(g)
prev = None
for e in sl:
if e == prev:
return e
prev = e
return None

def find_duplicate_set(g):
es = set()
for e in g:
if e in es:
return e
return None


The first one has O(n log n) complexity and the second one has O(1) complexity if we use a hash-based set. As the set-based approach also has a lower constant, we will use this approach in the rest of this post.

This algorithm works well if the sequences to be inspected for duplicates can be fit entirely in RAM, but in this case we have seen that tens of billions of elements must be inspected to have a high probability of finding a collision. In the next section we will analyse how this restriction can be evaded.

#### Finding collisions in big sets

Each of the subsets to be evaluated in this problem can be encoded using 70 bits and, to allow a simple and fast sum calculation algorithm to be used, this was rounded up to 10 bytes. Then, if we want to inspect 20 billion subsets of 35 elements to get a very high probability of not wasting the calculation time, we will need 200 GB to store the whole sequence. 200 GB of data cannot be stored in the RAM of an usual machine, but it’s quite easy to store this amount of data in a hard disk nowadays.

To allow a fast hash-set based RAM collision search while keeping the bulk of the data in disk, we can take the following approach:

1. Generate in RAM a big number of random subsets and sort them by their sums.
2. Find a vector of numbers (to be called “pivots”) aplitting the sorted subsets vectors by their sums in approximately equal-sized segments. (The segment n will be composed by all the subsets whose sum is between pivot n-1 and pivot n).
3. Generate in RAM a big number of random subsets and sort them by their sums.
4. Split the sorted subset vector in segments using the previously generated pivots and append each segment to an associated segment file (for example, append segment 3 to 0003.bin).
5. Repeat steps 3 and 4 until getting enough subsets.
6. Check each segment file at a time for collisions.

If we choose enough pivots, the size of each segment file will be small enough to allow easily doing step 6 with a hash-based set (each segment file won’t have the same size, as the generated subsets are random; but the law of large numbers ensures that their sizes won’t be very different).

#### Source code & parameters

The (ugly) C code that checked for collisions can be found in the repository associated with this blog. The chosen parameters were:

• Number of subsets: 2·1010.
• Number of subsets in RAM: 107.
• Number of elements in each subset: 35 (constant).
• Number of segments: 1000.
• Number of slots in the hash-set: 227.

#### Results

The first stage (segment file generation) elapsed time was approximately 41 hours, somewhat over my original estimation of 36 hours, and the segment file range ranged from 194827630 bytes to 206242980 bytes. The second stage (collision detection inside each segment file) lasted for 12-18 hours.

The output of the second stage (discarding files where no collisions were found) was:

Processing file 218...  Collision between identical elements.
DONE in 40.754850s.
Processing file 363...  Collision between different elements!!!
ed940f4f5710c6351a00
a35377a5a74a03961c00
DONE in 38.585990s.
Processing file 394...  Collision between different elements!!!
9ab5dfa6312e090e2a00
6558c9c8eab667233400
DONE in 35.570039s.
Processing file 409...  Collision between different elements!!!
2429d70f6ff500ab2c00
423cf8c758740ac73b00
DONE in 34.499926s.
Processing file 434...  Collision between different elements!!!
5a8bf5a532915f213100
496ebc281e9958713e00
DONE in 32.610608s.
Processing file 475...  Collision between different elements!!!
c35d3f427a7c488c0e00
4e37132b717a287a1900
DONE in 21.971667s.
Processing file 655...  Collision between different elements!!!
4653023676ea8e5f1900
6abb914e641ba45a3900
DONE in 21.514123s.
Processing file 792...  Collision between different elements!!!
8a4a63d85de377130600
853d3d45b15e0c471e00
DONE in 21.506716s.


Each set is represented as a byte string with bit number increasing when advancing through the byte string. For example, ed940f4f5710c6351a00 represents the bit string 10110111001010011111000011110010111010100000100001100011101011000101100000000000 and, consequently, the subset with indices 0, 2, 3, 5, 6, 7, 10, 12, 15, 16, 17, 18, 19, 24, 25, 26, 27, 30, 32, 33, 34, 36, 38, 44, 49, 50, 54, 55, 56, 58, 60, 61, 65, 67, 68. Its elements are

5213588008354709077 9115813602317594993
1796745334760975384 3579709154995762145
2312952310307873817 3627590721354606439
5763055694478716846 2730952213106576953
4868653895375087301 9737387704190733086
9262565481673300485 5968266991171521006
6752113930489992229 3772194655144630358
9029836747496103755 3318990262990089104
9205546877037990475 9849598364470384044
1376783969573792128 1108556560089425769
7820574110347009988 6951628222196921724
4776591302180789869 7999940522596325715
2290598705560799669 7835010686462271800
8998470433081591390 9131522681668251869
9096632376298092495 5295758362772474604
5953431042043343946 3151838989804308537
8643312627450063997 3624820335477016277
3782849199070331143


and its sum is 203743882076389458417.

In the same way, the set a35377a5a74a03961c00 has elements

5213588008354709077 9011219417469017946
3579709154995762145 3627590721354606439
5941472576423725122 4317696052329054505
2730952213106576953 5014371167157923471
9737387704190733086 9262565481673300485
5968266991171521006 5917882775011418152
5866436908055159779 9233099989955257688
3772194655144630358 3318990262990089104
9990105869964955299 2664344232060524242
1376783969573792128 1108556560089425769
7820574110347009988 9889945707884382295
7652184907611215542 8082643935949870308
4271233363607031147 6415171202616583365
6393762352694839041 2290598705560799669
7481066850797885510 5295758362772474604
5953431042043343946 9929081270845451034
7207546018039041794 3624820335477016277
3782849199070331143


and 203743882076389458417 as its sum, the same value as the previous different subset. 😀

#### JS puzzle

The previous post asked what happened when the JS code

(function(f){f(f)})(function(g){g(g)})


was executed. In first place we can observe that the two anonymous functions are really the same, as they only differ in the name of a dummy parameter. Let’s call this function u() and observe what happens when we apply u() over a function f():

u(f) -> f(f)


It’s clear that u() applies its argument over itself. As in this case we are applying u over itself, the result will be

u(u) -> u(u)


giving us an infinite recursion.

# Sumando subconjuntos – La respuesta

#### Solución “no constructiva”

(Esencialmente es la misma que Guillermo presentó en un comentario del post anterior, solo difiere en la forma de acotar 🙂 )

La solución es mucho más simple de lo que podría parecer a primera vista y se basa en el el crecimiento exponencial de la cantidad de subconjuntos respecto al tamaño del conjunto “base”.

Consideremos los posibles valores que puede tomar la suma de un subconjunto de los números: como todos los elementos son positivos, la mínima suma es 0 y corresponde al subconjunto vacío, mientras que la máxima corresponde al conjunto completo. Al tener todos los elementos 19 dígitos, si llamamos S a la suma de los elementos de un subconjunto cualquiera, tendremos:

0 ≤ S < 70·1019 < 1021.

Aprovechando la relación 103 < 210, vemos que

S < 1021 < 270.

Pero, como la cantidad de subconjuntos de un conjunto de 70 elementos es exactamente 270, acabamos de probar que hay más subconjuntos que posibles valores de la suma de sus elementos. Por lo tanto es inevitable que haya al menos dos subconjuntos con el mismo valor de la suma de sus elementos.

Un punto interesante es que la solución es la misma si agregamos el requerimiento de que los conjuntos sean disjuntos. Para ver esto, supongamos que tenemos dos subconjuntos distintos A y B con la misma suma:

sum(A) = sum(B).

Si analizamos los conjuntos A\B y B\A, observaremos que

sum(A\B) = sum(A) – sum(AB)

sum(B\A) = sum(B) – sum(AB)

sum(A\B) = sum(B\A).

Como obviamente A\B y B\A no pueden compartir elementos, entonces cada par de subconjuntos distintos con igual suma tiene asociado un par de subconjuntos disjuntos con igual suma.

#### Una solución constructiva?

Este problema es conocido como NP-completo, por lo que no cabe esperar una solución eficiente y general. Pero podría ser que con estos valores de parámetros fuera factible encontrar una solución particular por lo que realizaremos un análisis probabilístico.

Consideremos una distribución aleatoria uniforme en el intervalo (0, 1019): claramente tendrá una densidad de probabilidad constante e igual a 10-19. Esto implica que si elegimos un número aleatoriamente de esa distribución, considerándola ahora discreta, la probabilidad sería de 10-19, lo que concuerda con nuestra intuición.

Si analizamos ahora la distribución de probabilidad de la suma de 35 variables elegidas uniformemente de dicho intervalo, tendremos una distribución de probabilidad conocida como Irwin Hall. Pero para nuestros fines podemos aproximarla por una normal con media igual a 35 veces el promedio de los números, 2.08·1020, y desviación estándar igual a la desviación estándar de los números multiplicada por la raíz cuadrada de 35, 1.57·1019. Comparando la distribución normal con un histograma obtenido mediante simulación,

Comparación de un histograma de las sumas de combinaciones de 35 números con la aproximación normal.

se observa que, aunque la coincidencia de la media es muy buena, la desviación estándar de la normal es significativamente superior a la que se observa en el histograma. Esto se debe a que, a diferencia de lo que se supone en una distribución Irwin Hall, las variables sumadas no son independientes (podemos verlo más claramente pensando en el caso de 70 variables: la suma es única, pero la distribución de Irwin Hall nos daría una desviación estándar aún mayor que en el caso de 35 variables). Como este error nos lleva a sobrestimar la dificultad de encontrar una colisión, podemos ignorar esta discrepancia por el momento, ya que nos llevará a una aproximación conservadora.

Lo que necesitaríamos ahora es calcular la probabilidad de encontrar dos combinaciones con igual suma. Aunque queda claro que las probabilidades son mayores que en el caso de que las sumas estuvieran distribuidas uniformemente en el intervalo (0, 70·1019), no es tan simple calcular un valor exacto para esta probabilidad. Pero podemos hacer una estimación suponiendo a los valores distribuidos uniformemente en el intervalo ±1σ, lo que implica descartar menos del 35% de los valores.

Aplicando la expresión derivada en el “problema del cumpleaños” para una probabilidad de colisión p:

$\displaystyle n(p) \approx \sqrt{2 \cdot 1.57 \cdot 10^{19}\ln\left(\frac{1}{1-p}\right)}$.

Aceptando un 95% de probabilidad de éxito en la búsqueda de una colisión:

$\displaystyle n(0.95) \approx \sqrt{2 \cdot 1.57 \cdot 10^{19}\ln 20}$.

$\displaystyle n(0.95) \approx 9.7\cdot 10^{9}$.

Tomando en cuenta que se está desechando el 30% de los puntos, podemos concluir que examinar unos 2·1010 puntos nos debería dar una probabilidad elevada de encontrar una colisión.

Si bien es bastante complejo buscar una colisión entre más de 1010 valores, no es infactible. Pero sería mejor realizar algunas pruebas preliminares para ver si el resultado obtenido es razonable.

#### Prueba a escala

A modo de verificación podemos realizar una simulación a escala, empleando la suma de 30 números de 7 dígitos para disminuir la carga computacional. En este caso sigue estando garantizada la existencia de una colisión, ya que

0 ≤ S < 30·107 < 109 < 230.

Aplicando la aproximación normal, la suma de 15 números distribuidos seguirá una distribución de media 7.5·107 y varianza 15·1014/12. Por lo tanto la cantidad de sumas ensayadas para tener un 95% de probabilidad de encontrar una colisión (usando solo sumas dentro de 1σ) será:

$\displaystyle n(0.95) \approx \sqrt{2 \cdot 1.1 \cdot 10^7\ln\left(\frac{1}{1-0.95}\right)} \approx 8120$.

O sea que la probabilidad de requerir más de unos 14000 ( = 8200 / 0.6) intentos para encontrar una colisión debería ser pequeña. Y estos valores son mucho más simples de comprobar experimentalmente que los originales, ya que puede hacerse mediante un pequeño script:


import math, random

N = 1000

NUMBERS = [4401501, 1984319, 1064736, 6495167,
6402639, 7578056, 9301342, 6042069,
1144375, 5680006, 7450344, 9099174,
2011586, 8039920, 3010493, 1658370,
5927190, 3880633, 1318068, 9594698,
2877906, 1792394, 9120086, 7372216,
2141103, 7256943, 6603598, 8565018,
1698346, 3004879]

def run_until_collision():
combs_dict = {}
i = 0
while True:
comb = random.sample(NUMBERS, 15)
s = sum(comb)
if s in combs_dict:
return i, comb, combs_dict[s]
combs_dict[s] = comb
i += 1

def show_stats():
total = 0.0
sq_total = 0.0
print 'Please wait until the count reaches 100...'
for i in range(N):
tries, dum, dum = run_until_collision()
total += tries
sq_total += tries ** 2
if i % math.ceil(N / 100.0) == 0:
print '%d' % (i * 100.0 / N),
print '\n\nTries statistics'
print '  Average: %f' % (total / N)
print '  SD: %f' % math.sqrt(sq_total / N -
(total / N) ** 2)

if __name__ == '__main__':
show_stats()


La salida obtenida al ejecutar este script (descartando indicaciones de progreso) es:

Tries statistics
Average: 6123.049000
SD: 3360.260149


#### Conclusiones

Para poder encontrar una colisión en forma explícita en el problema original deberá ser necesario buscarla entre aproximadamente 1010 valores. Como cada subconjunto con su suma ocupará en el orden de 16 bytes y una PC suele tener menos que 100 GB de RAM, sería práctico buscar un método que permita buscar colisiones sin tener todos los valores previamente analizados en memoria (otra opción sería usar datos en disco, pero es menos elegante 😀 ).

En un próximo post analizaremos como implementar una solución a este problema y veremos si pueden encontrarse resultados.

# Sumando subconjuntos – Un problema

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.

# Ante la duda, aplicar Fourier

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

$\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

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