segunda-feira, 16 de julho de 2018

As Figurinhas dos Vingadores são Honestas? (parte 1)

O Brasil é um país muito estranho. Uma das inúmeras bizarrices brasileiras é que temos tabus temporários. Por exemplo, homem se vestir de mulher é tabu (vai acabar com a família, fim dos tempos, etc), mas no carnaval pode! Da mesma maneira, adulto comprando álbum de figurinhas é tabu (precisa crescer, largar os brinquedos, etc.), mas em época de Copa pode!


Sempre que chega época de Copa eu fico espantado com a dedicação que meus amigos tem para completar o álbum com os jogadores. E pensando um pouco no fenômeno, me veio a dúvida: quanto dinheiro você precisa para terminar um álbum desses?

Resolvi fazer um experimento para tentar calcular esse valor. Mas, ao invés de usar o álbum da Copa, escolhi completar o álbum dos Vingadores, já que é um universo que eu tenho mais afinidade. Os valores originais são os seguintes:

  • O álbum em si custa R$ 6.90, e você precisa de 180 figurinhas para completar. 
  • Cada pacotinho com 4 figurinhas custa R$ 1.20.
  • Existe um kit com 15 pacotinhos que você pode comprar com desconto por R$ 13.90. Mas é difícil de achar: eu só encontrei na loja Geek na av. Paulista.
  • Por fim, você sempre pode pedir figurinhas avulsas no site da editora (nesse caso você pode comprar até 40 figurinhas com preço individual de R$ 0.30, mais o frete de R$ 10.00).

Se você for uma pessoa sem imaginação, de cara já tem uma estratégia fácil para completar o álbum: é só pedir tudo pelo correio. Nesse caso, basta comprar o álbum e mais 180 figurinhas em 5 fretes, num total de R$110.90.

Tem como gastar menos que isso? Para descobrir, temos que analisar as figurinhas!

A Distribuição das Figurinhas


Para descobrir a distribuição das figurinhas, eu comprei ##t=106## pacotinhos (ou seja, ##n=106\times 4=424## figurinhas), e anotei o conteúdo de cada um deles (os dados originais estão no google docs). A primeira coisa que fiz foi traçar o histograma para ver o jeitão dos dados:

Pelo histograma dá para ver que faltam nove figurinhas para completar meu álbum (são os nove buracos no histograma). Por outro lado tem três figurinhas que repetiram 6 vezes, o que indica que possivelmente esse álbum não é honesto, já que algumas figurinhas são mais fáceis que outras.

Ou... talvez não? A intuição humana é notoriamente falha em intuir probabilidades, de repente essa impressão de que algumas figurinhas são mais fáceis é falsa. O jeito correto de decidir isso é fazendo um teste de hipótese. O teste em si é simples, mas eu vou fazer em câmera lenta caso alguém nunca tenha visto.

Para começar, precisamos de uma hipótese nula ##H_0## e um nível de significância ##\alpha##.

Nossa hipótese nula ##H_0## é que as figurinhas tem distribuição uniforme e são independentes. O teste vai nos falar se esse hipótese precisa ser rejeitada. O nível de significância nós podemos escolher, mas escolher o ##\alpha## é uma arte: você precisa levar em conta o número de amostras, e as taxas aceitáveis de falsos positivos e falsos negativos. Usualmente, ##\alpha##=5% funciona bem.

Qual o significado desses parâmetros? Fica simples de entender com um exemplo: suponha que todas as 424 figurinhas que eu comprei vieram repetidas. Em uma distribuição é uniforme, isso não é impossível, mas é extremamente improvável: a chance é de 1 em ##180^{424}\sim 1.7\times 10^{956}##. Esse valor é maior que o número de átomos no universo (muito maior!), então a chance de ter acontecido ao acaso é quase nula.

Se eu tirasse 424 figurinhas repetidas, eu teria bastante convicção de que as figurinhas não são uniformes. Quando fazemos um teste de hipótese, é essa ideia que estamos tentando capturar: os dados que coletamos são tão raros que não poderiam ter acontecido na prática? A chance de terem acontecido ao acaso é menor que ##\alpha##=5%? Se for, então talvez a hipótese ##H_0## seja falsa.

Existem vários testes de hipótese possíveis. Como queremos checar se nossas amostras seguem uma distribuição dada (uniforme), então o teste mais natural a se fazer é o teste do ##\chi^2## (pronuncia-se quiquadrado), que serve justamente para verificar se duas distribuições são iguais.

Você começa o teste do ##\chi^2## separando suas amostras em categorias. No nosso caso, o mais natural seria fazer cada figurinha ser uma categoria, mas isso não vai funcionar. O motivo é que o teste do ##\chi^2##, estritamente falando, só funciona com infinitas amostras. Como é impossível coletar infinitas amostras, nós assumimos que o resultado não muda muito de "infinitas amostras" para "muitas amostras".

Mas quantas amostras são "muitas amostras"? Em geral, a regra usada é que ##np>5##, onde ##n## é o número de amostras e ##p## é a menor probabilidade de uma amostra estar em uma categoria. Vamos fazer as contas: se cada figurinha é uma amostra, e as 180 figurinhas são uniformes, então ##p=1/180##. Como tenho ##n=424## figurinhas, então ##np=424/180\sim 2.35##. Oops, 2.35 é menor que 5, então 180 categorias é muito.

Vamos reduzir então. Eu vou escolher ##c=5## categorias, porque 180 divide certinho por 5. Cada categoria vai ter ##180/5=36## figurinhas, logo ##p=1/36## e ##np=424/36\sim 11.77##, agora sim! O histograma com cinco categorias é assim:


Agora nós vamos tabelar nossas amostras. A linha Observado mostra quantas amostras nós coletamos em cada categoria, a linha Esperado mostra quantas amostras uma distribuição ideal teria (no caso uniforme, todas são iguais a ##n/c=424/5=84.8##):

Categoria 1Categoria 2Categoria 3Categoria 4Categoria 5
Observado7590809584
Esperado84.884.884.884.884.8

Podemos calcular o ##\chi^2## agora. A intuição do teste é simples, vamos calcular qual o quadrado da diferença entre o observado e o esperado, e normalizar. Se esse número for muito grande, então tem algo errado. Calculando pela fórmula, temos:

$$\chi^2=\sum\frac{(Y_i-E_i)^2}{E_i}\sim 2.957$$
Esse valor, 2.957 é grande ou é pequeno? Para decidir isso precisamos primeiro calcular quantos graus de liberdade tem esse modelo. O número de graus de liberdade é o número de parâmetros que estamos estimando. Como nossa tabela tem 5 categorias, então o número de graus de liberdade é 4.

Opa! Se são 5 categorias, não tinha que ser 5 graus de liberdade? É verdade que estamos tentando estimar ##Y_1##, ##Y_2##, ##Y_3##, ##Y_4## e ##Y_5##. Mas nós sabemos que a soma de todos eles é ##n##, então se tivermos os quatro primeiros, o quinto é só calcular!

Em tempos idos esse era o momento em que você precisava sacar do bolso uma tabela de ##\chi^2## (como curiosidade, coloquei online a tabela que eu usava na Poli). Mas hoje em dia você pode ir direto no Mathematica para descobrir o valor de ##\chi^2## com ##\alpha##=5% e quatro graus de liberdade:
InverseCDF[ChiSquareDistribution[4], 1 - 0.05]
>> 9.48773

Como 2.957 é menor que 9.488, então não podemos rejeitar a hipótese de que as figurinhas são uniformes e independentes, e concluímos que sim, as figurinhas provavelmente são honestas.

Alguns de vocês irão chiar, eu sei. "Mas veja, algumas figurinhas saíram repetidas seis vezes, outras não apareceram nenhuma vez! Não tem como isso ser uniforme, claramente tem figurinhas que são mais fáceis!".

Bem, se você ainda estiver na dúvida, você pode comprar mais pacotinhos para deixar o teste mais robusto. Ao invés disso, eu vou usar só as figurinhas que eu tenho para ir direto na raiz da dúvida: afinal de contas, é normal tirar muitas figurinhas repetidas?

A Análise das Repetidas


A intuição nos diz que, em uma distribuição uniforme, todas as figurinhas tem que aparecer mais ou menos com a mesma frequência. Só que esse é mais um caso onde a intuição é falha! Quem lê o blog há mais tempo deve lembrar de um post sobre o paradoxo do aniversário, que diz exatamente isso: mesmo em distribuições uniformes, a chance de colisões é muito alta.

Podemos usar esse fato para tentar um novo teste de hipótese. Dessa vez, vamos contar quantas repetidas nós temos para cada figurinha, e separá-las por número de repetições. O histograma da quantidade de figurinhas repetidas fica assim:


O histograma foi traçado da quantidade de figurinhas repetidas:

Quantidade de repetições0123456
Figurinhas94150482453

O álbum tem 180 figurinhas, então a soma da segunda linha vale 180. Do total de figurinhas, nove delas não apareceram, três apareceram 6 vezes, e a maioria apareceu entre 1 e 3 vezes.

O formato do histograma parece familiar né? Se nós soubéssemos qual a é a fórmula exata para o número de figurinhas repetidas em uma distribuição uniforme, poderíamos usar o teste do ##\chi^2## para ver se a hipótese ainda se mantém. Só que agora não tem jeito, para descobrir a fórmula precisamos fazer as contas. Pule a caixa azul se você tem medo de combinatória:

O truque para usar combinatória analítica numerada (labelled) é sempre tentar achar um isomorfismo entre o seu problema e uma lista de inteiros. Para o nosso caso, imagine o seguinte experimento: temos ##n## bolinhas numeradas de ##1## a ##n##, e temos também ##m## caixinhas, cada caixa representando uma figurinha. A bolinha representa a ordem em que eu comprei a figurinha. Se eu atirar as bolinhas em caixas aleatórias, tenho um problema isomórfico a comprar figurinhas aleatórias.

Por exemplo: se a bolinha 5 cair na caixa 32, então significa que a quinta figurinha que eu comprei foi a de número 32. Se no fim do experimento a caixa 32 tiver seis bolinhas, então eu tenho a figurinha 32 repetida seis vezes.

No nosso experimento, as caixinhas são diferentes umas das outras, então vamos usar o operador ##\text{SEQ}_m(z)## para modelar as caixinhas. A ordem das bolinhas dentro da caixinha não importa (tanto faz se dissermos que a caixa 32 tem as bolinhas 3, 5, 9, 10 e 15; ou se tiver as bolinhas 9, 5, 3, 15, e 10, é a mesma coisa). Então dentro das caixinhas, usaremos o operador ##\text{SET}(z)##. O experimento completo é ##\text{SEQ}_m(\text{SET}(z))##.

Digamos que estamos interessados em saber quantas caixinhas vão terminar com ##s## bolinhas. Um jeito de medir isso é colocar um marcador nas caixinhas com ##s## elementos. Para isso, nós tiramos as das caixinhas os conjuntos de ##s## elementos, e colocamos eles de volta com o marcador, digamos que o marcador seja ##u##. Então o experimento para medir ##s## repetidas é ##\text{SEQ}_m(\text{SET}(z)-\text{SET}_s(z)+u \text{SET}_s(z))##.

Os experimentos vão mapear para funções geradoras exponenciais de duas variáveis. Das tabelas de combinatória analítica:
$$\begin{align*}
\text{SEQ}_m(z) &= z^m \\
\text{SET}(z) &= e^z \\
\text{SET}_s(z) &= \frac{z^s}{s!}\\
\text{SEQ}_m(\text{SET}(z)-\text{SET}_s(z)+u \text{SET}_s(z)) &=
\left(e^z-\frac{z^s}{s!}+u \frac{z^s}{s!}\right)^m \\
&= \left(e^z+ (u-1) \frac{z^s}{s!}\right)^m
\end{align*}$$
Se abrirmos essa função geradora em série, então o coeficiente do termo ##z^n u^k## indica quantas combinações possíveis de ##n## bolinhas vão ter ##k## caixinhas com ##s## bolinhas repetidas (como a função geradora é exponencial, você precisa multiplicar o coeficiente por ##n!## para ter o valor real). Se chamarmos esse coeficiente de ##R_s(n, k)##, então temos:
$$\sum_{n}\sum_{k}R_s(n,k) z^n u^k = \left(e^z+ (u-1) \frac{z^s}{s!}\right)^m $$
Vamos listar com cuidado o que temos. ##n! R_s(n,k)## é número de combinações possíveis com ##k## caixinhas contendo exatamente ##s## bolinhas. Logo, de todas as combinações possíveis de ##n## bolinhas, existem ##\sum k n! R_s(n,k)## caixinhas com ##s## bolinhas.

Por outro lado, o número total de combinações de ##n## bolinhas em ##m## caixinhas é ##m^n##. Como cada combinação dessas tem ##m## caixinhas, então o número total de caixinhas que podem ser preenchidas em todas as combinações é ##m^{n+1}##.

Juntando as duas contagens, a probabilidade de um caixinha tem ##s## bolinhas em ##n## arremessos é:
$$p_s(n)=\sum_{k>0}\frac{k \;n!}{m^{n+1}}R_s(n, k)$$
Mas olha só! Nós podemos conseguir essa somatória derivando a função geradora em relação a ##u##, e depois setando ##u=1##.
$$\left.\frac{\partial}{\partial u}\sum_{n}\sum_{k}R_s(n,k) z^n u^k \right|_{u=1}=\sum_n \sum_{k}k R_s(n,k) z^n$$
Agora é só fazer as contas:
$$\begin{align*}p_s(n)
&=\sum_k \frac{n!}{m^{n+1}}k R_s(n, k) \\
&=\frac{n!}{m^{n+1}} [z^n] \left.\frac{\partial}{\partial u}\sum_{n}\sum_{k}R_s(n,k) z^n u^k \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n] \left.\frac{\partial}{\partial u}\left(e^z+ (u-1) \frac{z^s}{s!}\right)^m  \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n] \left. m\left(e^z+ (u-1) \frac{z^s}{s!}\right)^{m-1} \frac{z^s}{s!} \right|_{u=1}\\
&=\frac{n!}{m^{n+1}} [z^n]  m\left(e^z\right)^{m-1} \frac{z^s}{s!} \\
&=\frac{n!}{s!\;m^n} [z^n]z^s e^{z(m-1)}\\
\end{align*}$$
Abrindo a exponencial em série:
$$\begin{align*}p_s(n)
&=\frac{n!}{s!\;m^{n}} [z^n]z^s\sum_{k\ge 0}\frac{(z(m-1))^k}{k!} \\ &=\frac{n!}{s!\;m^{n}} [z^n]\sum_{k\ge 0}\frac{z^{k+s}(m-1)^k}{k!} \\ &=\frac{n!}{s!\;m^{n}} [z^n]\sum_{p\ge s}\frac{z^p(m-1)^{p-s}}{(p-s)!} \\ &=\frac{n!}{s!\;m^{n}}\; \frac{(m-1)^{n-s}}{(n-s)!} \\ &=\frac{n!}{s!(n-s)!}\; \frac{(m-1)^{n-s}}{m^{n}} \\ &={{n}\choose{s}}\frac{(m-1)^{n-s}}{m^{n}} \\ \end{align*}$$
Essa é a solução exata, em forma fechada, que queríamos! Mas tem uma massagem que podemos fazer que vai ajudar a entender por que essa fórmula é assim. Em geral ##n## é bem maior que ##m##, então vamos reescrever ##m## de modo que ##n=\lambda m##, e calcular o limite quando ##n## é muito grande: $$\begin{align*}\lim_{n\to\infty}p_s(n) &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \frac{(n/\lambda-1)^{n-s}}{(n/\lambda)^{n}} \\ &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \left(\frac{n-\lambda}{\lambda}\right)^{n-s} \left(\frac{\lambda}{n}\right)^{n}\\ &=\lim_{n\to\infty}\frac{n!}{s!(n-s)!}\; \lambda^{s} \left(\frac{n-\lambda}{n}\right)^{n-s}\left(\frac{1}{n}\right)^{s} \\ &=\lim_{n\to\infty}\frac{n!\;\lambda^{s}}{s!(n-s)!\;n^{s}}\; \left(1-\frac{\lambda}{n}\right)^{n-s}\\ &=\lim_{n\to\infty}\left(\frac{\lambda^{s}}{s!}\right) \left(\frac{n!}{(n-s)!\;n^s}\right) \left(1-\frac{\lambda}{n}\right)^{n} \left(1-\frac{\lambda}{n}\right)^{-s} \\ \end{align*}$$ O limite do produto é o produto dos limites, desde que todos existam e sejam finitos, então vamos conferir um por um. O primeiro termo é constante, easy. No segundo, o numerador é um polinômio de grau ##n##, o denominador é um polinômio de grau ##n##, então a fração é constante e vale 1, já que os coeficientes iniciais ambos valem 1. O terceiro é a definição de ##e^{-\lambda}##. No quarto a segunda fração vai para zero, então o total tende a 1. Combinando tudo: $$\lim_{n\to\infty}p_s(n)= e^{-\lambda}\frac{\lambda^s}{s!}$$ Ahá, é uma distribuição de Poisson!

