Mostrando postagens com marcador caixaazul. Mostrar todas as postagens
Mostrando postagens com marcador caixaazul. Mostrar todas as postagens

segunda-feira, 5 de setembro de 2016

A Matemática do Pokémon GO

Você já capturou todos os pokémons do Pokémon GO? Eu ainda não consegui, e também não conheço ninguém que tenha conseguido. O motivo é simples: para cada Pikachu que você encontrar, vai passar por um centena de Zubats antes! Fica difícil pegar todos quando alguns são tão raros!

pokemon_zubat

A pergunta natural agora é: se eu quiser insistir e continuar jogando até pegar todos, quanto tempo vai levar? Ou, de outro modo, qual é o número esperado de capturas até completar a coleção? Esse problema é difícil de resolver no caso geral, mas fica mais simples usando Combinatória Analítica!

Combinatória Analítica Numerada


No último post do blog eu mostrei como funciona a Combinatória Analítica, mostrando uma aplicação de cálculo com sequências de cara ou coroa. Mas essa foi só metade da história! Aquele tipo de cálculo só funciona quando os objetos contados são indistinguíveis (um resultado cara não é diferente de outro resultado cara).

Se você quiser lidar com objetos que são todos diferentes entre si, então você precisa de outra abordagem. Por exemplo, quantas permutações de ##n## elementos existem? Nesse caso, os elementos precisam ser todos diferentes. Quando ##n=3##, existem ##6## permutações:
$$ \enclose{circle}{1} \enclose{circle}{2} \enclose{circle}{3} \longrightarrow \begin{cases} \enclose{circle}{1} \enclose{circle}{2} \enclose{circle}{3} \\ \enclose{circle}{1} \enclose{circle}{3} \enclose{circle}{2} \\ \enclose{circle}{2} \enclose{circle}{1} \enclose{circle}{3} \\ \enclose{circle}{2} \enclose{circle}{3} \enclose{circle}{1} \\ \enclose{circle}{3} \enclose{circle}{1} \enclose{circle}{2} \\ \enclose{circle}{3} \enclose{circle}{2} \enclose{circle}{1} \\ \end{cases}$$
Para esse tipo de cálculo nós precisamos da Combinatória Analítica Numerada. Existem duas diferenças principais entre as versões numeradas e não-numeradas da teoria. A primeira que é, ao invés de usar funções geradoras normais, na versão numerada nós usamos funções geradoras exponenciais (EGF), que são definidas como abaixo: 
$$EGF(F[n]) = \sum_{n\ge 0}{F[n]\frac{z^n}{n!}}$$
Quando a ordem dos objetos é importante, os resultados tendem a crescer muito rapidamente, então dividir por ##n!## termo a termo ajuda as funções a ficarem bem comportadas.

A outra diferença é em um dos operadores básicos. A adição continua representando a união de conjuntos, mas a multiplicação, que antes era concatenação, fica um pouco diferente. No caso numerado, você não pode simplesmente concatenar, porque os números precisam ser sempre de ##1## a ##n##, e a concatenação simples te deixaria com duplicatas:
$$ \enclose{circle}{1} \enclose{circle}{2} \star \enclose{circle}{2} \enclose{circle}{1} =\; ?$$
A solução é renomear as bolinhas após concatenar, mas você precisa renomear de todas as maneiras possíveis que deixem a ordem dentro dos subconjuntos originais intacta:
$$ \enclose{circle}{1} \enclose{circle}{2} \star \enclose{circle}{2} \enclose{circle}{1} = \begin{cases} \enclose{circle}{1} \enclose{circle}{2} \; \enclose{circle}{4} \enclose{circle}{3} \\ \enclose{circle}{1} \enclose{circle}{3} \; \enclose{circle}{4} \enclose{circle}{2} \\ \enclose{circle}{1} \enclose{circle}{4} \; \enclose{circle}{3} \enclose{circle}{2} \\ \enclose{circle}{2} \enclose{circle}{3} \; \enclose{circle}{4} \enclose{circle}{1} \\ \enclose{circle}{2} \enclose{circle}{4} \; \enclose{circle}{3} \enclose{circle}{1} \\ \enclose{circle}{3} \enclose{circle}{4} \; \enclose{circle}{2} \enclose{circle}{1} \\ \end{cases}$$
Com isso já conseguimos definir as permutações. Elas são formadas do conjunto vazio, ou de um elemento concatenado e renomeado com uma permutação:
$$P=\varepsilon + z\star P$$
Traduzimos direto para a EGF:
$$P=1 +z P\implies P=\frac{1}{1-z}$$
E extraímos os coeficientes por comparação com a definição de EGF:
$$\begin{align*} \sum_{n\ge 0}P[n]\frac{z^n}{n!} &= \frac{1}{1-z} \\ \sum_{n\ge 0}\frac{P[n]}{n!}z^n &= \sum_{n\ge 0}1\times z^n \\ \frac{P[n]}{n!} &= 1 \\ P[n] &= n! \end{align*}$$
A quantidade de permutações de tamanho ##n## é ##n!##, então o método funciona!

O operador SEQ


Também podemos definir o operador ##SEQ##, que enumera todas as sequências possíveis de um conjunto:
$$SEQ(B[z]) = \frac{1}{1-B[z]}$$
Note que as permutações também podem ser vistas como as sequências renomeadas de zero ou mais átomos, então a definição acima dá o mesmo resultado.

Outra definição útil é ##SEQ_k##, que enumera apenas as sequências com ##k## elementos:
$$SEQ_k(B[z]) = (B[z])^k$$

O operador SET


Todas as funções acima são para construções onde a ordem é importante. Quando a ordem não importa, você sempre pode voltar para a Combinatória Analítica Não-Numerada. Mas e se o seu problema envolver os dois tipos ao mesmo tempo?

Nesse caso, a solução é o operador ##SET##. Ele é um operador numerado que define um conjunto onde a ordem não importa. A definição é a abaixo:
$$SET(B[z]) = e^{B[z]}$$
Por exemplo, quantas sequências possíveis de números de ##1## a ##n## existem, se a ordem não for importante? Isso é dado por:
$$SET(z) = e^z$$
Agora basta abrir a exponencial em série de Taylor e comparar com a definição de EGF:
$$\begin{align*} \sum_{n\ge 0}F[n]\frac{z^n}{n!} &= e^z \\ \sum_{n\ge 0}\frac{F[n]}{n!}z^n &= \sum_{n\ge 0}\frac{1}{n!}z^n \\ \frac{F[n]}{n!} &= \frac{1}{n!} \\ F[n] &= 1 \end{align*}$$
Ou seja, quando a ordem não é importante, só existe uma única sequência de ##1## a ##n## que usa todos os números de ##1## a ##n##. Faz sentido né.

(É por causa desse caso que a EGF chama função geradora exponencial. A EGF de uma sequência de ##1##s é a função exponencial.)

Pokémons uniformes


Com as ferramentas em mãos, já conseguimos modelar o problema do Pokémon GO! Mas, antes de começar o problema real, vale a pena estudar um caso mais simples. Vamos supor que os pokémons são uniformes e resolver esse caso mais fácil primeiro, depois a gente generaliza para o caso onde alguns pokémons são mais raros que outros.

O truque para usar Combinatória Analítica Numerada é achar uma maneira de enxergar seu problema como uma sequência de ##1## de ##n##. No nosso caso podemos fazer da seguinte maneira: para cada pokémon nós temos uma caixinha, e cada vez que achamos um pokémon na rua, nós colocamos dentro da caixinha uma bola com um número indicando a ordem em que ele foi encontrado.

Por exemplo, suponha que andávamos pela rua e fizemos 6 capturas, nessa ordem: Bulbassauro, Squirtle, Bulbassauro de novo, Charmander, Squirtle de novo, Bulbassauro de novo. As caixinhas ficariam assim:

bulbassauro##\enclose{box}{\enclose{circle}{1} \,\enclose{circle}{3}\, \enclose{circle}{6}}##
squirtle##\enclose{box}{\enclose{circle}{2}\, \enclose{circle}{5} }##
charmander##\enclose{box}{\enclose{circle}{4}} ##

Se nenhuma caixinha está vazia, então em algum ponto nós pegamos todos os pokémons (nesse caso, foi na quarta captura). Então, para ##n## capturas, se calcularmos a quantidade de configurações em que nenhuma caixinha está vazia, dividida pela quantidade de configurações totais, nós vamos ter a probabilidade de que a coleção foi completada com ##n## capturas ou menos.

Vamos lá. Se nós temos ##k## pokémons e ##n## capturas, então a quantidade de configurações totais é ##k^n## (é só atribuir uma pokémon a cada captura).

Para a quantidade de configurações em que nenhuma caixinha está vazia, precisamos modelar o problema com os operadores que descrevemos antes. Os pokémons são diferentes entre si, e a ordem importa, então vamos usar o operador ##SEQ##. A ordem das bolinhas dentro da caixinha não importa (afinal, ##\enclose{circle}{1} \enclose{circle}{3} \enclose{circle}{6}## é a mesma coisa que ##\enclose{circle}{3} \enclose{circle}{6} \enclose{circle}{1}## nesse contexto: a ordem das capturas é determinada pelos números nas bolinhas, não pela posição delas). Então vamos usar o operador ##SET##. E não podemos esquecer de que nenhuma caixinha pode ser vazia. Então a construção combinatória que descreve o problema é:
$$C=SEQ_k(SET(z)-\varepsilon)$$
..ou seja, uma sequência de ##k## conjuntos não-vazios. A função geradora é, portanto: 
$$C[z]=(e^z-1)^k$$
Note que, se nós já temos a função geradora, então a parte combinatória do problema já acabou. Daqui em diante é só abrir a caixa azul e fazer as contas!

Nós queremos calcular o valor esperado do número de capturas até completar a coleção. Pela definição, isso é: $$E[X]=\sum_{n\ge 0}n P(X=n)$$ Onde ##P(X=n)## é a probabilidade de que nós completamos a coleção na captura ##n##. Mas nós não temos isso! O que temos é ##P(x\le n)##, a probabilidade de completamos a coleção com pelo menos ##n## capturas. Felizmente não tem problema, porque dá para converter um no outro: $$\begin{align*} E[X] &= \sum_{n\ge 0} n P(X=n) \\ &= \sum_{n\ge 0} \sum_{i\lt n} P(X=n) \\ &= \sum_{n\ge 0} \sum_{i\ge 0} P(X=n) [i\lt n] \\ &= \sum_{i\ge 0} \sum_{n\ge 0} P(X=n) [n\gt i] \\ &= \sum_{i\ge 0} \sum_{n\gt i} P(X=n) \\ &= \sum_{i\ge 0} P(X\gt i) \\ &= \sum_{i\ge 0} 1-P(X\le i) \\ \end{align*}$$ Esse truque funciona sempre que a distribuição é discreta, porque podemos abrir a probabilidade: $$P(X>n)=P(X=n+1)+P(X=n+2)+\cdots$$ Vamos substituir agora. A nossa probabilidade para ##X\le n## é a divisão da quantidade de configurações com caixinhas não-vazias (que é o coeficiente ##z^n## da função geradora, vezes ##n!## porque é uma função geradora exponencial) pela quantidade total (que é ##k^n##): $$P(X\le n)=\frac{n![z^n](e^z-1)^k}{k^n}$$ Logo a expectativa é: $$\begin{align*} E[X] &=\sum_{n\ge 0} 1-P(X\gt n) \\ &=\sum_{n\ge 0} 1-\frac{n![z^n](e^z-1)^k}{k^n} \\ &=\sum_{n\ge 0} 1-n![z^n](e^{z/k}-1)^k \end{align*}$$ Esse último passo usou uma propriedade de escala das funções geradoras. Se uma sequência ##f(n)## tem função geradora ##F(z)##, então: $$\begin{align*} \sum_{n\ge 0}f(n)z^n &= F(z) \\ \sum_{n\ge 0}k^n f(n)z^n &= \sum_{n\ge 0}f(n)(kz)^n = F(kz) \end{align*}$$ Voltando: nós vimos que a função geradora exponencial da sequência constante de ##1##s é ##e^z##, então dá pra jogar para dentro: $$\sum_{n\ge 0} 1-n![z^n](e^{z/k}-1)^k = \sum_{n\ge 0} n![z^n]\left(e^z -\left(e^{z/k}-1\right)^k\right)$$ Para prosseguir agora tem um truque muito bom, que funciona assim: se você tem uma função geradora ##F(z)## bem comportada, então vale que: $$\sum_{n\ge 0}n![z^n]F(z)=\int_0^\infty e^{-t}F(t)\,dt$$ Parece meio mágico, mas é só consequência da forma integral do fatorial: $$n!=\int_0^\infty t^n e^{-t}\,dt$$ Se nós chamarmos os coeficientes da função geradora de ##f_n##, então: $$\begin{align*} \sum_{n\ge 0}n![z^n]F(z) &=\sum_{n\ge 0}n! f_n = \sum_{n\ge 0} f_n n! \\ &=\sum_{n\ge 0}f_n\left(\int_0^\infty t^ne^{-t}\,dt\right) \\ &=\int_0^\infty \left(\sum_{n\ge 0}f_n t^n e^{-t}\right)\,dt\\ &=\int_0^\infty e^{-t}\left(\sum_{n\ge 0}f_n t^n \right)\,dt\\ &=\int_0^\infty e^{-t}F(t)\,dt\\ \end{align*}$$ Aplicando ao nosso caso: $$\begin{align*} \sum_{n\ge 0} n![z^n]\left(e^z -\left(e^{z/k}-1\right)^k\right) &= \int_0^\infty e^{-t} \left(e^t -\left(e^{t/k}-1\right)^k\right) \,dt \\ &= \int_0^\infty \left(1 -e^{-t}\left(e^{t/k}-1\right)^k\right) \,dt \\ &= \int_0^\infty \left(1 -\left(e^{-t/k}\right)^k\left(e^{t/k}-1\right)^k\right) \,dt \\ &= \int_0^\infty \left(1 -\left(1-e^{-t/k}\right)^k\right) \,dt \\ \end{align*}$$ Nesse ponto já daria para jogar num solver numérico e calcular a solução. Mas essa integral tem uma forma fechada! Para calcular, basta fazer a mudança de variável abaixo: $$\begin{align*} v &= 1-e^{-t/k} \\ 1-v &= e^{-t/k} \\ e^{t/k} &= \frac{1}{1-v} \end{align*}$$ Agora deriva para achar quanto vale ##dt##: $$\begin{align*} \frac{dv}{dt} &= \frac{e^{-t/k}}{k} \\ dt &= k \, e^{t/k}\, dv \\ dt &= \frac{k}{1-v} dv \end{align*}$$ Arruma os intervalos: $$\begin{align*} t=0 &\implies v=1-e^0=1-1=0 \\ t=\infty &\implies v=1-e^{-\infty}=1-0=1 \end{align*}$$ E por fim soca tudo na integral: $$\begin{align*} \int_0^\infty \left(1 -\left(1-e^{-t/k}\right)^k\right) \,dt &= \int_0^1 \left(1 -v^k\right) \frac{k}{1-v} \,dv \\ &= k\int_0^1 \frac{1-v^k}{1-v} \,dv \\ &= k\int_0^1 (1+v+v^2+\cdots+v^{k-1}) \,dv \\ &= k \left[v+\frac{v^2}{2}+\frac{v^3}{3}+\cdots+\frac{v^k}{k} \right]_0^1 \\ &= k \left[1+\frac{1}{2}+\frac{1}{3}+\cdots+\frac{1}{k} \right] \\ &= k H[k] \\ \end{align*}$$ Esse é o resultado final. O número esperado de capturas até você completar sua coleção de pokémons é ##k## vezes o harmônico de ##k##!

Qual seria o número esperado de capturas então? O jogo tem 151 pokémons, mas desses tem seis que ainda não podem ser pegos: Ditto, Articuno, Zapdos, Moltres, Mewtwo e Mew. Para pegar os 145 restantes, você precisaria de ##145\times H(145)\sim 806## capturas.

Mas isso só vale se os pokémons fossem uniformes! Como alguns são mais raros que outros, precisamos melhorar nossa estimativa.

Pokémons raros


A grande vantagem da Combinatória Analítica é que ela generaliza muito fácil. O problema de completar uma coleção com elementos não-uniformes é bem difícil com técnicas tradicionais, mas aqui fica simples: basta trocar os átomos ##Z_i## por ##p_iZ_i##, onde ##p_i## é a probabilidade de ocorrência dele.

Por que isso funciona? Suponha que você tem três pokémons, de probabilidades ##p_1=1/6##, ##p_2=1/3## e ##p_3=1/2##. Ao invés de fazer ##SEQ_3(SET(z)-\varepsilon)##, como no caso uniforme, você pode fazer:
$$(SET(Z_1)-\varepsilon)\star (SET(Z_2+Z_2)-\varepsilon)\star (SET(Z_3+Z_3+Z_3)-\varepsilon)$$
Isso vai dar um resultado seis vezes maior que o esperado, mas na proporção correta. Se você jogar esse valor ##1/6## para dentro, a conta fica:
$$\begin{align*} &(SET(\frac{Z_1}{6})-\varepsilon)\star (SET(\frac{2Z_2}{6})-\varepsilon)\star (SET(\frac{3Z_3}{6})-\varepsilon)=\\ &(SET(p_1Z_1)-\varepsilon)\star (SET(p_2Z_2)-\varepsilon)\star (SET(p_3Z_3)-\varepsilon) \end{align*}$$
Exatamente como mencionamos! Isso mostra que o resultado é válido para qualquer probabilidade racional, e você pode fazer um sanduíche se quiser reais também.

Daqui em diante as contas são as mesmas do caso uniforme, mudando somente o interior da integral: 
$$\begin{align*} E[X] &=\int_0^\infty\left( 1-\left(1-e^{-p_1t}\right)\left(1-e^{-p_2t}\right)\cdots\left(1-e^{-p_kt}\right) \right)\,dt\\ &=\int_0^\infty\left( 1- \prod_{1\le i\le k}\left(1-e^{-p_i t}\right) \right)\,dt\\ \end{align*}$$
Infelizmente essa integral não tem forma fechada, precisa ser calculada numericamente mesmo. Usando uma tabela com as probabilidades estimadas de cada pokémon, eu calculei a integral, e o resultado é que você precisaria de 51568 capturas para completar sua coleção!

Como de costume, eu fiz uma simulação para conferir os valores. Os scripts estão no github, e os resultados são:

Distribuição uniforme, valor teórico 805.8, computacional 805.1
Distribuição não-uniforme, valor teórico 51568, computacional 50931

Os resultados bateram como previsto! Vale a pena também estimar o tempo: se você pegar um pokémon a cada 3 minutos, então para fazer 51568 capturas você precisaria jogar por 107 dias seguidos, sem parar para dormir!

Na prática você pode completar em menos tempo que isso usando outros recursos do jogo, como evoluções e ovos, mas não tem muito segredo, para ser um mestre pokémon você precisa se dedicar bastante :)