Não foi à toa que aquele histograma das figurinhas repetidas parecia familiar, elas seguem uma distribuição de Poisson! Em retrospecto até faz sentido. Se a distribuição das figurinhas é uniforme, então em um longo período de tempo, cada figurinha deve ter uma taxa de aparição constante (igual a ##\lambda=n/m##). Se a taxa é constante ao longo do tempo, então a quantidade de figurinhas em um pequeno intervalo tem que seguir a distribuição de Poisson.

Agora falta pouco para conseguir usar o teste do ##\chi^2## nas repetidas. Precisamos primeiro separar as amostras em categorias. O natural seria usar uma categoria por quantidade de repetidas, mas em teoria, um cara muito azarado poderia ter 424 figurinhas repetidas, e nesse caso o produto ##np## seria muito, muito menor que 5. Vamos então aglutinar todas as categorias de 6 repetidas em diante em um grupo só. A tabela fica assim:

Quantidade de repetições0123456 ou mais
Figurinhas, medido94150482453
Figurinhas, teórico17.140.247.437.221.910.35.9

Agora podemos calcular o ##\chi^2## medido e o teórico, lembrando que temos 6 graus de liberdade e ##\alpha## valendo 5%: $$\chi^2=\sum\frac{(Y_i-E_i)^2}{E_i}\sim 11.5275$$
InverseCDF[ChiSquareDistribution[6], 1 - 0.05]
>> 12.59
Conclusão: como 11.53 é menor que 12.59, então novamente não podemos rejeitar a hipótese de que as figurinhas são uniformes e independentes, e concluímos que sim, as figurinhas provavelmente são honestas. É verdade que tem muitas repetidas, mas o número de repetidas está dentro do esperado.

Quanto gastar para completar?


No pŕoximo post, vamos finalmente calcular qual o minimo de dinheiro que você precisa gastar para completar o álbum. Stay tuned!

segunda-feira, 19 de março de 2018

A Marcha dos Lemmings

Lemmings são pequenos roedores que vivem na tundra ártica. Eles possuem pelos longos e macios, rabo curto, são herbívoros e explodem.

Sim, eles explodem! Os lemmings possuem a capacidade de ligar um contador luminoso sobre suas cabeças. A cada passo que dão, o contador decrementa, até que finalmente explodem quando o contador chega a zero.


O que você faria se soubesse que vai explodir em seis passos? Parece que não tem muitos passeios possíveis com exatamente seis passos, mas a intuição é meio falha aqui. Na verdade, o número de passeios possíveis é exponencial no número de passos! Como fazer um algoritmo que calcula o número de caminhos possíveis que o lemming pode fazer antes de explodir?

As Regras do Jogo


Para modelar o problema, vamos supor que os passos que lemming pode fazer são descritos por um grafo orientado. Por exemplo, no grafo abaixo o lemming começa no nó A. Do A ele pode ir para B ou C, de B ele pode ir para C, e de C ele pode ir para A.


Se o nosso lemming pode dar seis passos, então existem sete caminhos diferentes que ele pode fazer:

ABCABCA
ABCACAB
ABCACAC
ACABCAB
ACABCAC
ACACABC
ACACACA

Logo ##p(6)=7##, e podemos até tabelar ##p(n)## para todos os ##n## de 0 a 6:

$$p(0..6)=1,2,2,3,4,5,7$$
Dado ##n##, como fazer um algoritmo que calcule ##p(n)##?

Eu vou mostrar três algoritmos possíveis usando a linguagem do Wolfram Mathematica. Para isso, o grafo precisa ser descrito de alguma maneira que os algoritmos entendam. Eu vou usar uma matriz de adjacências: para um grafo com ##k## nós, a matriz é ##k\times k## e o elemento na linha ##i## da coluna ##j## é 1 se existe uma ligação saindo do nó ##j## e indo em direção ao nó ##i##, e 0 em caso contrário. Para o nosso exemplo, a matriz é:

$$M=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right]$$

Força Bruta


Sempre que eu tenho que fazer um algoritmo do zero, prefiro começar com a força bruta. É fácil e rápido de implementar, e embora não seja eficiente, serve para conferir o resultado dos algoritmos mais avançados.

Nesse caso, para implementar a força bruta você passa o grafo, o nó em que o lemming está, quantos passos faltam para terminar, e o resto é uma recursão simples. Como toda recursão, não esqueça de fazer o caso base! Para esse problema, o caso base é ##p(0)=1##, ou seja, tem um único caminho possível sem dar passo nenhum, que é não sair do lugar.

No fim, a solução em força bruta fica assim:
count[graph_, pos_, 0] := 1
count[graph_, pos_, size_] := Sum[
    If[graph[[i, pos]] == 1, count[graph, i, size - 1], 0],
    {i, 1, Length[graph]}]
Vamos conferir o resultado:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}
Perfeito! Mas qual é a complexidade desse algoritmo? Bem, ele precisa visitar cada um dos caminhos possíveis, e o número total de caminhos é exponencial, logo a complexidade vai ser do tipo ##O(c^n)##, onde a constante ##c## depende do grafo. Ou seja, você não quer usar esse algoritmo quando o ##n## é grande.

Multiplicação de Matrizes


Como melhorar então? Tem uma maneira calcular ##p(n)## mais rápido usando multiplicação de matrizes. É mais fácil mostrar como funciona por exemplos: no estado inicial, nós temos um lemming no nó A, então vamos montar um vetor correspondente: ##V=[1 \;0\; 0]^T##. Se multiplicarmos ##M## por ##V##, qual o resultado?

$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}1\\0\\0\end{matrix}\right]=\left[\begin{matrix}0\\1\\1\end{matrix}\right]$$
Ou seja, um lemming que está em A vai parar em B ou C. Continuando:
$$M V=\left[\begin{matrix}0 && 0 && 1 \\ 1 && 0 && 0 \\ 1 && 1 && 0\end{matrix}\right] \left[\begin{matrix}0\\1\\1\end{matrix}\right]=\left[\begin{matrix}1\\0\\1\end{matrix}\right]$$
Dado um lemming em B ou C, o resultado é um lemming em A ou C. Se o vetor de entrada mostra quantos lemmings tem em cada nó, multiplicar pela matriz M vai mostrar onde eles podem estar depois de dar um passo.

Agora é só estender o raciocínio. Se multiplicar por M correponde a um passo, então multiplicar ##n## vezes por M é o equivalente a ##n## passos. E se você começou o processo com um vetor unitário, então a soma do elementos de V é o número de caminhos possíveis partindo do nó inicial!

Uma implementação possível é a abaixo:
count[graph_, pos_, size_] :=
  Total[MatrixPower[graph, size] . UnitVector[Length[graph], pos]]
Conferindo:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}
Bateu certinho novamente! E qual é a complexidade? Se você implementou a multiplicação e a exponenciação de matrizes do jeito mais simples, vai ser da ordem de ##O(k^3 n)##, que certamente é melhor que exponencial mas ainda não é boa o suficiente para usar na prática.

No entanto, você pode melhorar isso com implementações melhores! A multiplicação de matrizes pode usar o algoritmo de Strassen, e a exponenciação pode usar o método binário, aí o total fica ##O(k^{2.8} \log_2 n)##, bem melhor; e se a matriz for esparsa dá para ficar ainda mais rápido.

Função Geradora


Mas tem um método que supera o anterior, que é usar funções geradoras! Nós começamos tudo com um grafo orientado, e todo grafo orientado é equivalente a um autômato finito. Por sua vez, todo autômato finito é equivalente a uma gramática regular, então deve ter um jeito de descrever os caminhos do lemming usando uma gramática. Para o exemplo dado, a gramática fica assim:

$$\begin{align*}\alpha &\to A \;|\; A \beta \;|\; A \gamma \\ \beta &\to B \;|\;  B \gamma \\ \gamma &\to C \;|\; C \alpha\end{align*}$$
Todos os caminhos possíveis produzem strings que são aceitas por essa gramática, e a gramática não aceita nenhuma string que não seja um caminho possível.

O truque agora é usar combinatória analítica, que para o problema de enumerar strings de uma gramática é bem simples: basta trocar cada token por um ##z##, e cada OR por uma adição. A gramática vira um sistema de equações:

$$\begin{align*}\alpha &= z + z\beta +z \gamma \\ \beta &= z+z \gamma \\ \gamma &= z+z \alpha\end{align*}$$
Resolvendo para ##\alpha##, chegamos na função geradora:

$$\alpha=\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}$$
E para que serve a função geradora? Olhe só o que acontece quando abrimos a função em série de potências:

$$\frac{z + 2 z^2 + z^3}{1 - z^2 - z^3}=z+2z^2+2z^3+3z^4+4z^5+5z^6+7z^7+\dots$$
Ahá! Os coeficientes de ##z## são exatamente o número de caminhos no grafo! Para calcular ##p(6)##, basta olhar o coeficiente de ##z^7## (você precisa somar um porque a função geradora está contando o número de tokens na strings, e nós queremos o número de passos).

Podemos ir além agora. Quando a função geradora é uma função racional, você sempre consegue inverter e conseguir uma fórmula explícita para ##p(n)##:

$$p(n)=0.957\times 1.325^n+0.426\times 0.869^n\cos(1.469+ 2.438 n)$$
Essa fórmula é exata (ou seria, se eu não tivesse truncado os números para caber na tela). Note que o resultado sempre vai ser um inteiro, e que a parte oscilante da fórmula sempre vai ser menor que 0.4, então dá para simplificar mais ainda:

$$p(n)=\left\lfloor 0.5+0.957\times 1.325^n\right\rfloor$$
Basta pegar a parte exponencial e arredondar para o inteiro mais próximo. Podemos observar duas coisas aqui. Primeiro, agora nós temos um algoritmo onde o número de operações independe do valor de ##n##, então efetivamente esse algoritmo é ##O(1)##. Na prática você vai implementar usando bignum, aí não fica ##O(1)## de verdade; mas os outros algoritmos mostrados também vão ficar mais lentos na mesma proporção.

Segundo, nós provamos que aquele primeiro algoritmo era de fato ##O(c^n)##! Conseguimos até calcular quem é o tal do ##c##: é aproximadamente 1.325; ou, se você quiser o valor exato com radicais:

$$c=\frac{\sqrt[3]{9+\sqrt{69}}+\sqrt[3]{9-\sqrt{69}}}{\sqrt[3]{18}}$$
Para implementar esse método computacionalmente, é só reescrever aquele sistema de equações em forma matricial:

$$\left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right] = z \left[\begin{matrix}0 && 1 && 1\\0 && 0 && 1\\ 1 &&0 && 0\end{matrix}\right] \left[\begin{matrix}\alpha\\ \beta \\ \gamma\end{matrix}\right]  + z\left[\begin{matrix}1\\ 1\\1\end{matrix}\right] $$
A matriz ali no meio é a transposta da matriz de adjacências! Agora já podemos implementar:
count[graph_, pos_, size_] := With[{n = Length[graph]},
  SeriesCoefficient[LinearSolve[
      IdentityMatrix[n] - z Transpose[graph],
      ConstantArray[z, n]][[pos]],
    {z, 0, size + 1}]]
Como esperado, a resposta é correta:
graph := {{0, 0, 1}, {1, 0, 0}, {1, 1, 0}}
Print[Table[count[graph, 1, n], {n, 0, 6}]]

{1,2,2,3,4,5,7}

Um Problema em Aberto


Nós chegamos em um algoritmo que é ##O(1)##, dá para otimizar mais? Não! Esse é um limite teórico.

Quer dizer, existem algoritmos que são ##O(0)##. Por exemplo, "crie um algoritmo que, dado uma circunferência e um diâmetro, calcule a razão entre os dois". Esse é algoritmo é ##O(0)## porque a resposta é sempre ##\pi##, independe do input. Mas, se a resposta não é constante, ##O(1)## é o melhor algoritmo possível. (A não ser que você crie um algoritmo que dê a resposta antes do input chegar. Nesse caso me ensine como faz, porque você inventou a viagem no tempo).

Mas existe uma pequena variação no problema dos lemmings que o torna muito, muito mais difícil. Suponha que agora o lemming não pode passar duas vezes pelo mesmo nó. A força bruta para essa variação ainda é ##O(c^n)##. Mas ninguém bolou um algoritmo que seja melhor que a força bruta nesse caso!

Quando o grafo é finito, você pode usar a força bruta para enumerar todas as soluções, já que eventualmente o grafo acaba. Mas se o grafo for infinito, aí danou-se. Para esse caso, tem gente procurando a solução desde a década de 50. Já são setenta anos procurando o algoritmo e ninguém achou nada para um grafo genérico, nem mesmo um assintótico! O melhor resultado até hoje é um de 2011 onde acharam o ##c## para um grafo específico, o lattice hexagonal. Para todos os outros grafos, o melhor que temos são heurísticas e aproximações numéricas.

Se você estiver sem nada o que fazer no próximo fim de semana, essa pode ser uma boa maneira de ocupar o tempo!

sexta-feira, 19 de janeiro de 2018

O dia que o Knuth ganhou meu cheque

No último dia 10, o Donald Knuth completou 80 anos de vida. Para comemorar, ele fez uma festinha na remota e gelada cidadezinha de Piteå, no norte da Suécia, e eu tive a sorte de poder participar. A festa foi em formato de simpósio de matemática: ao invés de presentes, cada amigo apresentou uma palestra sobre um tema relacionado ao trabalho do Knuth. (E os amigos dele são os caras mais famosos da computação, estavam lá o Sedgewick, o Karp, o Tarjan, entre outros).



Eu também fiz minha contribuição, presenteei o Knuth com vários artigos contendo idéias para futuras revisões dos livros dele. Mas em um dos artigos aconteceu uma coisa engraçada: o Knuth notou que uma das minhas contas podia ser simplificada, e como agradecimento pela sugestão eu tive a oportunidade de inverter expectativas, e dar um cheque para o Knuth :D


A Soma dos Quadrados


O artigo em questão era uma sugestão para o Concrete Mathematics, que é o meu livro de matemática predileto. Sou fã desse livro desde os 17 anos, quando o usava como referência para estudar para a Olimpíada de Matemática.

Quem já leu sabe que o livro contém mais de dez demonstrações diferentes da fórmula para a soma dos quadrados. Para ##n\ge 0##, a fórmula é:
$$\sum_{1\le x\le n} x^2=\frac{n(n+1)(2n+1)}{6}$$
As demonstrações vão desde as mais simples (usando indução finita), até as mais complexas (usando a fórmula de Euler-MacLaurin). Porém, uns anos atrás, eu bolei uma demonstração nova que é provavelmente a mais curta demonstração de todas! Ela está na caixa azul abaixo:


Testando valores pequenos, notamos que a fórmula é verdadeira para ##n=0##, ##n=1##, ##n=2##, ##n=3##. Logo, é verdadeira para qualquer ##n##.

Sim, eu sei o que você está pensando: "Ricbit, você tá doido? Só porque funciona para alguns números não quer dizer que funciona para todos os números!".

A preocupação é válida: o que mais tem na matemática são fórmulas que funcionam para números pequenos mas falham para números grandes. Um exemplo famoso é o polinômio sortudo de Euler, ##k^2-k+41##, que produz apenas números primos para ##k=1##, ##k=2##, ##k=3##, etc., mas falha lá na frente quando ##k=41##.

Porém, nesse caso, testar números pequenos é suficiente para demonstrar a fórmula! Para entender o motivo precisamos de dois fatos sobre cálculo e polinômios.


O Cálculo Discreto