PS: Como quase não tem literatura em português sobre o tema, as traduções ainda não são padrão. Eu traduzi Labelled Analytic Combinatorics como Combinatória Analítica Numerada, porque não achei nenhuma outra tradução que não fosse feia demais.

sábado, 30 de julho de 2016

Totorial de Combinatória Analítica

Se você jogar Cara ou Coroa várias vezes em sequência, qual a chance de nunca aparecer duas caras seguidas?

Esse é um problema que você resolve com ferramentas da Combinatória. Se você estudou para o vestibular, certamente deve lembrar que esse tipo de problema usualmente se resolvia com permutações, arranjos ou combinações.

Infelizmente, isso é verdade no vestibular, mas não na vida real. Os problemas avançados de Combinatória resistem a essas ferramentas básicas. Em geral, você vai precisar de recorrências complexas e somatórias complicadas; e até recentemente não tinha nenhum método que fosse simples e geral, para usar em qualquer situação.

A boa notícia é que isso mudou! Nos últimos anos, surgiu um novo ramo da matemática, chamado de Combinatória Analítica, que permite calcular esses problemas com facilidade. Como ainda não existe nenhum texto introdutório em português, resolvi escrever um totorial para mostrar como as coisas ficam bem simples com esse método!
totoro

Um exemplo simples


Como funciona esse método novo? A analogia mais simples é com desenho geométrico. Com as técnicas tradicionais, você precisa de régua, compasso e uma boa dose de insight para resolver os problemas. Ao invés disso, você pode usar geometria analítica: o problema é transformado em uma série de equações, aí você não precisa pensar, só resolver (e em muitos casos nem precisa resolver, é só jogar num sistema de computação automática como o Wolfram Alpha que ele resolve para você).

Com combinatória analítica é parecido: você vai transformar a descrição do seu problema de contagem em uma equação, e aí existem técnicas padrão para resolvê-la. Para mostrar como funciona, vamos pegar um problema simples: Quantas strings binárias de tamanho ##n## existem?

Para ##n=3##, por exemplo, existem oito strings: ##000##, ##001##, ##010##, ##011##, ##100##, ##101##, ##110## e ##111##.

A primeira coisa a ser feita é caracterizar o conjunto de todas as strings. Podemos recursivamente construir as strings binárias válidas com três casos:

Caso 1: A string vazia (##\varepsilon##), ou...

Caso 2: Uma string que começa com ##0##, seguida de uma string válida, ou...

Caso 3: Uma string que começa com ##1##, seguida de uma string válida.

Pronto, agora podemos agora escrever a equação que descreve o conjunto ##B## das strings binárias. Em combinatória analítica, a regra é que a descrição "ou" vira adição, e a operação "seguida de" vira multiplicação. Vou também substituir o caractere ##0## por ##Z_0## e ##1## por ##Z_1## só para deixar claro que nesse contexto são átomos, e não números:
$$B=\varepsilon+Z_0\times B+Z_1\times B$$
Essa é a construção combinatória correspondente às strings binárias. Agora o pulo do gato: eu vou trocar ##\varepsilon## pelo número ##1##, cada átomo individual pela variável ##z##, e resolver a equação:
$$\begin{align*} B&=\varepsilon+Z_0\times B+Z_1\times B\\ B&=1+zB+zB\\ B&=\frac{1}{1-2z} \end{align*}$$
Chegamos então no elemento principal da combinatória analítica: a função geradora do conjunto das strings binárias. E por que ela é tão importante? É que essa função simples tem escondida dentro dela a solução do problema para qualquer ##n##! Basta expandir a função em série (nesse caso, basta reconhecer que essa é a fórmula da soma infinita de uma PG com razão ##2z##): 
$$B=\frac{1}{1-2z}=1+2z+4z^2+8z^3+16z^4+\cdots$$
Olha só, a solução aparece na série! O termo ##8z^3## significa que, para ##n=3##, a solução é ##8##, exatamente como tínhamos visto enumerando as possbilidades.

O operador SEQ


O método acima foi simples né? Mas dá para ficar mais fácil ainda! O pessoal que estuda combinatória analítica faz bastante tempo inventou uma série de operadores que deixam a solução mais fácil de escrever.

O primeiro deles é o operador ##SEQ##, que serve para definir sequências. Ao invés da definição recursiva com três casos, podemos simplesmente dizer que uma string binária é formada por uma sequência de zero ou mais átomos, onde os átomos possíveis são ##0## e ##1##. Daí a construção combinatória sai direto:
$$B=SEQ(Z_0+Z_1)$$
E como transformar isso em função geradora? Se você sabe que um conjunto ##A## tem função geradora ##A(z)##, então a função geradora da sequência de ##A## é: 
$$SEQ(A)=\frac{1}{1-A(z)}$$ 
No nosso caso, a função geradora das strings binárias é:
$$SEQ(Z_0+Z_1)=SEQ(z+z)=\frac{1}{1-2z}$$
Pronto, a função geradora saiu direto, sem precisar resolver nenhuma equação!

O teorema de transferência


Ainda resta o problema de abrir a função geradora em série. Nesse caso foi fácil porque conseguimos ver de cara que a função era a soma de uma PG, mas e em casos mais complicados?

Bem, resolver de maneira exata continua sendo um problema difícil até hoje. Mas na maioria das vezes você não precisa da solução exata, um assintótico para ##n## grande já é suficiente. E para isso a combinatória analítica possui uma série de teoremas de transferência, que dão a resposta assintótica da função geradora sem precisar abrir a série!

O teorema de transferência usado nesse caso é simples de usar, embora sua definição seja meio assustadora:

Se a sua função geradora for uma função racional ##f(z)/g(z)##, onde ##f(z)## e ##g(z)## são relativamente primas, ##g(0)\neq 0##, ##g(z)## tem um único pólo ##1/\beta## de módulo mínimo, e ##\nu## é a multiplicidade desse pólo, então um assintótico para a função é dado por: $$[z^n]\frac{f(z)}{g(z)}=\nu\frac{(-\beta)^\nu f(1/\beta)}{g^{(\nu)}(1/\beta)}\beta^n n^{\nu-1}$$

Eu juro que é simples, apesar de assustador. Para o nosso caso, ##f(z)=1##, ##g(z)=1-2z##, ##g'(z)=-2##, ##\beta=2##, ##f(1/2)=1##, ##g'(1/2)=-2## e ##\nu=1##. Substituindo os valores, temos:
$$[z^n]B(z)=1\times\frac{(-2)^1\times 1}{-2}\times(2)^n n^{1-1} =2^n $$
Portanto, existem ##2^n## strings binárias de tamanho ##n##, o que bate com nosso resultado de ##8## para ##n=3##.

Outro exemplo


E se eu quiser calcular quantas strings binárias existem, mas só as que não possuem dois zeros seguidos? Esse problema é complicado com combinatória tradicional, mas com combinatória analítica fica simples!

Primeiro, vamos caracterizar essas strings. Podemos descrevê-las como uma sequência de ##1## ou ##01##, seguida de ##0## ou vazio:
$$C=SEQ(Z_1+Z_0Z_1)\times(Z_0 + \varepsilon)$$
Traduzindo para função geradora:
$$C=SEQ(z+z^2)\times(z+1)=\frac{z+1}{1-z-z^2}$$
Agora aplicamos o teorema de transferência. ##f(z)=z+1## e ##g(z)=1-z-z^2##. As raízes de ##g(z)## são ##0.618## e ##-1.618##, a de menor módulo é ##0.618## com multiplicidade ##1##. Então:
$$\begin{align*} \beta &= 1/0.618 = 1.618 \\ f(1/\beta)&=f(0.618)=1.618 \\ g'(z) &=-2z-1 \\ g'(1/\beta) &=g'(0.618)=-2.236 \\ C[n]&\sim \frac{-1.618\times 1.618}{-2.236}\times 1.618^n \\ C[n] &\sim 1.171\times 1.618^n \end{align*}$$
Prontinho, transformamos raciocínio, que é díficil, em álgebra, que qualquer chimpanzé bêbado consegue fazer!

A pergunta original


Podemos agora voltar à pergunta original: se você jogar Cara ou Coroa várias vezes em sequência, qual a chance de nunca aparecer duas caras seguidas?

Ora, essa probabilidade é exatamente igual ao número de strings binárias que não possuem dois zeros seguidos, dividida pelo número total de strings binárias. Esses são os dois casos que analisamos, logo:
$$p\sim\frac{1.171\times1.618^n}{2^n} = 1.171\times 0.809^n$$
Easy! E funciona mesmo? Eu fiz uma simulação computacional para comparar com o nosso assintótico, eis o resultado:

Para ##n=5##, simulação ##40.78\%##, assintótico ##40.57\%##.

Para ##n=10##, simulação ##14.17\%##, assintótico ##14.06\%##.

Para ##n=15##, simulação ##4.84\%##, assintótico ##4.87\%##.

Script com a simulação no github

Ou seja, funcionou super bem, e a precisão vai ficando melhor conforme o ##n## aumenta.

Para saber mais


Se você gostou do que viu, a referência definitiva sobre o tema é o livro Analytic Combinatorics, do Flajolet, que foi o inventor da técnica. Não é o livro mais didático do mundo, mas está inteiramente disponível na web. No Coursera, de tempos em tempos, o Sedgewick faz um curso online, que eu super recomendo, foi onde eu aprendi inclusive. E na wikipedia... o tópico é tão novo que não tem na wikipedia em português ainda! Fica como exercício para o leitor criar uma página na wikipedia sobre o tema :)

terça-feira, 19 de maio de 2015

Livros de Colorir são Equivalentes a Sudoku

De tempos em tempos a mídia impressa inventa uma modinha para vender mais papel. A modinha da vez são os livros de colorir. Eles são uma sensação! Os livros não param nas lojas, existem grupos e grupos falando disso na internet, rola até competição para ver quem pinta melhor!

Mas, se você tem memória boa, deve lembrar que mesma coisa aconteceu dez anos atrás. Ao invés de livros de colorir, a modinha era Sudoku. Era uma sensação! Haviam livros e livros nas lojas, todo mundo falava disso na internet, rolava até competição para ver quem resolvia mais rápido!

À primeira vista, os dois passatempos são bem distintos: o livro de colorir precisa de sensibilidade artística, o sudoku de raciocínio lógico. Mas hoje eu vou provar que, na verdade, os dois são equivalentes! E tudo começa com a regra não escrita dos livros de colorir.


A Regra dos Livros de Colorir


Livros de colorir tem uma regra implícita. Todo mundo concorda que isso aqui é uma pintura válida, certo?

colorircerto (1)   

Mas isso aqui é trapacear:

colorirerrado 

A regra implícita dos livros de colorir é que pintar duas regiões adjacentes com a mesma cor não vale. Afinal, tem uma linha preta ali, e ela está ali por uma razão, que é separar regiões de cores diferentes. Se você ignorar a linha preta, não está aproveitando todos os detalhes da figura.

Se você for a uma papelaria atualmente, verá que todas as caixas de lápis de cor com 48 cores estão esgotadas. O motivo é que os livros de colorir tem muitas regiões adjacentes, então quanto mais cores, mais fácil de pintar.

E qual é o mínimo de cores que você precisa? Um estojo de 12 cores é suficiente? Incrivelmente, doze cores dá e sobra. Na verdade, com apenas 4 cores você consegue pintar qualquer desenho, sem deixar regiões adjacentes com a mesma cor! Essa proposição é o lendário Teorema das Quatro Cores, que ficou famoso por ter sido o primeiro grande teorema da matemática provado por computador.

Além disso, o problema de pintar um desenho com a menor quantidade de cores possíveis não é apenas teórico. Ele tem várias utilidades práticas! Por exemplo, compiladores usam pintura para fazer alocação de registros, e companhias aéreas usam pintura para fazer alocação de aviões em rotas aéreas.

O problema da pintura mínima também serve para resolver Sudokus. Mas antes de mostrar como isso pode ser feito, vamos resolver um problema mais simples, o Bidoku.

Bidoku


O Bidoku é a versão fácil do Sudoku. Você tem um grid 2x2, e precisa preenchê-lo com 0 e 1, de modo não tenha dois números iguais na mesma linha e na mesma coluna:

bidoku 

Esse puzzle é bem fácil, dá para resolver de cabeça né? Mas vamos ver como é o raciocínio usando pintura. Para isso, vamos dar um nome para cada uma das células do puzzle:

bidoku-letter 

Cada uma das letras pode ser 0 ou 1. Vamos enumerar as possiblidades: A pode ser 0 ou 1, B pode ser 0 ou 1, e assim por diante:

bidokugraph1 

Agora vamos ligar com uma linha todos os nós que não podem acontecer ao mesmo tempo. Começando do B, nós sabemos que ou o B é zero, ou o B é um. Então coloco uma linha entre B0 e B1. Além disso, se B for 0, então nem A e nem C podem ser 0. Da mesma maneira, se B for 1, então nem A e nem C podem ser 1. Desenhando:

bidokugraph2 

Agora é só repetir o que fizemos com o B para todas as outras letras:

bidokugraph3 

Este é o grafo equivalente do Bidoku. Com o grafo em mãos, podemos achar a figura de colorir equivalente, basta achar um desenho onde as regiões adjacentes são exatamente aquelas indicadas pelas linhas do grafo. Qualquer desenho que satisfaça essa condição serve, pode ser esse aqui, por exemplo:

bidokudrawing 

Com o desenho em mãos, podemos resolver o Bidoku através da coloração!. Começa assim: você pega todos os números que foram dados, e pinta de verde. No Bidoku exemplo acima, era dado que o A valia 1, então nós pintamos de verde o A1:

bidokugreen 

Agora é só pintar o restante com o menor número de cores possível, respeitando a regra de que regiões adjacentes precisam ter cores distintas. Se você fizer a pintura, vai notar que é possível pintar com apenas duas cores, de uma única maneira:

bidokufull 

Agora você pega todas as regiões verdes e retorna para o Bidoku original:

bidokusolved 

Funciona mesmo! O Bidoku foi resolvido pela coloração mínima!

Uma dúvida pertinente é o que acontece com a coloração se você tentar pintar um Bidoku insolúvel. Por exemplo, o Bidoku abaixo não tem solução (o B não pode ser 0 e nem 1):

bidokuunsolvable 

Se você tentar a coloração mínima, vai notar que são necessárias no mínimo três cores para respeitar as regras:

bidokuthree 

Essa regra vale sempre, se a sua coloração mínima tem mais de duas cores, então o Bidoku original era insolúvel.

Agora que já dominamos o Bidoku, podemos passar para a solução do Sudoku.

A Coloração do Sudoku


O Sudoku usa o mesmo raciocínio do Bidoku, mas o problema é muito maior. No Bidoku a célula A tinha duas possibilidades, então nós criávamos dois nós: A0 e A1. No Sudoku, ela tem nove possibilidades, então a célula A gera os nós A1, A2, A3, .., até A9. E você precisa adicionar linhas ligando os números, as linhas, as colunas, e as caixas (eu chamo de caixas aqueles blocos 3x3). O grafo final do Sudoku, levando em conta tudo isso, fica assim:

sudokugraph
Clique na imagem para ver o SVG original

São 729 nós e 11664 linhas, só um pouquinho maior que o Bidoku, que tinha 8 nós e 12 linhas. Além disso, no Bidoku nós conseguimos fazer um desenho equivalente ao grafo. No Sudoku isso não é possível, porque esse grafo não é planar (ou melhor, você pode fazer um desenho, mas vai ser um desenho multidimensional em uma geometria possivelmente não-euclideana). Isso não afeta nossa proposição, já você pode colorir diretamente os nós do grafo, ao invés de colorir as regiões de um desenho equivalente.

E como fazemos isso? Da mesma maneira que o Bidoku, você parte de um Sudoku que você quer resolver, por exemplo esse aqui, que eu peguei no The Guardian:

SUD1235E_2704 

Nós vamos pintar de verde todas as células correspondentes no grafo do Sudoku (eu não vou mostrar as linhas dessa vez, assim fica mais visível):

sudoku1
Clique na imagem para ver o SVG original

Pois bem, a minha proposição é que se você fizer a pintura mínima nesse grafo, o resultado do Sudoku vai aparecer nas células que forem pintadas de verde. Nós vimos um exemplo de que isso funciona no Bidoku, mas um exemplo não é uma demonstração, pode ter sido coincidência. Eu vou provar abaixo que isso sempre é verdade. Se você tem medo de demonstrações, pule a caixa azul para ver o Sudoku pintado.

Antes de mais nada, quantas cores precisamos para pintar esse grafo? O instinto é falar que precisamos de 4 cores, por causa do Teorema das Quatro Cores. Mas esse teorema só funciona se o grafo for planar! E o grafo do Sudoku claramente não é planar, já que ele admite o K5 como subgrafo, e, portanto, viola o teorema de Kuratowski para grafos planares.

Vamos tentar achar qual o número de cores então. Pegue uma célula qualquer, por exemplo uma célula A. Nós temos 9 nós correspondentes, de A1 até A9, e nenhum desses nós pode ter a mesma cor que o outro. Então o grafo da célula A é um clique de tamanho 9. Como o número cromático de um grafo é sempre maior ou igual ao grau do maior clique, então sabemos que o grafo tem pelo menos 9 cores.

Mas, se esse grafo foi construído a partir de um Sudoku que tem solução, então nós podemos construir uma coloração com 9 cores para ele. Você começa pintando todos os nós da solução com a cor 0. Depois, para cada tripla (i, j, digito) que foi pintada com a cor 0, você pinta a tripla (i, j, digito + 1) com a cor 1 (o incremento no dígito é mod 9). Essa pintura não viola o grafo, porque ele é simétrico nos dígitos. Repetindo o processo para as cores 2..9, você tem um grafo totalmente pintado e sem nenhuma violação.

Se o número mínimo de cores é 9, e nós sempre conseguimos construir uma pintura com 9 cores, então o número cromático é 9. Done.

Ainda falta um detalhe: ninguém garante que existe uma única pintura possível com 9 cores (e na prática vai ter mais de uma mesmo). Mas todas as pinturas possíveis resolvem o Sudoku! Como o número cromático é 9, então cada célula (i,j) vai ter um elemento pintado com cada uma das cores possíveis (já que a célula é um clique de tamanho 9). Como cada célula tem todas as cores, então cada célula vai ter um elemento com a cor 0; e se todas as células tem a cor 0, então toda pintura possível embute uma solução do Sudoku. QED.