Quem já estudou cálculo com certeza sabe algumas integrais de cabeça. Por exemplo, a integral do monômio:
$$\int x^n\;dx=\frac{x^{n+1}}{n+1}$$
O que você talvez não saiba é que o cálculo que a gente aprende nas faculdades de exatas é só um tipo de cálculo: o Cálculo Contínuo. Existem outros tipos, como o Cálculo Discreto, que lida com somatórias ao invés de integrais. Esse cálculo possui várias fórmulas análogas ao cálculo mais comum. Em especial, ele tem uma fórmula análoga à integral do monômio:
$$\sum x^{\underline{n}}\;\delta x=\frac{x^{\underline{n+1}}}{n+1}$$
Nesse fórmula o monômio não usa a potência simples, ele usa a potência fatorial, que é definida da seguinte maneira:
$$x^{\underline n}=x(x-1)(x-2)\dots(x-n+1)$$
É como se fosse um fatorial, mas ao invés de ir até o fim, você pega só os ##n## primeiros termos. (Para pronunciar ##x^{\underline n}##,  você fala "x elevado a n caindo").

Para converter uma potência fatorial em uma potência normal, você pode abrir a definição e multiplicar todos os termos, mas isso dá trabalho. É mais fácil ter à mão uma tabela com os números de Stirling (que vêm em dois tipos, o primeiro tipo ##{n\brack k}##, e o segundo tipo ##{n\brace k}##). Esses números são fáceis de calcular porque eles obedecem uma regra similar ao triângulo de Pascal.  Tendo a tabela, as fórmulas abaixo fazem a conversão:
$$\begin{align*}
x^{\underline{n}}&=\sum_k{n\brack k}(-1)^{n-k}x^k\\
x^{n}&=\sum_k{n\brace k}x^{\underline{k}}\end{align*}$$
Usando as fórmulas acima e alguma paciência, você consegue demonstrar a fórmula da soma dos quadrados (converta ##x^2## em potências fatoriais, use a fórmula de somatória do cálculo discreto, depois converta de volta as potências fatoriais em potências tradicionais).

Mas isso dá muito trabalho. Ao invés disso, note que as fórmulas acima tem uma consequência interessante: com elas é possível ver que a somatória de um polinômio de grau ##n## sempre vai ser um polinômio de grau ##n+1##, e em especial a somatória de ##x^2## vai ser um polinômio de grau 3, a gente só não sabe qual polinômio ainda. Mas aí temos outro truque nas mangas!

Interpolação de polinômios


Certo, eu não sei qual é o polinômio de grau 3 que resolve a somatória. Vamos então escrever esse polinômio na forma ##P(n)=an^3+bn^2+cn+d##. Apesar de não conhecermos o polinômio, é fácil descobrir os pontos por onde ele passa, basta calcular a somatória!

Para valores pequenos de ##n##, podemos tabelar ##P(0)=0##, ##P(1)=1##, ##P(2)=1+4=5##, ##P(3)=1+4+9=14##, agora podemos achar os coeficientes usando álgebra linear:
$$\left[\begin{array}{c}0\\1\\5\\14\end{array}\right]=
\left[\begin{array}{cccc}0^3&0^2&0^1&1\\1^3&1^2&1^1&1\\2^3&2^2&2^1&1\\3^3&3^2&3^1&1\end{array}\right]
\left[\begin{array}{c}a\\b\\c\\d\end{array}\right]$$
Novamente, dá muito trabalho. Ao invés disso, notamos que um polinômio de grau ##n## fica completamente determinado por ##n+1## pontos, e agora podemos finalmente entender a demonstração inicial!

Olhe novamente para o que queremos demonstrar:
$$\sum_{1\le x\le n} x^2=\frac{n(n+1)(2n+1)}{6}$$
De um lado, temos uma somatória de polinômios de grau 2, que sabemos que é um polinômio de grau 3. Do outro lado, temos um polinômio cujo grau também é 3. Para provar que esses dois polinômios são iguais, é suficiente testar os dois polinômios em quatro pontos. Escolhemos os mais fáceis que são ##n=0##, ##n=1##, ##n=2##, ##n=3##, como a fórmula funciona para esses quatro valores, então funciona para todos os valores, QED.

O cheque do Knuth


Quando eu presenteei o artigo com essa demonstração para o Knuth, em segundos após a leitura ele comentou: "Você podia ter usado -1 aqui!".

É verdade! Essa demonstração é tão curta que dá para fazer de cabeça, mas as contas para ##n=3## ficam meio chatinhas. Ao invés disso, você pode usar -1 no lugar de 3, e aí fica bem mais fácil: do lado esquerdo a soma é vazia e dá 0, do lado direito o termo ##(n+1)## zera, logo o resultado dá zero também. Você precisa expandir o domínio da fórmula, de ##n\ge 0## para ##n\ge -1##, o que significa que a fórmula provada com -1 é até melhor que a fórmula provada com 0!

Tradicionalmente, o Knuth oferece $0x1 hexadólar para quem acha um erro em um artigo dele, e $0x0.20 para sugestões. Pela simplificação que ele sugeriu, eu achei pertinente oferecer um cheque também!


Aqui o detalhe do cheque. O Knuth emite cheques pelo banco fictício de San Serriffe (eu tenho $0x5.80 lá atualmente). Os meus cheques eu resolvi emitir pelo Banco de Dados:


O Knuth adorou o cheque! Ele ficou impressionado em como eu consegui fazer tão rápido um cheque que parece mesmo um cheque, mas o segredo é simples: quem fez fez a arte foi a Ila Fox, eu só imprimi no hotel onde estava, usando papel reciclado.

Sem dúvida foi a melhor festa de aniversário que já participei! O Knuth já tinha feito uma festa dessas quando completou 64 anos, e para manter a simetria ele disse que vai fazer outra quando completar 96 anos. Eu certamente estarei lá mais uma vez :)

segunda-feira, 25 de julho de 2016

A Busca Por Bozo Binário

Se você trabalhar por tempo o suficiente com Computação, eventualmente uma das suas atribuições vai ser entrevistar candidatos para a sua equipe. Entrevistar é uma arte: existem inúmeros truques e macetes para descobrir se é esse sujeito na sua frente é mesmo o cara que você quer como seu parceiro. O Joel tem uma boa introdução sobre o tema que vale a pena ler.

Uma das habilidades que você precisa desenvolver rápido é aprender que o candidato não pensa como você. Em especial, quando ele responde uma pergunta de um jeito diferente daquele que você esperava, não quer dizer necessariamente que ele esteja errado. E esse caso eu já vi de perto, com o algoritmo que eu apelidei de Busca Por Bozo Binário.

Essa história aconteceu faz algum tempo. Eu estava na cantina conversando sobre entrevistas, quando alguém comentou que rejeitou um candidato porque ele não conseguiu implementar uma busca binária. Poxa, isso é grave mesmo! E qual o algoritmo errado que ele fez? A resposta foi essa:

"Sorteie um número de 1 até o tamanho do vetor, e verifique o elemento que tinha esse índice. Se for igual você achou, então retorna. Se for menor, repete para o lado esquerdo; se for maior, repete para o lado direito."

Isso é errado, me disseram, porque o pior caso é ##O(n)##. De fato, se o elemento que você quer achar é o primeiro, e o gerador de aleatórios retornar sempre o último, ele vai passar pelo vetor inteiro.

Mas nesse ponto eu tive que interromper: "Peraí! Você explicitou que queria ##O(\log n)## no pior caso? Se você não falou nada no enunciado, ele pode ter entendido que era ##O(\log n)## no caso médio, e aí o algoritmo dele está certo!"

Na hora ninguém botou fé que esse algoritmo de Busca por Bozo Binário era de fato ##O(\log n)##, então eu tive que voltar para casa e escrever a demonstração. Se você tem medo de Matemática Discreta, pule a caixa azul:

Para a Busca por Bozo Binário, nós queremos calcular qual é o valor médio do número de comparações que ele realiza. E como você começa uma análise de caso médio?

Primeiro você define as características da entrada: eu vou assumir que todos os casos estão distribuídos uniformemente. Depois, é só usar a fórmula do valor esperado: você soma o número de comparações em cada caso, multiplicado pela probabilidade daquele caso ocorrer. $$ F[x] = \sum_i p[i]C[i] $$ Antes de começar, vejamos os casos extremos. Quando o vetor é vazio, não precisa de nenhuma comparação, então ##F[0]=0##. Se o vetor tiver tamanho um, então só precisa de uma comparação, ##F[1]=1##. É sempre bom ter uns casos pequenos para conferir no final.

Vamos ver agora o caso de um vetor de tamanho ##n##. Eu não sei qual número vai ser escolhido para cortar o vetor em dois, é o gerador de números aleatórios que escolhe. Então eu vou tirar a média sobre todas as escolhas possíveis. Como os ##n## valores são equiprováveis, então a probabilidade individual é ##1/n##: $$F[n] = \sum_{0\le k\lt n}\frac{1}{n}F[k,n] = \frac{1}{n}\sum_{0\le k\lt n}F[k,n] $$ Na fórmula acima, ##F[n]## é o número médio de comparações para um vetor de tamanho ##n##, e ##F[k,n]## é o número médio de comparações para o vetor de tamanho ##n##, no caso em que o corte foi feito no elemento ##k##.

Vamos calcular ##F[k,n]## agora. Esse valor depende de onde está o número que estamos procurando. Eu não sei onde está o número; como estamos calculando o caso médio, precisamos tirar a média sobre todos as posições possíveis. $$F[k,n]=\sum_{0\le i\lt n}\frac{1}{n}F[i,k,n]=\frac{1}{n}\sum_{0\le i\lt n}F[i,k,n]$$ Agora ##F[i,k,n]## significa o valor médio para um vetor de tamanho ##n##, cortado na posição ##k##, e com número final na posição ##i##.

Nós temos três casos agora. Se ##i=k##, então você só precisa comparar uma vez, a comparação vai ser verdadeira e o algoritmo termina. Se ##i\lt k##, então você compara uma vez e faz recursão para o lado esquerdo. Se ##i\gt k##, então você compara uma vez e faz recursão para o lado direito. Resumindo: $$F[i,k,n]=\begin{cases} 1+F[k] &(i \lt k) \\ 1 &(i = k) \\ 1+F[n-k-1] &(i \gt k) \\ \end{cases}$$ Podemos sintetizar tudo que vimos em uma única recorrência. $$\begin{align*} F[0] &= 0 \\ F[n] &= \frac{1}{n} \sum_{0\le k \lt n}\left( \frac{1}{n}\sum_{0\le i\lt k}\left(1+ F[k]\right) + \frac{1}{n} + \frac{1}{n}\sum_{k\lt i\lt n}\left(1+ F[n-k-1]\right) \right) \end{align*}$$ Agora é só resolver!

Antes de mais nada, vamos coletar quem é constante em relação ao índice da somatória. Todos os ##1/n## coletam, e as somatórias mais internas viram constantes. $$\begin{align*} F[n] &= \frac{1}{n^2} \sum_{0\le k \lt n}\left( \sum_{0\le i\lt k}\left(1+ F[k]\right) + 1 + \sum_{k\lt i\lt n}\left(1+ F[n-k-1]\right) \right) \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} k \left(1+ F[k]\right) + 1 + \left(n-k-1\right)\left(1+ F[n-k-1]\right) \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} k + k F[k] + 1 + n-k-1 +\left(n-k-1\right)F[n-k-1] \\ &= \frac{1}{n^2} \sum_{0\le k \lt n} n + k F[k] +\left(n-k-1\right)F[n-k-1] \\ &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] \right) \\ \end{align*}$$ Vamos focar naquela segunda somatória. Primeiro fazemos ##q=n-k-1##: $$ \sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] =\sum_{0\le n-q-1 \lt n} qF[q] $$ Agora é só manipular os índices: $$\begin{align*} 0\le n-&q-1 \lt n \\ 1-n\le -&q \lt 1 \\ -1\lt &q \le n-1 \\ 0\le &q \lt n \\ \end{align*}$$ Olha só, a segunda somatória é igual à primeira! $$ \sum_{0\le n-q-1 \lt n} qF[q] =\sum_{0\le q \lt n} qF[q] =\sum_{0\le k \lt n} kF[k] $$ Substituindo, chegamos em uma fórmula curta para ##F[n]##: $$\begin{align*} F[n] &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} \left(n-k-1\right)F[n-k-1] \right) \\ &= 1+\frac{1}{n^2} \left(\sum_{0\le k \lt n} k F[k] +\sum_{0\le k \lt n} kF[k] \right) \\ &= 1+\frac{2}{n^2} \sum_{0\le k \lt n} k F[k] \\ \end{align*}$$ A somatória ainda está atrapalhando, o ideal seria sumir com ela. Uma maneira de fazer isso é isolando: $$\begin{align*} F[n] &= 1+\frac{2}{n^2} \sum_{0\le k \lt n} k F[k] \\ \sum_{0\le k \lt n} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) \end{align*}$$ Essa última fórmula vale para todo ##n##, em especial vale também para ##n-1##: $$\begin{align*} \sum_{0\le k \lt n} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) \\ \sum_{0\le k \lt n-1} k F[k] &= \frac{(n-1)^2}{2}\left(F[n-1]-1\right) \\ \end{align*}$$ Ahá! Agora é só subtrair uma da outra que a somatória desaparece! $$\begin{align*} \sum_{0\le k \lt n} k F[k] - \sum_{0\le k \lt n-1} k F[k] &= \frac{n^2}{2}\left(F[n]-1\right) - \frac{(n-1)^2}{2}\left(F[n-1]-1\right)\\ (n-1) F[n-1] &= \frac{n^2}{2}\left(F[n]-1\right) - \frac{(n-1)^2}{2}\left(F[n-1]-1\right)\\ (n-1) F[n-1] &= \frac{n^2}{2}F[n]-\frac{n^2}{2} - \frac{(n-1)^2}{2}F[n-1]+\frac{(n-1)^2}{2}\\ \frac{n^2}{2}F[n] &= \left(\frac{(n-1)^2}{2}+(n-1)\right)F[n-1]+\left(\frac{n^2}{2} -\frac{(n-1)^2}{2}\right)\\ n^2F[n] &= \left(n^2-1\right)F[n-1]+(2n-1)\\ \end{align*}$$ Chegamos finalmente em uma recorrência sem somatórias. Melhor ainda, essa recorrência está no ponto certo para usar a técnica do summation factor! Como esse é um truque muito útil de se conhecer, eu vou fazer em câmera lenta. Você pode usar um summation factor sempre que sua recorrência for da forma abaixo: $$a_n F[n] = b_n F[n-1]+c_n$$ Esse é o nosso caso, se usarmos a correspondência abaixo: $$\begin{align*} a_n &= n^2 \\ b_n &= n^2-1 \\ c_n &= 2n-1 \end{align*}$$ O segredo do summation faction é tentar achar uma sequência mágica ##s_n## que tenha a seguinte propriedade: $$s_n b_n=s_{n-1}a_{n-1}$$ E por que isso é bom? Olha só o que acontece quando você multiplica a recorrência original por ##s_n##: $$\begin{align*} a_n F[n] &= b_n F[n-1]+c_n \\ s_n a_n F[n] &= s_n b_n F[n-1]+s_n c_n \\ s_n a_n F[n] &= s_{n-1}a_{n-1}F[n-1]+s_n c_n \\ \end{align*}$$ Agora, se você substituir ##T[n]=s_n a_n F[n]##, temos: $$\begin{align*} s_n a_n F[n] &= s_{n-1}a_{n-1}F[n-1]+s_n c_n \\ T[n] &= T[n-1]+s_n c_n \end{align*}$$ Opa, essa é fácil! Dá pra resolver de cabeça, é praticamente a definição de somatória! $$\begin{align*} T[n] &= T[n-1]+s_n c_n\\ T[n] &= T[0] + \sum_{1\le k \le n} s_k c_k \\ s_n a_n F[n] &= s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \\ F[n] &= \frac{1}{s_n a_n} \left(s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \right) \\ \end{align*}$$ Pronto, agora não é mais uma recorrência, é uma fórmula simples. A dificuldade do método é encontrar o tal do ##s_n##. Para achá-lo você pode usar a intuição, chutar, fazer um sacrifício aos deuses; ou então você pode usar o método do porradão. Nesse método você começa com ##s_n=s_{n-1}\left(a_{n-1}/b_n\right)## e abre recursivamente, chegando no monstrinho abaixo: $$s_n = \frac{a_{n-1}a_{n-2}\ldots a_2 a_1}{b_n b_{n-1}\ldots b_3 b_2}$$ Essa fórmula parece grande, e é grande mesmo. Mas na maioria das vezes tem um monte de cancelamentos, o que deixa o resultado final pequeno. Vamos ver no nosso caso como fica: $$\begin{align*} s_n &= \frac{a_{n-1}a_{n-2}\ldots a_2 a_1}{b_n b_{n-1}\ldots b_3 b_2} \\ &= \frac{(n-1)^2 (n-2)^2 \ldots 2^2 1^2} {\left(n^2-1\right) \left(\left(n-1\right)^2-1\right) \ldots (3^2-1) (2^2-1)} \\ &= \frac{(n-1)^2 (n-2)^2 \ldots 2^2 1^2} {\left((n+1)(n-1)\right) \left((n-1+1)(n-1-1)\right) \ldots (3+1)(3-1) (2+1)(2-1)} \\ &= \frac{(n-1) (n-2)(n-3) \ldots3\cdot 2\cdot 1} {(n+1) (n) (n-1)\ldots (4+1) (3+1) (2+1)} \\ &= \frac{2} {(n+1) (n) } \\ \end{align*} $$ Cancelou mesmo! Bastou lembrar que ##(a^2-b^2)=(a+b)(a-b)## que o resto saiu. Agora é só lembrar que ##F[0]=0## e substituir na fórmula: $$\begin{align*} F[n] &= \frac{1}{s_n a_n} \left(s_0 a_0 F[0] + \sum_{1\le k \le n} s_k c_k \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{2(2k-1)}{k(k+1)}\\ \end{align*}$$ Essa é a solução da recorrência. Dá para achar uma forma fechada? Nesse caso dá, desde que você tope uma forma fechada em função dos números harmônicos ##H[n]##, que são definidos como: $$H[n]=\sum_{1\le k \le n}\frac{1}{k}$$ Você começa abrindo a fração no somatório: $$\begin{align*} \frac{2(2k-1)}{k(k+1)} &= \frac{A}{k} + \frac{B}{k+1}\\ \frac{4k-2}{k(k+1)} &= \frac{A(k+1)+B k}{k(k+1)}\\ 4k-2 &= Ak+A +B k \\ 4k-2 &= (A+B)k+A \\ A &= -2 \\ B &= 6 \end{align*}$$ Agora você divide o somatório em dois: $$\begin{align*} F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{2(2k-1)}{k(k+1)}\\ F[n] &= \frac{n(n+1)}{2n^2} \sum_{1\le k \le n} \frac{6}{k+1} - \frac{2}{k}\\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - \sum_{1\le k \le n} \frac{2}{k} \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2\sum_{1\le k \le n} \frac{1}{k} \right) \\ F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2H[n] \right)\\ \end{align*}$$ O primeiro somatório precisa de um pouco de massagem: $$\begin{align*} \sum_{1\le k \le n} \frac{6}{k+1} &= \sum_{1\le q-1 \le n} \frac{6}{q} \\ &= \sum_{2\le q \le n+1} \frac{6}{q} \\ &= -6+\frac{6}{n+1}+6\sum_{1\le q \le n} \frac{1}{q} \\ &= -6+\frac{6}{n+1}+6H[n] \\ \end{align*}$$ Por fim: $$\begin{align*} F[n] &= \frac{n(n+1)}{2n^2} \left( \sum_{1\le k \le n} \frac{6}{k+1} - 2H[n] \right)\\ F[n] &= \frac{n(n+1)}{2n^2} \left( -6+\frac{6}{n+1}+6H[n] - 2H[n] \right)\\ F[n] &= -3 +\frac{2(n+1)}{n}H[n]\\ \end{align*}$$ Ufa, deu trabalho! E isso é ##O(\log n)## mesmo? Bem, podemos verificar usando um assintótico para o harmônico, por exemplo: $$ H[n] = \ln n + \gamma + O(1/n) $$ Agora você substitui: $$\begin{align*} F[n] &= -3 +\frac{2(n+1)}{n}H[n]\\ &= -3 +2H[n] +\frac{H[n]}{n} \\ &= -3 +2\ln n +2\gamma +O(1/n)+ 2\frac{\ln n}{n}+2\frac{\gamma}{n}+O(1/n^2) \\ &= 2\ln n + O\left(\frac{\ln n}{n}\right) \\ &= O(\log n) \end{align*}$$ Como queríamos demonstrar!

Ou seja, realmente a Busca por Bozo Binário é ##O(\log n)##. Em demonstrações grandes assim eu sempre gosto de fazer uma simulação por Monte Carlo para conferir as contas. A minha simulação está no github:

Busca por Bozo Binário, simulação por Monte Carlo

Simulando 50000 vezes em um vetor de tamanho 3000, o número médio de comparações foi de ##14.1769##. O valor teórico previsto pela fórmula é de ##14.1732##, nada mau!

Eu também fiz uma simulação da busca binária tradicional como comparação. Com os mesmos parâmetros, a busca tradicional usou ##10.6356## comparações no caso médio, ou seja, foi mais eficiente que o Bozo Binário.

Por que isso acontece? Analisando o assintótico fica claro. Para ##n## grande, o valor médio da Busca por Bozo é de ##2 \ln n## (logaritmo natural), enquanto que na busca binária é ##\log_2 n## (logaritmo de base 2). Então, em média, esperamos que Bozo Binário seja mais lento por um fator de:

$$\begin{align*} \frac{2\ln n}{\log_2 n} &= {2\frac{\log n}{\log e}} \div {\frac{\log n}{\log 2}} \\ &= 2\frac{\log n}{\log e} \times \frac{\log 2}{\log n} \\ &= 2\frac{\log 2}{\log e} \\ &= 2\ln 2 \\ &\simeq 1.386 \end{align*}$$ Ou seja, mais ou menos uns 38% mais lento, o que bate com a simulação.

Depois da análise do algoritmo, a dúvida que resta é: e o candidato? Bem, se fosse eu a entrevistar, precisaria coletar mais dados para saber se ele realmente sabia do tempo médio, ou se estava perdido e implementou a primeira coisa que veio na cabeça.

Um possível follow-up é "Certo, esse método é mais complicado que a busca binária tradicional, mas funciona. Você pode me falar alguma situação onde ele é preferível sobre o tradicional?". Curiosamente, existe sim um caso onde a Busca por Bozo Binário é melhor que a busca binária tradicional! Mas esse caso fica para a próxima :)

sexta-feira, 24 de julho de 2015

Retrocomputação de Verdade #1

Quando eu falo que gosto de retrocomputação, a maioria das pessoas assume que eu estou falando de TK90X e MSX. Mas para mim esses aí não são retrocomputação, são só computação. Eu usei essas máquinas quando elas eram o que havia de melhor (na minha faixa de preço). Antes delas, porém, teve uma época mágica onde computadores eram grandes, barulhentos e ocupavam salas inteiras!

Sempre que eu estou viajando por aí, aproveito a oportunidade para visitar museus de ciência e conhecer esses grandes retrocomputadores do passado. Mas só ver de perto eu não acho tão legal, o que completa a experiência é entender como eles funcionam! Minha idéia original era falar sobre todos os retrocomputadores que já conheci, mas, como ficaria muito grande, vou começar só com os retrocomputadores mecânicos:

O Mecanismo de Antikhytera, 200-100 AC.

Museu Arqueológico Nacional, Atenas, Grécia


A maior prova da engenhosidade grega, o Mecanismo de Antikhytera é muitas vezes chamado de o primeiro computador. Quer dizer, é um computador analógico não-reprogramável, hoje em dia estaria mais próximo do que a gente chama de calculadora.

E para que ele serve? O Mecanismo tem a função de calcular intervalos de tempo: você gira uma manivela que aciona o ponteiro dos anos, e ele ajusta sozinho os outros ponteiros, que são muitos. Ele calcula o mês, a fase da Lua, a posição dos planetas, os eclipses mais próximos, a data da Olimpíada, a constelação do Zodíaco e assim por diante.