Se você tiver paciência de fazer a pintura mínima naquele grafo gigante, o resultado será esse aqui:

sudoku3

Clique na imagem para ver o SVG original

Transcrevendo as células verdes para o diagrama original, temos a solução do Sudoku, como prometido!

solved 

Na prática, esse é provavelmente um dos piores métodos para resolver o Sudoku. Mesmo um algoritmo no estado da arte em graph coloring vai levar um bom tempo para achar uma solução através de pintura mínima. Mas pelo menos você pode dizer que Sudoku é equivalente a um livro de pintura, desde que você não ligue de fazer pinturas mínimas em livros multidimensionais de geometria possivelmente não-euclideana :)

  ricbit_ilafox_pintando2  

segunda-feira, 2 de junho de 2014

A Equação da Torre Eiffel

No mês de maio eu completei cinco anos de casado! Quando eu casei, aproveitei a lua-de-mel em Paris para falar sobre a ciência da Torre Eiffel, então nada mais apropriado que aproveitar o aniversário para falar sobre outro mistério da torre: a equação que descreve o seu formato!
  mini_ilafox_ricbit_paris  

A torre mais alta do mundo


O desafio proposto pelo Gustave Eiffel era grandioso: construir a torre mais alta do mundo. Na época, o recorde era o Monumento de Washington, com 169 metros. O Eiffel pretendia erguer sua torre com 324 metros, quase o dobro do recorde anterior.

Para entender o desafio, vale a pena comparar a Torre Eiffel com outras torres da antiguidade. O desenho abaixo está em escala, e mostra a Torre Eiffel comparada com a Grande Pirâmide de Giza no Egito, e com o Campanário de São Marcos em Veneza:

compare
Imagem: Wikipedia (CC-SA)

O primeiro desafio ao fazer uma torre é garantir que ela consiga suportar seu próprio peso. Os egípcios resolveram isso com força bruta: usaram blocos de pedra maciços e fizeram a base bem maior que o topo. O resultado final certamente funcionou, a pirâmide está aí de pé 4500 anos depois. Mas ela é um desperdício de material, com tecnologia moderna você não precisa de tanta matéria-prima para chegar nessa altura.

O Campanário de São Marcos é o extremo oposto. Ele foi construído em 1514, e nessa época a tecnologia já permitia criar uma torre com a largura da base igual à do topo. O problema é que a torre caiu! Em 1902 ela não resistiu ao próprio peso e colapsou (nenhuma vida humana foi perdida, mas o gato do zelador morreu). A torre que você visita hoje em Veneza é uma reconstrução.

Como fazer, então, uma torre alta que suporte o próprio peso? O Eiffel sabia a resposta porque tinha experiência no assunto. Antes de criar os ícones pelo qual é mais lembrado, ele construía pontes. Existem pontes do Eiffel por todo o mundo, desde a França e a Romênia, até o Vietnã e o Chile.

Para entender a tecnologia do Eiffel, vale a pena comparar com uma ponte antiga: a Pont du Gard, um aqueduto romano construído dois mil anos atrás. Ele cobre uma distância de 275 metros a uma altura de 48 metros. Para isso, ela precisa de cinco apoios no chão e três níveis de arcos para aguentar seu próprio peso:

Pont_du_gard_v1_082005
Imagem: Wikipedia (CC-SA)

Agora, compare com uma das mais belas pontes construídas pelo Eiffel: a Ponte Dona Maria, na Cidade do Porto, em Portugal. Ela cobre 353m a 60m de altura, em ambos os casos maior que a ponte romana. E faz isso usando só dois apoios no chão! O segredo do Eiffel para conseguir esse feito é a magia da treliça. Ao usar uma treliça, você consegue suportar mais peso com menos material.

Ponte_Maria_Pia_-_Porto Imagem: Wikipedia (CC-SA)

Todo mundo sabia que o Eiffel era expert no assunto, e que ele saberia construir uma torre que suportasse o próprio peso. Mas, ainda assim, na época ele ouviu muitas críticas. Uma coisa é você construir uma ponte na horizontal, outra é construir uma torre na vertical. Estruturas na vertical tem um problema adicional: será que o Eiffel sabia como fazer a torre resistir ao vento?

O método geométrico do Eiffel


Quando a torre ficou pronta, o Eiffel calou os críticos. Quem já visitou a torre sabe que no terceiro piso venta. Venta muito mesmo! Porém, a torre não se abala. Fica lá, imóvel, mesmo com os ventos mais fortes.

E qual foi o truque do Eiffel para deixar a torre tão estável com o vento? Até pouco tempo atrás, ninguém sabia. As plantas da torre Eiffel são conhecidas, você pode ver todas as especificações em um livrão publicado pela Taschen (todas mesmo, tem até o diâmetro e comprimento de cada parafuso). Mas as plantas não falam o truque da estabilidade, tudo que dizem é que o Eiffel determinou o formato da torre de forma a minimizar a ação do vento, e fez isso usando "métodos geométricos baseados em observações empíricas".

 DSC09174 

O formato da torre foi um mistério por mais de um século, até que Patrick Weidman, professor da Universidade de Colorado, matou o problema, e da maneira mais simples possível. Ele foi atrás dos diários do Eiffel e achou lá a dica que resolve o problema! O trecho do diário dizia:

"Suponha que as faces de uma treliça simples formam uma parede resistente às forças cisalhantes do vento, cujas componentes horizontais são P', P'', P''', P''''. Sabemos que, para calcular as forças atuantes no corte pelo plano MN, precisamos determinar a resultante P de todas as forças exteriores atuando acima da secção, e decompor essa resultante nas três forças passando pelas peças cortadas. Se a forma do sistema for tal que, para cada corte horizontal MN, a extensão das peças externas da treliça se intersectam na força exterior P, então as forças na barra interna serão zero e podemos excluir esse membro."

O diário tinha esse desenho acompanhando:

trelica_original