Todos esses cálculos são baseados em engrenagens dentadas. Se você tem uma engrenagem com x dentes ligadas a outra com y dentes, a segunda vai rodar com velocidade angular igual a x/y da original. Então você só precisa saber a razão entre a duração de um mês e um ano para fazer um ponteiro que indica o mês.

Mas os gregos eram realmente espertos, e sabiam que algumas relações você não consegue fazer só com razões. A posição da Lua, por exemplo. Eles sabiam que a Lua anda mais devagar em alguns trechos que em outros (por causa do que hoje conhecemos como segunda Lei de Kepler), e, para simular esse efeito, eles fazem uma das engrenagens não ter o centro fixo, assim ela tem velocidade angular variável com a posição.

E como sabemos tanto sobre ela, dado que só temos pedaços de um exemplar que foi achado em um naufrágio? É que um dos pedaços era o manual!

Olhe as inscrições em grego

Essa é a parte legal de visitar museus e ver essas coisas de perto, dá para prestar atenção nos detalhes!

O Motor de Diferenças, 1822

Museu da História da Computação, Mountain View, USA


O século XIX foi a época do Maxwell, Faraday, Tesla, Edison. A engenharia elétrica estava bombando! Mas, para fazer os circuitos funcionarem, você precisa saber resolver equações diferenciais bem complicadas. Essas contas eram resolvidas com tábuas de logaritmos que simplificavam bastante o processo. Mas quem criava as tábuas de logaritmos?