Lendo com cuidado dá para entender direitinho a motivação do formato da torre. Vamos recriar o raciocínio em câmera lenta. Nós começamos supondo que o vento é uma carga horizontal atuando sobre toda a torre (que ele representou como as forças P', P'', etc). O truque agora é pegar um corte horizontal arbitrário (que ele chamou de MN), e tirar toda a parte de cima da torre. As forças na estrutura não mudam se substituirmos a parte de cima por um ponto material de mesma massa, localizado no centróide do que foi tirado.

   diagrama 

Vamos então ver qual é o efeito do vento. Eu tenho uma força P horizontal, que é a resultante do vento. Tenho ainda três forças F1, F2, F3, que são as forças atuando sobre as vigas da treliça. Quando uma treliça está em equilíbrio, as forças sobre ela são sempre de tração ou compressão, ou seja, estão na mesma direção de suas vigas respectivas.

A torre está em equilíbrio estático, logo a soma das forças é zero. Além disso, ela também não está rotacionando; logo a soma dos torques também é zero. O torque é o produto da força pela distância; mas não é qualquer distância, é a distância na direção perpendicular à força.

Agora o truque do Eiffel é claro. Ele construiu as laterais da torre de modo que as tangentes às vigas mais externas sempre passam pelo centróide. Quando você faz isso, a direção das forças P, F1 e F3 apontam para o centróide, então o torque delas é zero. Só a força F2 está aplicando torque no centróide!

diagrama2

Podemos calcular a soma dos torques agora. A distância associada à força F2 é a distância d, em verde na figura. Note que ela é estritamente positiva.
$$\begin{align*} P\times 0 +F_1\times 0+F_2\times d+F_3\times 0&=0\\ F_2\times d&=0\\ F_2&=0 \end{align*}$$
Ahá! Se as tangentes externas apontam para o centróide, então a força na viga interna é zero! Se o vento horizontal fosse a única força atuando sobre a torre, então nós poderíamos até tirar a viga interna que não faria diferença!

É claro que na prática não podemos tirar a viga. Ela não ajuda a resistir ao vento, mas ajuda a sustentar o peso. Por outro lado, essa viga pode ser mais fraca que as outras. Com a idéia do Eiffel, só as vigas externas precisam ser fortes. Nessa foto que eu tirei do pilar oeste isso é bem visível: as vigas externas são largas e feitas de material rígido, enquanto que as internas são finas e feitas de mini-treliças, o que deixa a estrutura como um todo mais leve. E quanto mais leve, mais alto você pode subir.

ouest

A equação da Torre Eiffel


Agora que nós sabemos o truque do Eiffel, podemos calcular qual é o lugar geométrico dos pontos cujas tangentes passam pelo centróide, isso vai nos dar a equação que descreve o fomato da torre. Aparentemente o Eiffel nunca precisou fazer isso, já que ele usou "um método geométrico baseado em observações empíricas" (para mim isso é um jeito de dizer que ele calculou uma solução numérica ao invés de resolver a equação de fato).

Mas ele poderia ter resolvido a equação se quisesse. Vamos fazer as contas na caixa azul para conferir, pule se você tem medo de equações diferenciais:

Começaremos supondo que o formato da torre pode ser descrito por uma função f(x). Para simplificar as contas (e, acredite, simplifica mesmo), nós vamos rotacionar o gráfico de modo a colocar o eixo x na vertical, orientado para baixo, e o f(x) na horizontal. Vamos ainda supor que f(x) é real, contínua e pode ser diferenciada muitas vezes. grapheiffel
No gráfico acima, f(x) está desenhada em azul. Por um ponto qualquer x eu traço as retas tangentes (em vermelho), elas vão se encontrar em um ponto xC, que nós estamos supondo que sempre será o centróide da parte superior da torre (ou seja, do trecho que vai do ponto escolhido x, até o topo da torre xT). A posição do centróide é dado diretamente pela definição: $$x_C=\frac{\int_{x_T}^{x}t f(t)\,dt}{\int_{x_T}^{x} f(t)\,dt}$$ A tangente podemos calcular de duas maneiras. Nós sabemos que a tangente que passa por x é dada diretamente por f'(x). Por outro lado, ela também é a tangente do ângulo formado pelo eixo x e pela reta vermelha, que você calcula direto como cateto oposto pelo adjacente: $$f'(x)=\frac{\Delta y}{\Delta x}=\frac{f(x) -0}{x-x_C}=\frac{f(x)}{x-x_C}$$ Podemos isolar xC nessa última equação: $$\begin{align*} x-x_C&=\frac{f(x)}{f'(x)}\\ x_C&=\frac{x f'(x)-f(x)}{f'(x)} \end{align*}$$ E depois igualar as duas equações: $$\begin{align*} \frac{x f'(x)-f(x)}{f'(x)} &= \frac{\int_{x_T}^{x}t f(t)\, dt}{\int_{x_T}^{x} f(t)\, dt}\\ \left(x f'(x)-f(x)\right)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ \left(f'(x)\int_{x_T}^{x}x f(t)\, dt \right)-\left(f(x)\int_{x_T}^{x} f(t)\, dt\right) &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ f(x)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Na última linha, veja que o x passou para dentro da integral, você pode fazer isso porque ele é constante em relação à variável de integração, que é t. Antes de prosseguir, vamos dar nomes às integrais, chamando-as de g(x) e h(x). Note que essas integrais são funções de x. $$\begin{align*} g(x)&=\int_{x_T}^{x} f(t)\, dt \\ h(x)&= \int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Vamos substituir as integrais pelas funções que nomeamos, e depois derivar dos dois lados. Quando você deriva dos dois lados, você está comendo uma constante, então no final podem aparecer soluções espúrias. Depois de acabar tudo, temos que pegar as soluções e voltar para o começo, conferindo se elas são realmente válidas. $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x) \end{align*}$$ Como calcular g'(x) e h'(x)? Vamos começar por g'(x). Note, primeiro, que você pode usar o teorema fundamental do cálculo para escrever g(x) em função da primitiva de f(x), que iremos chamar de F(x): $$g(x)=\int_{x_T}^x f(t)\,dt=F(x) -F(x_T) $$ Já que F(xT) é uma constante, sua derivada é zero. F(x) é uma primitiva, e a derivada da primitiva de uma função é a própria função (desde que ela seja contínua): $$\begin{align*} g(x)&=F(x) -F(x_T)\\ g'(x)&=F'(x)-0=f(x) \end{align*} $$ Então g'(x)=f(x), essa foi fácil. Para calcular h'(x) é mais complicado. Lembremos a regra de integração por partes para integrais definidas: $$\begin{align*} \int_c^d a(t)b'(t)\,dt&=[a(t)b(t)]_c^d-\int_c^d a'(t)b(t)\,dt \end{align*} $$ Vamos escolher quem são a(t) e b(t): $$\begin{align*} a(t)=x-t&\implies a'(t)=-1\\ b(t)=F(t)&\implies b'(t)=f(t)\\ \end{align*} $$ Agora podemos quebrar h(x) em duas partes: $$\begin{align*} h(x)&=\int_{x_T}^x(x-t)f(t)\,dt\\ &=[(x-t)F(t)]_{x_T}^x-\int_{x_T}^x-F(t)\,dt\\ &=(x-x)F(x)-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ \end{align*} $$ No último passo, por falta de notação melhor, eu chamei de P(x) a primitiva de F(x). Agora podemos derivar h(x). Note que, se P(x) é a primitiva de F(x), então F(x) é a derivada de P(x). Além disso, como P(xT) é constante, sua derivada é zero: $$\begin{align*} h(x)&=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ h'(x)&=-F(x_T)+F(x)-0\\ h'(x)&=\int_{x_T}^xf(t)\,dt\\ h'(x)&=g(x)\\ \end{align*} $$ Concluímos então que h'(x) é o mesmo g(x) que tínhamos visto antes! Antes de continuar, vale a pena resumir o que nós temos até agora: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x)\\ f(x)&=g'(x)\\ f'(x)&=g''(x)\\ f''(x)&=g'''(x)\\ h'(x)&=g(x) \end{align*}$$ Note que dá para reescrever quase tudo usando só g(x) e suas derivadas. Fica faltando h(x), mas isso não é problema porque temos duas equações usando h(x), então é só isolar. Da primeira equação temos: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ h(x)&=\frac{f(x)g(x)}{f'(x)}\\ h(x)&=\frac{g'(x)g(x)}{g''(x)}\\ \end{align*}$$ Da segunda equação, temos: $$\begin{align*} f'(x)g(x)+f(x)g'(x)&=f''(x)h(x)+f'(x)h'(x)\\ f''(x)h(x)&=f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)\\ h(x)&=\frac{f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)}{f''(x)}\\ h(x)&=\frac{g''(x)g(x)+g'(x)g'(x)-g''(x)g(x)}{g'''(x)}\\ h(x)&=\frac{g'(x)^2}{g'''(x)}\\ \end{align*}$$ Igualando as duas expressões para h(x): $$\begin{align*} \frac{g'(x)g(x)}{g''(x)} &= \frac{g'(x)^2}{g'''(x)} \\ \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)} \end{align*}$$ Tudo que falta agora é resolver a equação para g(x) e finalmente temos o formato da torre. Mas essa equação é uma equação diferencial não-linear de terceira ordem. E isso é os que os matemáticos usam para assustar os filhos: "se você não comer, vai aparecer o bicho-papão, o homem do saco e a equação diferencial não-linear de terceira ordem para te pegar!" Felizmente, essa equação pode ser resolvida se você tiver a visão além do alcance. O que acontece quando você tenta diferenciar log g(x)? $$\begin{align*} \frac{d}{dx}\log(g(x))=\frac{1}{g(x)}\times g'(x)=\frac{g'(x)}{g(x)} \end{align*}$$ Ahá! Essa é exatamente a forma da equação. Fazendo o análogo no outro lado, podemos reescrever a equação original: $$\begin{align*} \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)}\\ \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x) \end{align*}$$ Agora nós podemos integrar dos dois lados (lembrando que vai aparecer uma constante), e depois tirar a exponencial dos dois lados: $$\begin{align*} \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x)\\ \log g(x) &= \log g''(x) + C\\ g(x) &= k g''(x) \\ \end{align*}$$ Agora ficou fácil né? A equação medonha virou uma equação diferencial linear de segunda ordem, essa é bem-comportada, dá a patinha e finge de morta. A solução geral é da forma abaixo: $$g(x)=A e^{B x}+C e^{-B x}$$ Como f(x)=g'(x), então f(x) também tem a forma acima, só os coeficientes que serão diferentes. A natureza da função depende dos coeficientes; por exemplo, se A=0.5, C=0.5 e B=i, então f(x)=cos(x). Mas nem todas as combinações de coeficientes são soluções válidas! Para determinar os coeficientes corretos, precisamos voltar o f(x) lá na equação original: $$f(x)\int_{x_T}^{x} f(t)\, dt = f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt $$ Daqui em diante é só fazer as contas. Se você substituir tudo, vai notar que o único jeito de conseguir uma solução real válida para qualquer x é fazendo C=0 e xT=-∞; nesse caso A pode ser um real qualquer, e B pode ser um real positivo qualquer. Estritamente falando, a torre teria que ser infinitamente alta para satisfazer a essa equação, mas na prática você pode aproximar "infinitamente alta" por "muito, muito alta". Concluindo, a solução final para o formato da torre Eiffel é: $$f(x)=A e^{B x} $$ Ou seja, a torre tem um perfil exponencial!

Matemática e Engenharia


Agora que sabemos qual o formato teórico da torre, podemos usar os dados reais da torre para estimar os coeficientes da exponencial. Eu rodei um best fit e cheguei em A=3.34323 e B=0.0102592:

bestfit

A aproximação é boa, mas... parece que não encaixa direitinho? É meio frustrante fazer esse monte de contas, e no fim a curva real não ser tão parecida assim com a equação deduzida.

Mas tem um motivo para isso, e o motivo é que o Eiffel era engenheiro, não matemático. Quando você plota uma exponencial em um gráfico semi-log, o resultado tem que ser uma linha. Mas olha só o que acontece com os dados reais da torre: eles formam duas linhas!

semilog 

Coisas que você aprende na sua primeira aula de engenharia: se você faz o projeto de um elevador, e determina que o cabo do elevador precisa aguentar no mínimo 10 toneladas, qual o cabo que você coloca no elevador? Um matemático diria que é um cabo de 10 toneladas, mas o engenheiro vai falar que o mínimo são 14 toneladas.

Isso é o resultado do fator de segurança. Por mais que o seu modelo matemático seja correto, sempre tem alguma coisa que você não levou em conta (pode ser que o aço que você usou tenha impurezas, ou que o solo não é tão firme quanto você achava, e assim por diante). Por isso, é costume sempre multiplicar o valor final por um fator de segurança para levar em contas esses imprevistos.

O Eiffel sabia disso, e no diário ele explica que imaginava que a base da torre estaria sujeita a torques maiores que o topo, e por isso usou um fator de segurança diferente na base e no topo. Quando você faz o best fit com duas exponenciais ao invés de uma só, o gráfico fica bem melhor!

twofit 

Agora finalmente podemos dar o mistério por resolvido: o formato da Torre Eiffel é descrito por duas exponenciais, escolhidas para minimizar o material necessário para resistir à carga do vento.

A Torre Eiffel é o literalmente o maior monumento à ciência construído, da próxima vez que passar por Paris aproveite para fazer uma reverência aos cientistas do passado :)

eiffeltower

segunda-feira, 31 de março de 2014

Matemática x Tanques de Guerra

Quando eu era criança pequena lá na década de 80, eu tinha um brinquedo de tanque de guerra. Ele chamava Panzer, e era fabricado pela brinquedos Mimo. Para a época era bem divertido: tinha controle remoto, e, além de andar, ainda disparava mísseis!

panzer

Mas eu não percebia a ironia do brinquedo. O Panzer não era só um brinquedo, ele existiu de verdade. Na vida real, o Panzer era uma máquina nazista de destruição. A série Panzer de tanques alemães foi responsável pela morte direta de incontáveis aliados. Ver um Panzer surgindo no horizonte era como ver a morte vindo em sua direção. Era um mimo esse Panzer.

Por outro lado, o Panzer era forte, mas não era invencível. Os aliados tinham armas capazes de deter o Panzer, o problema era descobrir quantas dessas armas eles deveriam levar para o front. Poucas semanas antes do Dia D, os alemães começaram a usar o novo modelo Panzer V, e os aliados não tinham como estimar quantas armas levar sem saber quantos Panzer V os inimigos fabricaram.

Foi então que os aliados usaram a arma definitiva contra o Panzer: a Matemática! Usando um truque muito esperto, eles estimaram que os alemães estavam fabricando 270 Panzers por mês. Depois que a guerra acabou, eles foram nos arquivos alemães conferir os números, e na verdade estavam fazendo 276 Panzers por mês. A estimativa foi muito boa!


A Quantidade de Tanques


O truque dos aliados se baseava em um erro dos alemães e em uma hipótese estatística. O erro dos alemães era que os tanques tinham número de série sequencial. Quando saíam da fabrica, eram numerados como 1, 2, 3, e assim por diante. Ao fazer isso, eles deixaram uma brecha para os hackers aliados.

Suponha que você capturou um Panzer alemão, e o serial dele era 14. Você consegue estimar quantos tanques eles produziram baseado só nesse número? No caso geral não, mas você pode fazer uma hipótese simplificadora, que é supor que o tanque que você capturou é aleatório (ou seja, os alemães estavam distribuindo igualmente os tanques por toda a Europa).

Com essa hipótese a intuição é clara. Vamos supor, hipoteticamente, que os alemães fizeram 200 tanques. Qual chance de você capturar um tanque de serial menor ou igual a 14 nessa situação? É 14/200 né? E se os alemães fizeram só 20 tanques? Aí a chance é 14/20, dez vezes maior que 14/200. Com a hipótese de tanque aleatório, é mais provável que os alemães tenham feito 20 tanques do que 200 tanques.

Você pode formalizar esse argumento usando o teorema de Bayes, e calcular o valor esperado do número de tanques. Vamos fazer as contas na caixa azul:

Suponha que conseguimos capturar k tanques. Dentre esses k tanques, o maior serial encontrado tem o valor m. Com essas informações, vamos calcular qual é o valor esperado do número de tanques produzidos, que vamos chamar de n. Podemos abrir esse valor esperado pela definição:
$$E[n\;|\;m,k] = \sum_{n} n \; p(n\; |\; m,k)$$
O problema então é achar a probabilidade condicional de n, dados m e k. Nós ainda não sabemos quanto vale essa probabilidade, mas podemos abri-la usando o teorema de Bayes:
$$p(n\;|\;m,k)=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)}$$
Essas três probabilidades nós conseguimos calcular! A primeira delas, p(m,k|n), sai por combinatória. Dado um vetor de tamanho n, de quantas maneiras conseguimos escolher k elementos, de modo que o maior deles seja m?

Um dos elementos precisa necessariamente ser o m, então só precisamos calcular a posição dos outros k-1 restantes. E eles precisam ser menores que m, então só podem estar nas m-1 posições menores que m. Por isso, essa quantidade vale binomial(m-1k-1). Note, entretanto, que isso só vale quando n>=m, caso contrário a probabilidade é zero (não tem como o máximo ser m se o vetor tem um tamanho menor que m).

Para calcular a probabilidade desejada, temos que achar de quantas maneiras podemos escolher k elementos dentre n, sem considerar a restrição do maior deles ser m. Mas isso é o próprio binomial(n, k). Então, a probabilidade é:
$$p(m,k\;|\;n) = \frac{\displaystyle{m-1\choose k-1}} {\displaystyle{n\choose k}}[n\ge m]$$
Vamos agora para a segunda, p(n). Eu não tenho nenhuma informação sobre quantos tanques foram produzidos, então não sei qual é a probabilidade de ter n tanques. Mas eu posso chutar que a distribuição é uniforme. A chance de ter 15 tanques, ou de ter 500 tanques, a priori, deve ser a mesma. Por isso, vamos fazer p(n) ser uma constante independente de n:
$$p(n)=c$$
A última, p(m,k), é a probabilidade do máximo ser m em k escolhas, independente do valor de n. Para isso, podemos usar a lei da probabilidade total
$$p(m,k) = \sum_n p(m,k\;|\;n) \;p(n)$$ 
Essas duas já sabemos calcular! É só substituir:
 $$\begin{align*} p(m,k) &= \sum_n p(m,k\;|\;n) \;p(n) \\ &= \sum_n \frac{{m-1\choose k-1}} {{n\choose k}}[n\ge m]\;c\\ &= \sum_{n\ge m} c {m-1\choose k-1} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ \end{align*} $$
A somatória tem uma forma fechada, que você pode achar com o algoritmo de Gosper:
$$\sum_{n=m}^{\infty}{n\choose k}^{-1}=\frac{m}{k-1}{m\choose k}^{-1}$$ 
Substituindo:
$$\begin{align*} p(m,k) &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \times \frac{m}{k-1}{m\choose k}^{-1} \\ &= \frac{cm}{k-1} {m-1 \choose k-1} {m\choose k}^{-1} \\ \end{align*} $$
Agora podemos voltar e substituir as três probalidades na fórmula de Bayes:
$$ \begin{align*} p(n\;|\;m,k)&=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)} \\ &= \frac{\displaystyle{m-1\choose k-1}{n\choose k}^{-1}[n\ge m]\times c} {\displaystyle\frac{cm}{k-1}{m-1\choose k-1}{m\choose k}^{-1}}\\ &= \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m] \end{align*}$$
Por fim, podemos substituir na fórmula do valor esperado:
 $$\begin{align*} E[n\;|\;m,k] &= \sum_{n} n \; p(n\; |\; m,k)\\ &= \sum_n n \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m]\\ &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ \end{align*}$$
Essa somatória também tem uma forma fechada por Gosper:
$$ \sum_{n=m}^{\infty} n{n\choose k}^{-1}=\frac{m(m-1)}{k-2}{m\choose k}^{-1}$$
Chegamos então na substituição final:
$$\begin{align*} E[n\;|\;m,k] &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ &= \frac{k-1}{m}{m\choose k}\times\frac{m(m-1)}{k-2}{m\choose k}^{-1} \\ &= \frac{(k-1)(m-1)}{k-2} \\ \end{align*}$$

A Estimativa na Prática


A fórmula final é super simples, mas ela funciona mesmo na prática? Podemos testá-la fazendo uma simulação numérica. Para isso, eu fiz dez mil simulações. Em cada uma eu sorteio o número de tanques fabricados, e escolho aleatoriamente 3 deles. Aí eu calculo a estimativa pela fórmula, e normalizo o erro. O código está no github, e o histograma resultante é o abaixo. A normalização é tal que o valor 0 é uma estimativa totalmente correta:

germantank3

A estimativa não ficou muito boa não! Ao invés de ter uma gaussiana em torno do zero, o estimador tem um bias muito forte para os negativos. Felizmente, esse problema só acontece porque capturamos apenas 3 tanques. Se, ao invés disso, tivéssemos capturado 30 tanques, então o histograma seria bem melhor. Compare o histograma de 30 tanques sobreposto ao de 3 tanques:

germantank30

Bem melhor né? Agora quase todas as estimativas estão bem próximas do zero. Esse método é bastante sensível ao número de tanques capturados. Os aliados sabiam disso, por isso que a estimativa deles foi tão boa (estimaram 270 tanques/mês, quando o valor real era 276). Sabe quantos tanques eles capturaram para fazer essa estimativa? Só dois tanques!

Claro que tem um truque. A estimativa com dois tanques é muito ruim, mas eles notaram que os alemães colocavam serial em todas as peças do tanque. Em especial, cada tanque tinha 48 rodas, cada uma com um serial único. Por isso, eles conseguiram usar a fórmula com k=96, o que deu a precisão alta que queriam. E no fim a matemática ganhou a guerra :)

sábado, 1 de março de 2014

O Objetivo Final da Matemática

Quem acompanhou as contas na caixa azul do post sobre o Pokèmon Twitch deve ter notado que eu pulei um pedaço. Em certo ponto da demonstração, eu magicamente tirei do chapéu uma fórmula nada óbvia:

presto (1)

E como foi que eu cheguei nessa fórmula? Simples, eu calculei no Wolfram Alpha:

wolframalpha

"Ricbit, você trapaceou! Que roubalheira! Um matemático de verdade jamais faria isso!"

Mas será que é trapaça mesmo? Nesse post, eu vou mostrar como se resolve essa conta sem usar o Wolfram Alpha; e, se você acha que é trapaça, talvez acabe mudando de idéia :)

A progressão hipergeométrica


Para resolver essa somatória sem usar o Wolfram Alpha, nós precisamos relembrar as progressões que aprendemos no colégio. A primeira delas foi a progressão aritmética, onde a diferença entre dois números consecutivos é constante:

   

A segunda delas é a progressão geométrica, onde a razão entre dois números consecutivos é constante:

   

Uma progressão extremamente útil que não ensinam no colégio é a progressão hipergeométrica. Ela parece com a progressão geométrica, porque é definida a partir da razão entre dois elementos consecutivos. Mas agora, ao invés dessa razão ser constante, nós permitimos que a razão seja uma função racional (que é só um jeito chique de dizer que ela é o quociente entre dois polinômios):

 

Naquela somatória que queremos resolver, os termos sendo somados formam uma progressão hipergeométrica. Confira:

   

O numerador é um polinômio em q, o denominador é um polinômio em q. Então confere com a nossa definição de progressão hipergeométrica.

A soma da progressão hipergeométrica


Certamente vocês lembram que as progressões do colégio tinham uma fórmula para a soma de seus termos. A soma dos n primeiros termos da progressão aritmética é:


A soma dos n primeiros termos da progressão geométrica é:

 

A soma dos n primeiros termos da progressão hipergeométrica é... oops! Não existe uma fórmula geral para a soma de termos hipergeométricos.

Antigamente, achar uma fórmula para uma soma hipergeométrica era tão raro que elas até ganhavam nomes especiais. Por exemplo, a igualdade abaixo tem um nome, é o Binômio de Newton:

 

Outras somas hipergeométricas não tem fórmula, mas são tão úteis, que a própria soma ganha um nome e passa a valer como se fosse uma fórmula:


Por muitos séculos esse foi o estado da arte. Alguém acha uma fórmula aqui, alguém nomeia uma soma ali. Mas isso mudou no meio do século 20.

O algoritmo de Gosper


O primeiro progresso na direção de um método geral de resolução de somas hipergeométricas foi feito pela Mary Celine Fasenmyer em 1945. Depois de se formar em matemática, ela se juntou à ordem religiosa Sisters of Mercy (sim, a mesma que inspirou o nome da banda), e lá ela criou o Método da Irmã Celine, um algoritmo que consegue resolver automaticamente alguns tipos de somas hipergeométicas.

Em 1978 o método foi aprimorado por Bill Gosper e se tornou o algoritmo de Gosper (muito embora eu desconfie que, se a Irmã Celine fosse homem, o algoritmo seria de Gosper-Celine). O método de Gosper se baseia em um teorema surpreendente: A solução de uma soma hipergeométrica é outro hipergeométrico, se e somente se ela puder ser escrita na forma abaixo, onde R[j] é uma função racional.


Essa somatória parece esquisita. Se j é o índice da somatória, então o que ele está fazendo no lado direito da equação? E cadê os limites da somatória?

A resposta é que essa é uma somatória indefinida. Funciona do mesmo jeito que no cálculo: lá você tem integrais definidas e integrais indefinidas, aqui nós temos somatórias definidas (com limites), e somatórias indefinidas (sem limites). Para uma somatória indefinida, sempre vale a propriedade abaixo:

 

Essas duas equações são tudo que nós precisamos para usar o algoritmo de Gosper. Vamos à caixa azul calcular a nossa somatória original:

Começamos substituindo uma equação na outra:


Podemos então dividir tudo por h[j]:

 

Opa, a razão entre os h nós já calculamos lá em cima para a nossa somatória!

 

Falta achar R[j]. Eu sei que ele é um quociente de polinômios, mas não sei qual o grau desses polinômios. Eu vou chutar que o grau é 1, se der errado depois eu volto e aumento o grau. Pois bem, se o grau é 1, então R[j] é da forma abaixo:

 

Agora eu vou substituir. Aperte o cinto que lá vem turbulência:

   

O resultado é um polinômio gigante em j que é identicamente nulo. Se ele é identicamente nulo, então cada um dos seus coeficientes precisa ser zero. Com isso podemos montar um sistema de equações:

 

Agora é "só" resolver o sistema! Vou poupá-los de páginas e páginas, o resultado é o seguinte:

   

Se o sistema não tivesse solução, eu teria que voltar lá em cima e aumentar o grau do polinômio (se as contas foram titânicas com grau 1, imagine com grau 2). Mas o processo não é infinito, existem métodos que permitem achar um limite superior para o grau do polinômio. Se você atingir esse limite sem achar nenhuma solução, então você provou que sua somatória não pode ser escrita como um hipergeométrico.

Voltando à nossa equação original, agora nós sabemos quanto vale a somatória indefinida:

   

Nós queremos a soma de 1 até m, então precisamos converter a soma indefinida em uma soma definida, usando o teorema fundamental do cálculo discreto:

 

Aplicando o teorema, temos:

   

Pronto, esse é o resultado final. Eu acho muito importante que a pessoa resolva o algoritmo de Gosper uma vez na vida, para nunca mais cometer a sandice de repetir o processo. O método é mecânico, você só precisa seguir as regras e não errar as contas. O computador faz isso melhor que você. Use o Wolfram Alpha e economize grafite.

O objetivo final da Matemática


No final do século 19, o Bertrand Russell e o Alfred Whitehead tiveram uma idéia. Ele bolaram um plano para axiomatizar toda a matemática, de uma maneira completa e livre de contradições. Se o plano funcionasse, então toda a Matemática seria reduzida à manipulação de símbolos, que é um processo puramente mecânico e automatizável. O plano começou muito bem, a ponto do Whitehead fazer a seguinte declaração:

"A Civilização avança estendendo o número de operações importantes que podemos fazer sem ter que pensar."

Baseado nessa citação, o Concrete Mathematics do Knuth vai além e conclui:

"O objetivo final da matemática é eliminar a necessidade de todo pensamento inteligente."

É verdade, o objetivo sempre foi esse. Para os babilônios, resolver uma equação de segundo grau era dificílimo, só os grandes mestres conseguiam. Hoje em dia, nós transformamos a equação de segundo grau em uma fórmula que é ensinada para crianças de 14 anos. Você não precisa ser inteligente para usar a fórmula, é só substituir os números e fazer as contas.

Imagine se fosse possível fazer isso com toda a matemática, ao invés de só com equações de segundo grau! Você poderia usar computadores para fazer as contas, no lugar de crianças de 14 anos. A Matemática seria um problema resolvido, e aí você pode se concentrar no que interessa, que são as aplicações. Um físico acabou de descobrir uma propriedade quântica nova, será que ela permite criar um aparelho de teletransporte? Eu monto a equação, jogo no computador, e tcharam! Está aqui a resposta!

Toda a saga por trás do plano do Russell e do Whitehead é contada em quadrinhos no fabuloso Logicomix. Infelizmente, o final da história é que o plano falhou. E falhou espetacularmente. Em 1931, o Gödel provou que nenhum sistema axiomático é capaz de provar todas as sentenças verdadeiras que é capaz de descrever, acabando de uma vez por todas com o sonhado plano da axiomatização completa.

E o século 20 desceu ladeira abaixo. O Church provou que, mesmo se a axiomatização completa tivesse tido sucesso, seria impossível criar um método que prove um teorema qualquer a partir dos axiomas. O Turing provou que é impossível criar um método que analise um programa de computador e diga se ele termina ou se entra em loop infinito. O Matiyasevich provou que é impossível criar um método que resolva todas as equações diofantinas. O Richardson provou que é impossível criar um método que prove todas as identidades trigonométricas. Foi o fim da matemática automatizada... ou não?

Claro que não! Todos esses teoremas provam a impossibilidade do caso geral. Mas você sempre pode automatizar um caso específico. Não existe uma fórmula geral para a soma hipergeométrica, mas o algoritmo de Gosper prova que pelo menos um caso especial pode ser automatizado: o caso da soma hipergeométrica cujo resultado também é um termo hipergeométrico. E esse caso é imensamente útil na prática, especialmente se você estuda computação ou combinatória.

E no fim das contas, usar o Wolfram Alpha para resolver uma somatória é trapacear ou não? Por dentro, tudo que o Wolfram Alpha faz é usar o algoritmo de Gosper. Se você acredita que usar um método computacional é roubar, então por consistência precisa abrir de mão de todos os outros métodos computacionais. Não pode mais resolver a equação de segundo grau usando a fórmula de Bhaskara, não pode mais calcular o máximo divisor comum com o algoritmo de Euclides, não pode nem calcular quanto é 463 vezes 367 sem antes decorar a tabuada do 367.

Para mim isso não faz sentido. E por isso eu uso o Wolfram Alpha sem medo de ser feliz :)