Em geral elas eram feitas por um batalhão de pessoas que calculavam logaritmos na mão, com seis ou sete dígitos significativos. É um trabalho repetitivo e propenso a erros. Hoje em dia nós sabemos que a melhor maneira de realizar cálculos repetitivos e propenso a erros é um computador, mas o primeiro a ter essa idéia foi o Babbage, em 1822. Nessa época ele conseguiu o financiamento para construir o Motor de Diferenças, que tem uma única finalidade: calcular polinômios. Assim ele teria como aproximar os logaritmos por séries de Taylor equivalentes.

Depois de sugar dinheiro do governo por vinte anos sem entregar nada, o financiamento foi cortado. (O governo percebeu que com aquela dinheirama toda dava para contratar gente o suficiente para calcular os polinômios e conferir as contas). Mas foi uma pena, porque o projeto do Babbage era sólido. Se tivesse sido terminado, nunca mais alguém precisaria fazer esse tipo de conta.

Mas será que o Babbage teria conseguido terminar, se não tivesse perdido o contrato? Em 1985, cientistas resolveram seguir o projeto original e montar um Motor de Diferenças de verdade. Eles descobriram que o projeto tinha pequenos erros, mas os erros eram tão bobos que certamente foram introduzidos de propósito, como forma de proteção contra cópia (a técnica de introduzir erros como proteção contra pirataria foi usada desde a Renascença, quem costumava fazer isso era o Leonardo da Vinci). Corrigindo os erros propositais, o Motor funciona! Atualmente tem duas réplicas em funcionamento, uma em Londres e outra na Califórnia.

Por dentro, o Motor de Diferenças também funciona com engrenagens como a Mecanismo de Antikhytera, mas com uma diferença fundamental. O Mecanismo codificava a informação como um ângulo da engrenagem, se você quisesse mais precisão, precisava aumentar fisicamente o tamanho da engrenagem. O Babbage queria precisões de sete dígitos significativos, aí esse método é inviável. A solução foi abandonar a codificação analógica e usar engrenagens digitais. Cada engrenagem marca um dígito de 0 a 9, e com sete engrenagens em sequência você codifica os sete dígitos desejados.

Duas engrenagens encaixadas implementam adição ou subtração (já que virando uma delas, a outra vira junto). Mas, se são engrenagens digitais, então falta um detalhe importante: o vai-um. Eu gravei o vídeo abaixo mostrando o carry propagando, garanto que é a coisa mais bonita que você vai ver hoje!

Propagação mecânica de carry

A Bomba de Turing, 1940

Bletchley Park, UK




Sobre a história desse computador eu não preciso falar muito, já que ela foi transformada em um filme: O Jogo da Imitação (tem no Netflix). Eles precisaram transformar a história do Turing num romance heterossexual para ficar palatável para as massas, mas tá valendo, pelo menos agora as pessoas sabem quem foi o Turing.

A máquina que ele criou era conhecida entre eles como a Bomba. Ela tinha uma única função: quebrar o código de criptografia da máquina alemã Enigma, então para entender a Bomba você precisa entender o Enigma antes. As primeiras versões do Enigma era assim:

O Enigma

O coração do Enigma são os rotores (as engrenagens à esquerda, na foto). Ao longo do rotor está sequência alfabética misturada. Quando inserida no slot apropriado, o rotor faz uma cifra de permutação, obedecendo uma tabela como essa:


Ou seja, a letra A no original vira Z, a letra B vira P, e assim por diante. Quebrar esse código é fácil, qualquer criança esperta consegue. Mas o Enigma é mais seguro que isso. Ao invés de um único rotor, ele tem três rotores em sequência, e os três são diferentes: a letra A vira Z, depois no rotor seguinte você pega esse Z e vira F, depois no último o F pode virar Q. E, para complicar mais, cada vez que você aperta uma tecla o rotor gira, ou seja, dá um offset de uma letra no primeiro rotor. Quando o primeiro rotor gira por completo, induz o segundo a girar uma posição, e assim por diante.

Bastante complicado né? Mas isso não é seguro o suficiente em um contexto de guerra mundial. A senha desse sistema é a configuração inicial dos rotores. Você pode encaixar os rotores no slot em qualquer ordem (6 combinações), e iniciando em qualquer ângulo inicial (26^3 combinações). Isso dá um total de 105456 senhas possíveis. Ninguém consegue testar todas essas senhas em um dia. Mas, também, ninguém precisa testar sozinho! É só pegar 105456 soldados e mandar cada um testar uma senha diferente. Um deles vai acertar de primeira! (E 105 mil soldados não é muito, é só lembrar que no Maracanã cheio cabem 80 mil pessoas. 105 mil soldados é um Maracanã e meio só).

Por isso os alemães colocaram na frente da máquina uns plugues adicionais, que servem para permutar letras na entrada. Se você liga um plugue da posição A para a R, então a letra A vira R e a letra R vira A. Eles usavam em média uns dez plugues, então o número de senhas sobe de 105 mil para 100 bilhões. Parece uma solução esperta, mas isso, na verdade, foi uma burrice dos alemães! Desse jeito, é garantido que a letra A nunca vai virar A de novo. Então eles aumentaram o número de chaves, mas diminuíram o espaço de busca.

Agora já dá para entender a Bomba de Turing. A Bomba, na verdade, é um emulador de Enigma (o Turing sempre foi fã de emuladores, seu conceito de máquina universal é basicamente um emulador genérico). Ele tem todos os rotores iguais aos do Enigma, e um motor que gira os rotores em todas as posições possíveis, bem rapidamente. Atrás de cada rotor tem um contato elétrico: se você colocar uma palavra codificada de um lado, e uma palavra traduzida do outro, quando os rotores giraram na posição correta, dá o contato e ele avisa que achou a senha. E na prática ele não precisa girar todas as combinações possíveis. Se em alguma rotação tiver uma letra que foi mapeada nela mesma, você pode cortar a árvore de busca, porque essa senha com certeza está errada.

Então tudo que resta para quebrar o Enigma é achar um par de palavras, uma codificada, e sua correspondente não-codificada. Isso é o que hoje em dia chamamos de known-plaintext attack. Parece difícil, mas existiam vários truques para conseguir pares assim.

Um deles era baseado no tamanho da mensagem máxima dos alemães. Quando estourava o limite de tamanho, eles quebravam a mensagem em duas, e a segunda parte começava com "continuando, ...". Então, quando encontravam pacotes de várias mensagens longas, uma seguida da outra, certeza que uma delas começaria com "Fortsetzung".

Outra maneira era colocando minas explosivas no meio do mar. Os aliados colocavam as minas de um jeito fácil de encontrar, e em um lugar conhecido. Aí era só esperar os alemães acharem as minas. Quando eles achavam, transmitiam as coordenadas pelo rádio, então era só colocar as coordenadas da mina no computador, que eventualmente ele iria achar uma mensagem correspondente.

Quem viu o filme do Turing deve lembrar que no final eles jogam tudo numa fogueira. Isso aconteceu mesmo, quase tudo sobre o projeto foi destruído, e o que não foi era classificado como ultra-secreto. Só na década de 70 é que o mundo descobriu a existência das Bombas de Turing, e hoje em dia temos informação suficiente para recriar uma delas, como essa que está exposição no Bletchley Park.

A Bomba de Turing por dentro

Ainda tem mais computadores que eu conheci de perto, mas os outros ficam para a parte 2 :)